IMP logo
IMP Reference Guide  develop.d247202c1c,2019/10/23
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-2019 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/types.h>
14 #include <IMP/check_macros.h>
15 #include <IMP/exception.h>
16 #include <IMP/utility.h>
17 #include <IMP/InputAdaptor.h>
18 #include <IMP/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 the basis vector for the given coordinate
235 template <int D>
236 inline VectorD<D> get_basis_vector_kd(int Di, unsigned int coordinate) {
237  IMP_USAGE_CHECK(D == Di, "D must be equal");
238  IMP_UNUSED(Di);
239  return get_basis_vector_d<D>(coordinate);
240 }
241 
242 //! Return a dynamically sized basis vector
243 template <>
244 inline VectorD<-1> get_basis_vector_kd(int D, unsigned int coordinate) {
245  IMP_USAGE_CHECK(D > 0, "D must be positive");
246  IMP_USAGE_CHECK(coordinate < static_cast<unsigned int>(D),
247  "There are only " << D << " basis vectors");
248  boost::scoped_array<double> vs(new double[D]);
249  for (int i = 0; i < D; ++i) {
250  if (i == static_cast<int>(coordinate))
251  vs[i] = 1;
252  else
253  vs[i] = 0;
254  }
255  return VectorD<-1>(vs.get(), vs.get() + D);
256 }
257 
258 //! Return a dynamically sized basis vector
259 inline VectorD<-1> get_basis_vector_kd(int D, unsigned int coordinate) {
260  return get_basis_vector_kd<-1>(D, coordinate);
261 }
262 
263 //! Return a vector of zeros
264 template <int D>
266  IMP_USAGE_CHECK(D > 0, "D must be positive");
267  VectorD<D> ret;
268  for (int i = 0; i < D; ++i) {
269  ret[i] = 0;
270  }
271  return ret;
272 }
273 
274 //! Return a dynamically sized vector of zeros
275 template <int D>
277  IMP_USAGE_CHECK(D == Di, "D must be equal");
278  IMP_UNUSED(Di);
279  return get_zero_vector_d<D>();
280 }
281 
282 //! Return a dynamically sized vector of zeros
283 template<>
284 inline VectorD<-1> get_zero_vector_kd(int D) {
285  IMP_USAGE_CHECK(D > 0, "D must be positive");
286  Floats vs(D, 0);
287  return VectorD<-1>(vs.begin(), vs.end());
288 }
289 
290 //! Return a dynamically sized vector of zeros
291 inline VectorD<-1> get_zero_vector_kd(int D) {
292  return get_zero_vector_kd<-1>(D);
293 }
294 
295 //! Return a vector of ones (or another constant)
296 template <int D>
297 inline VectorD<D> get_ones_vector_d(double v = 1) {
298  IMP_USAGE_CHECK(D > 0, "D must be positive");
299  VectorD<D> ret;
300  for (unsigned int i = 0; i < D; ++i) {
301  ret[i] = v;
302  }
303  return ret;
304 }
305 
306 //! Return a vector of ones (or another constant)
307 /** Di must equal D. */
308 template <int D>
309 inline VectorD<D> get_ones_vector_kd(unsigned int Di, double v = 1) {
310  IMP_USAGE_CHECK(D == Di, "D must be equal");
311  IMP_UNUSED(Di);
312  return get_ones_vector_d<D>(v);
313 }
314 
315 //! Return a vector of ones (or another constant)
316 template <>
317 inline VectorD<-1> get_ones_vector_kd(unsigned int D, double v) {
318  IMP_USAGE_CHECK(D > 0, "D must be positive");
319  boost::scoped_array<double> vv(new double[D]);
320  for (unsigned int i = 0; i < D; ++i) {
321  vv[i] = v;
322  }
323  return VectorD<-1>(vv.get(), vv.get() + D);
324 }
325 
326 //! Return a dynamically sized vector of zeros
327 inline VectorD<-1> get_ones_vector_kd(unsigned int D, double v = 1) {
328  return get_ones_vector_kd<-1>(D, v);
329 }
330 
331 #ifndef SWIG
332 
333 /** \name Norms
334  We define a number of standard, \f$L^p\f$, norms on VectorD.
335  - \f$L^1\f$ is the Manhattan distance, the sum of the components
336  - \f$L^2\f$ is the standard Euclidean length
337  - \f$L^{\inf}\f$ is the maximum of the components
338  @{
339 */
340 
341 template <int D>
342 inline double get_l2_norm(const VectorD<D> &v) {
343  return v.get_magnitude();
344 }
345 
346 template <int D>
347 inline double get_l1_norm(const VectorD<D> &v) {
348  double n = std::abs(v[0]);
349  for (unsigned int i = 1; i < v.get_dimension(); ++i) {
350  n += std::abs(v[i]);
351  }
352  return n;
353 }
354 
355 template <int D>
356 inline double get_linf_norm(const VectorD<D> &v) {
357  double n = std::abs(v[0]);
358  for (unsigned int i = 1; i < v.get_dimension(); ++i) {
359  n = std::max(n, std::abs(v[i]));
360  }
361  return n;
362 }
363 
364 /** @} */
365 
366 #ifndef IMP_DOXYGEN
367 
368 template <int D>
369 struct SpacesIO {
370  const VectorD<D> &v_;
371  SpacesIO(const VectorD<D> &v) : v_(v) {}
372 };
373 
374 template <int D>
375 struct CommasIO {
376  const VectorD<D> &v_;
377  CommasIO(const VectorD<D> &v) : v_(v) {}
378 };
379 template <int D>
380 inline std::ostream &operator<<(std::ostream &out, const SpacesIO<D> &s) {
381  s.v_.show(out, " ", false);
382  return out;
383 }
384 template <int D>
385 inline std::ostream &operator<<(std::ostream &out, const CommasIO<D> &s) {
386  s.v_.show(out, ", ", false);
387  return out;
388 }
389 
390 //! Use this before outputting to delimited vector entries with a space
391 /** std::cout << spaces_io(v);
392  produces "1.0 2.0 3.0"
393  \see VectorD
394  */
395 template <int D>
396 inline SpacesIO<D> spaces_io(const VectorD<D> &v) {
397  return SpacesIO<D>(v);
398 }
399 
400 //! Use this before outputting to delimited vector entries with a comma
401 /** std::cout << commas_io(v);
402  produces "1.0, 2.0, 3.0"
403  \see VectorD
404  */
405 template <int D>
406 inline CommasIO<D> commas_io(const VectorD<D> &v) {
407  return CommasIO<D>(v);
408 }
409 #endif // doxygen
410 
411 #endif // swig
412 /** 1D vector typedef for swig */
414 /** 1D vectors typedef for swig */
416 /** 2D vector typedef for swig */
418 /** 2D vectors typedef for swig */
420 /** 3D vector typedef for swig */
422 /** 3D vectors typedef for swig */
424 /** 4D vector typedef for swig */
426 /** 4D vectors typedef for swig */
428 /** 5D vector typedef for swig */
430 /** 5D vectors typedef for swig */
432 /** 6D vector typedef for swig */
434 /** 6D vector typedef for swig */
436 /** KD vector typedef for swig */
437 typedef VectorD<-1> VectorKD;
438 /** KD vectors typedef for swig */
439 typedef Vector<VectorD<-1> > VectorKDs;
440 
441 #ifndef SWIG
442 /** \see VectorD \genericgeometry */
443 template <class C>
444 inline const VectorD<C::DIMENSION> &get_vector_geometry(const C &g) {
445  return g;
446 }
447 /** \see VectorD \genericgeometry */
448 template <class C, class E>
449 inline void set_vector_geometry(C &g, const E &v) {
450  g = v;
451 }
452 #endif
453 
454 //! Return the vector that is the elementwise product of the two.
455 /** \see VectorD */
456 template <int D>
458  const algebra::VectorD<D> &b) {
459  VectorD<D> ret(a);
460  for (unsigned int i = 0; i < ret.get_dimension(); ++i) {
461  ret[i] *= b[i];
462  }
463  return ret;
464 }
465 
466 //! Return the vector that is the elementwise product of the two.
467 /** \see VectorD */
468 template <int D>
470  const algebra::VectorD<D> &b) {
471  IMP_USAGE_CHECK(a.size() == b.get_dimension(), "Dimensions don't match,");
472  VectorD<D> ret(b);
473  for (unsigned int i = 0; i < ret.get_dimension(); ++i) {
474  ret[i] *= a[i];
475  }
476  return ret;
477 }
478 
479 IMPALGEBRA_END_NAMESPACE
480 
481 #endif /* IMPALGEBRA_VECTOR_D_H */
Basic types used by IMP.
VectorD< 6 > Vector6D
Definition: VectorD.h:433
VectorD< 2 > Vector2D
Definition: VectorD.h:417
VectorD< 5 > Vector5D
Definition: VectorD.h:429
Vector< VectorD< 1 > > Vector1Ds
Definition: VectorD.h:415
Vector< VectorD< 3 > > Vector3Ds
Definition: VectorD.h:423
double get_squared_distance(const VectorD< D > &v1, const VectorD< D > &v2)
Compute the squared distance between two vectors.
Definition: VectorD.h:201
Vector< VectorD< 5 > > Vector5Ds
Definition: VectorD.h:431
Exception definitions and assertions.
VectorD< 4 > Vector4D
Definition: VectorD.h:425
Vector< VectorD< 2 > > Vector2Ds
Definition: VectorD.h:419
const VectorD< C::DIMENSION > & get_vector_geometry(const C &g)
Definition: VectorD.h:444
Basic types used by IMP.
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:469
int compare(const VectorD< D > &a, const VectorD< D > &b)
lexicographic comparison of two vectors
Definition: VectorD.h:179
VectorD< D > get_zero_vector_d()
Return a vector of zeros.
Definition: VectorD.h:265
Vector< VectorD<-1 > > VectorKDs
Definition: VectorD.h:439
VectorD< D > get_basis_vector_d(unsigned int coordinate)
Return the basis vector for the given coordinate.
Definition: VectorD.h:222
VectorD< 1 > Vector1D
Definition: VectorD.h:413
void set_vector_geometry(C &g, const E &v)
Definition: VectorD.h:449
Vector< VectorD< 6 > > Vector6Ds
Definition: VectorD.h:435
For backwards compatibility.
Vector< VectorD< 4 > > Vector4Ds
Definition: VectorD.h:427
VectorD<-1 > get_ones_vector_kd(unsigned int D, double v)
Return a vector of ones (or another constant)
Definition: VectorD.h:317
Simple D vector class.
VectorD< D > get_ones_vector_d(double v=1)
Return a vector of ones (or another constant)
Definition: VectorD.h:297
Exception definitions and assertions.
VectorD< 3 > Vector3D
Definition: VectorD.h:421
VectorD<-1 > get_basis_vector_kd(int D, unsigned int coordinate)
Return a dynamically sized basis vector.
Definition: VectorD.h:244
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
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:437
Various helper macros.
VectorD<-1 > get_zero_vector_kd(int D)
Return a dynamically sized vector of zeros.
Definition: VectorD.h:284