IMP  2.3.0
The Integrative Modeling Platform
VectorD.h
Go to the documentation of this file.
1 /**
2  * \file IMP/algebra/VectorD.h \brief Simple D vector class.
3  *
4  * Copyright 2007-2014 IMP Inventors. All rights reserved.
5  *
6  */
7 
8 #ifndef IMPALGEBRA_VECTOR_D_H
9 #define IMPALGEBRA_VECTOR_D_H
10 
11 #include <IMP/algebra/algebra_config.h>
12 #include "VectorBaseD.h"
13 #include <IMP/base/types.h>
14 #include <IMP/base/check_macros.h>
15 #include <IMP/base/exception.h>
16 #include <IMP/base/utility.h>
17 #include <IMP/base/InputAdaptor.h>
18 #include <IMP/base/random.h>
19 #include "algebra_macros.h"
20 #include <boost/random/variate_generator.hpp>
21 #include <boost/random/normal_distribution.hpp>
22 #include <boost/range.hpp>
23 #include "internal/vector.h"
24 
25 #include <limits>
26 #include <cmath>
27 #include <boost/random/normal_distribution.hpp>
28 #include <boost/static_assert.hpp>
29 
30 #if IMP_HAS_CHECKS >= IMP_USAGE
31 #define IMP_ALGEBRA_VECTOR_CHECK check_vector()
32 #define IMP_ALGEBRA_VECTOR_CHECK_INDEX(i) check_index(i)
33 #define IMP_ALGEBRA_VECTOR_CHECK_COMPATIBLE(o) \
34  check_compatible_vector(o); \
35  o.check_vector()
36 #else
37 #define IMP_ALGEBRA_VECTOR_CHECK
38 #define IMP_ALGEBRA_VECTOR_CHECK_INDEX(i)
39 #define IMP_ALGEBRA_VECTOR_CHECK_COMPATIBLE(o)
40 #endif
41 
42 IMPALGEBRA_BEGIN_NAMESPACE
43 //! A Cartesian vector in D-dimensions.
44 /** Store a vector of Cartesian coordinates. It supports all expected
45  mathematical operators, including using * for the dot product.
46  \see Vector3D
47  \see Vector2D
48 
49  \geometry
50  */
51 template <int D>
52 class VectorD : public VectorBaseD<D> {
53  /* implementing it via a specialization is in order to get swig
54  to only use the right constructors as well as C++. */
55  IMP_ALGEBRA_VECTOR_METHODS(D);
56 };
57 
58 template <>
59 class VectorD<-1> : public VectorBaseD<-1> {
60  std::vector<double> get_range(double x0, double x1, double x2, double x3,
61  double x4, double x5) {
62  IMP_USAGE_CHECK(x0 != std::numeric_limits<double>::max(),
63  "Bad init values");
64  std::vector<double> ret;
65  ret.push_back(x0);
66  if (x1 != std::numeric_limits<double>::max()) ret.push_back(x1);
67  if (x2 != std::numeric_limits<double>::max()) ret.push_back(x2);
68  if (x3 != std::numeric_limits<double>::max()) ret.push_back(x3);
69  if (x4 != std::numeric_limits<double>::max()) ret.push_back(x4);
70  if (x5 != std::numeric_limits<double>::max()) ret.push_back(x5);
71  return ret;
72  }
73 
74  public:
75  IMP_ALGEBRA_VECTOR_METHODS(-1);
76  explicit VectorD(double x0, double x1 = std::numeric_limits<double>::max(),
77  double x2 = std::numeric_limits<double>::max(),
78  double x3 = std::numeric_limits<double>::max(),
79  double x4 = std::numeric_limits<double>::max(),
80  double x5 = std::numeric_limits<double>::max())
81  : P(get_range(x0, x1, x2, x3, x4, x5)) {}
82 };
83 
84 template <>
85 class VectorD<1> : public VectorBaseD<1> {
86  public:
87  IMP_ALGEBRA_VECTOR_METHODS(1);
88 
89  //! Initialize the 1-vector from its value.
90  VectorD(double x) { P::operator[](0) = x; }
91 };
92 template <>
93 class VectorD<2> : public VectorBaseD<2> {
94  public:
95  IMP_ALGEBRA_VECTOR_METHODS(2);
96 
97  VectorD(double x, double y) {
98  P::operator[](0) = x;
99  P::operator[](1) = y;
100  }
101 };
102 template <>
103 class VectorD<3> : public VectorBaseD<3> {
104  public:
105  IMP_ALGEBRA_VECTOR_METHODS(3);
106 
107  VectorD(double x, double y, double z) {
108  P::operator[](0) = x;
109  P::operator[](1) = y;
110  P::operator[](2) = z;
111  }
112 };
113 template <>
114 class VectorD<4> : public VectorBaseD<4> {
115  public:
116  IMP_ALGEBRA_VECTOR_METHODS(4);
117 
118  VectorD(double x0, double x1, double x2, double x3) {
119  P::operator[](0) = x0;
120  P::operator[](1) = x1;
121  P::operator[](2) = x2;
122  P::operator[](3) = x3;
123  }
124 };
125 template <>
126 class VectorD<5> : public VectorBaseD<5> {
127  public:
128  IMP_ALGEBRA_VECTOR_METHODS(5);
129 
130  VectorD(double x0, double x1, double x2, double x3, double x4) {
131  P::operator[](0) = x0;
132  P::operator[](1) = x1;
133  P::operator[](2) = x2;
134  P::operator[](3) = x3;
135  P::operator[](4) = x4;
136  }
137 };
138 template <>
139 class VectorD<6> : public VectorBaseD<6> {
140  public:
141  IMP_ALGEBRA_VECTOR_METHODS(6);
142 
143  VectorD(double x0, double x1, double x2, double x3, double x4, double x5) {
144  P::operator[](0) = x0;
145  P::operator[](1) = x1;
146  P::operator[](2) = x2;
147  P::operator[](3) = x3;
148  P::operator[](4) = x4;
149  P::operator[](5) = x5;
150  }
151 };
152 
153 #ifndef IMP_DOXYGEN
154 
155 template <int D>
156 inline std::ostream &operator<<(std::ostream &out, const VectorD<D> &v) {
157  v.show(out);
158  return out;
159 }
160 
161 #ifndef SWIG
162 // SWIG 3.0.0 cannot parse operator>>
163 template <int D>
164 inline std::istream &operator>>(std::istream &in, VectorD<D> &v) {
165  for (unsigned int i = 0; i < D; ++i) {
166  in >> v[i];
167  }
168  return in;
169 }
170 #endif
171 
172 #endif
173 
174 //! lexicographic comparison of two vectors
175 /** Note that this is not very reliable and probably should not be used.
176  \see VectorD
177  */
178 template <int D>
179 inline int compare(const VectorD<D> &a, const VectorD<D> &b) {
180  IMP_USAGE_CHECK(a.get_dimension() == b.get_dimension(),
181  "Dimensions don't match.");
182  for (unsigned int i = 0; i < a.get_dimension(); ++i) {
183  if (a[i] < b[i])
184  return -1;
185  else if (a[i] > b[i])
186  return 1;
187  }
188  return 0;
189 }
190 
191 /** \see VectorD */
192 template <int D>
193 inline VectorD<D> operator*(double s, VectorD<D> o) {
194  return o *= s;
195 }
196 
197 //! Compute the squared distance between two vectors
198 /** \see VectorD
199  */
200 template <int D>
201 inline double get_squared_distance(const VectorD<D> &v1, const VectorD<D> &v2) {
202  return (v1 - v2).get_squared_magnitude();
203 }
204 
205 //! Compute the distance between two vectors
206 /** \see VectorD
207  */
208 template <int D>
209 inline double get_distance(const VectorD<D> &v1, const VectorD<D> &v2) {
210  return std::sqrt(get_squared_distance(v1, v2));
211 }
212 
213 //! Return the basis vector for the given coordinate
214 /** Return the unit vector pointing in the direction of the requested
215  coordinate. That is
216  \code
217  get_basis_vector_d<3>(2)== Vector3D(0,0,1);
218  \endcode
219  \see VectorD
220  */
221 template <int D>
222 inline VectorD<D> get_basis_vector_d(unsigned int coordinate) {
223  IMP_USAGE_CHECK(coordinate < D, "There are only " << D << " basis vectors");
224  double vs[D];
225  for (unsigned int i = 0; i < D; ++i) {
226  if (i == coordinate)
227  vs[i] = 1;
228  else
229  vs[i] = 0;
230  }
231  return VectorD<D>(vs, vs + D);
232 }
233 
234 //! Return a dynamically sized basis vector
235 inline VectorD<-1> get_basis_vector_kd(int D, unsigned int coordinate) {
236  IMP_USAGE_CHECK(D > 0, "D must be positive");
237  IMP_USAGE_CHECK(coordinate < static_cast<unsigned int>(D),
238  "There are only " << D << " basis vectors");
239  boost::scoped_array<double> vs(new double[D]);
240  for (int i = 0; i < D; ++i) {
241  if (i == static_cast<int>(coordinate))
242  vs[i] = 1;
243  else
244  vs[i] = 0;
245  }
246  return VectorD<-1>(vs.get(), vs.get() + D);
247 }
248 
249 //! Return a vector of zeros
250 template <int D>
252  IMP_USAGE_CHECK(D > 0, "D must be positive");
253  VectorD<D> ret;
254  for (int i = 0; i < D; ++i) {
255  ret[i] = 0;
256  }
257  return ret;
258 }
259 
260 //! Return a dynamically sized vector of zeros
261 template <int D>
263  IMP_USAGE_CHECK(D == Di, "D must be positive");
264  IMP_UNUSED(Di);
265  return get_zero_vector_d<D>();
266 }
267 
268 //! Return a dynamically sized vector of zeros
269 inline VectorD<-1> get_zero_vector_kd(int D) {
270  IMP_USAGE_CHECK(D > 0, "D must be positive");
271  Floats vs(D, 0);
272  return VectorD<-1>(vs.begin(), vs.end());
273 }
274 
275 //! Return a vector of ones (or another constant)
276 template <int D>
277 inline VectorD<D> get_ones_vector_d(double v = 1) {
278  IMP_USAGE_CHECK(D > 0, "D must be positive");
279  VectorD<D> ret;
280  for (unsigned int i = 0; i < D; ++i) {
281  ret[i] = v;
282  }
283  return ret;
284 }
285 
286 //! Return a vector of ones (or another constant)
287 /** Di must equal D. */
288 template <int D>
289 inline VectorD<D> get_ones_vector_kd(unsigned int Di, double v = 1) {
290  IMP_USAGE_CHECK(D == Di, "D must be equal");
291  IMP_UNUSED(Di);
292  return get_ones_vector_d<D>(v);
293 }
294 
295 //! Return a vector of ones (or another constant)
296 inline VectorD<-1> get_ones_vector_kd(unsigned int D, double v = 1) {
297  IMP_USAGE_CHECK(D > 0, "D must be positive");
298  boost::scoped_array<double> vv(new double[D]);
299  for (unsigned int i = 0; i < D; ++i) {
300  vv[i] = v;
301  }
302  return VectorD<-1>(vv.get(), vv.get() + D);
303 }
304 
305 #ifndef SWIG
306 
307 /** \name Norms
308  We define a number of standard, \f$L^p\f$, norms on VectorD.
309  - \f$L^1\f$ is the Manhattan distance, the sum of the components
310  - \f$L^2\f$ is the standard Euclidean length
311  - \f$L^{\inf}\f$ is the maximum of the components
312  @{
313 */
314 
315 template <int D>
316 inline double get_l2_norm(const VectorD<D> &v) {
317  return v.get_magnitude();
318 }
319 
320 template <int D>
321 inline double get_l1_norm(const VectorD<D> &v) {
322  double n = std::abs(v[0]);
323  for (unsigned int i = 1; i < v.get_dimension(); ++i) {
324  n += std::abs(v[i]);
325  }
326  return n;
327 }
328 
329 template <int D>
330 inline double get_linf_norm(const VectorD<D> &v) {
331  double n = std::abs(v[0]);
332  for (unsigned int i = 1; i < v.get_dimension(); ++i) {
333  n = std::max(n, std::abs(v[i]));
334  }
335  return n;
336 }
337 
338 /** @} */
339 
340 #ifndef IMP_DOXYGEN
341 
342 template <int D>
343 struct SpacesIO {
344  const VectorD<D> &v_;
345  SpacesIO(const VectorD<D> &v) : v_(v) {}
346 };
347 
348 template <int D>
349 struct CommasIO {
350  const VectorD<D> &v_;
351  CommasIO(const VectorD<D> &v) : v_(v) {}
352 };
353 template <int D>
354 inline std::ostream &operator<<(std::ostream &out, const SpacesIO<D> &s) {
355  s.v_.show(out, " ", false);
356  return out;
357 }
358 template <int D>
359 inline std::ostream &operator<<(std::ostream &out, const CommasIO<D> &s) {
360  s.v_.show(out, ", ", false);
361  return out;
362 }
363 
364 //! Use this before outputting to delimited vector entries with a space
365 /** std::cout << spaces_io(v);
366  produces "1.0 2.0 3.0"
367  \see VectorD
368  */
369 template <int D>
370 inline SpacesIO<D> spaces_io(const VectorD<D> &v) {
371  return SpacesIO<D>(v);
372 }
373 
374 //! Use this before outputting to delimited vector entries with a comma
375 /** std::cout << commas_io(v);
376  produces "1.0, 2.0, 3.0"
377  \see VectorD
378  */
379 template <int D>
380 inline CommasIO<D> commas_io(const VectorD<D> &v) {
381  return CommasIO<D>(v);
382 }
383 #endif // doxygen
384 
385 #endif // swig
386 /** 1D vector typedef for swig */
388 /** 1D vectors typedef for swig */
390 /** 2D vector typedef for swig */
392 /** 2D vectors typedef for swig */
394 /** 3D vector typedef for swig */
396 /** 3D vectors typedef for swig */
398 /** 4D vector typedef for swig */
400 /** 4D vectors typedef for swig */
402 /** 5D vector typedef for swig */
404 /** 5D vectors typedef for swig */
406 /** 6D vector typedef for swig */
408 /** 6D vector typedef for swig */
410 /** KD vector typedef for swig */
411 typedef VectorD<-1> VectorKD;
412 /** KD vectors typedef for swig */
414 
415 #ifndef SWIG
416 /** \see VectorD \genericgeometry */
417 template <class C>
418 inline const VectorD<C::DIMENSION> &get_vector_geometry(const C &g) {
419  return g;
420 }
421 /** \see VectorD \genericgeometry */
422 template <class C, class E>
423 inline void set_vector_geometry(C &g, const E &v) {
424  g = v;
425 }
426 #endif
427 
428 //! Return the vector that is the elementwise product of the two.
429 /** \see VectorD */
430 template <int D>
432  const algebra::VectorD<D> &b) {
433  VectorD<D> ret(a);
434  for (unsigned int i = 0; i < ret.get_dimension(); ++i) {
435  ret[i] *= b[i];
436  }
437  return ret;
438 }
439 
440 //! Return the vector that is the elementwise product of the two.
441 /** \see VectorD */
442 template <int D>
444  const algebra::VectorD<D> &b) {
445  IMP_USAGE_CHECK(a.size() == b.get_dimension(), "Dimensions don't match,");
446  VectorD<D> ret(b);
447  for (unsigned int i = 0; i < ret.get_dimension(); ++i) {
448  ret[i] *= a[i];
449  }
450  return ret;
451 }
452 
453 IMPALGEBRA_END_NAMESPACE
454 
455 #endif /* IMPALGEBRA_VECTOR_D_H */
Basic types used by IMP.
VectorD< 6 > Vector6D
Definition: VectorD.h:407
VectorD< 2 > Vector2D
Definition: VectorD.h:391
VectorD< D > get_zero_vector_kd(int Di)
Return a dynamically sized vector of zeros.
Definition: VectorD.h:262
VectorD< 5 > Vector5D
Definition: VectorD.h:403
base::Vector< VectorD< 1 > > Vector1Ds
Definition: VectorD.h:389
base::Vector< VectorD< 2 > > Vector2Ds
Definition: VectorD.h:393
double get_squared_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the squared distance between two vectors.
Definition: VectorD.h:201
base::Vector< VectorD< 5 > > Vector5Ds
Definition: VectorD.h:405
Exception definitions and assertions.
VectorD< 4 > Vector4D
Definition: VectorD.h:399
VectorD< D > get_ones_vector_kd(unsigned int Di, double v=1)
Return a vector of ones (or another constant)
Definition: VectorD.h:289
VectorD<-1 > get_basis_vector_kd(int D, unsigned int coordinate)
Return a dynamically sized basis vector.
Definition: VectorD.h:235
const VectorD< C::DIMENSION > & get_vector_geometry(const C &g)
Definition: VectorD.h:418
base::Vector< VectorD< 3 > > Vector3Ds
Definition: VectorD.h:397
Basic types used by IMP.
base::Vector< VectorD< 4 > > Vector4Ds
Definition: VectorD.h:401
A Cartesian vector in D-dimensions.
Definition: VectorD.h:52
#define IMP_UNUSED(variable)
VectorD< D > get_elementwise_product(const Ints &a, const algebra::VectorD< D > &b)
Return the vector that is the elementwise product of the two.
Definition: VectorD.h:443
int compare(const VectorD< D > &a, const VectorD< D > &b)
lexicographic comparison of two vectors
Definition: VectorD.h:179
base::Vector< VectorD< 6 > > Vector6Ds
Definition: VectorD.h:409
VectorD< D > get_zero_vector_d()
Return a vector of zeros.
Definition: VectorD.h:251
VectorD< D > get_basis_vector_d(unsigned int coordinate)
Return the basis vector for the given coordinate.
Definition: VectorD.h:222
Various general useful functions for IMP.
VectorD< 1 > Vector1D
Definition: VectorD.h:387
void set_vector_geometry(C &g, const E &v)
Definition: VectorD.h:423
base::Vector< VectorD<-1 > > VectorKDs
Definition: VectorD.h:413
Simple D vector class.
VectorD< D > get_ones_vector_d(double v=1)
Return a vector of ones (or another constant)
Definition: VectorD.h:277
Exception definitions and assertions.
VectorD< 3 > Vector3D
Definition: VectorD.h:395
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:170
VectorD< D > operator*(double s, VectorD< D > o)
Definition: VectorD.h:193
double get_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the distance between two vectors.
Definition: VectorD.h:209
Random number generators used by IMP.
A Cartesian vector in D-dimensions.
Definition: VectorBaseD.h:50
VectorD<-1 > VectorKD
Definition: VectorD.h:411
Various important macros for implementing geometry.