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