IMP logo
IMP Reference Guide  develop.ba5ef9f91e,2024/02/26
The Integrative Modeling Platform
VectorBaseD.h
Go to the documentation of this file.
1 /**
2  * \file IMP/algebra/VectorBaseD.h \brief Simple D vector class.
3  *
4  * Copyright 2007-2022 IMP Inventors. All rights reserved.
5  *
6  */
7 
8 #ifndef IMPALGEBRA_VECTOR_BASE_D_H
9 #define IMPALGEBRA_VECTOR_BASE_D_H
10 
11 #include <IMP/algebra/algebra_config.h>
12 #include "GeometricPrimitiveD.h"
13 #include <IMP/check_macros.h>
14 #include <IMP/exception.h>
15 #include <IMP/random.h>
16 #include <IMP/utility.h>
17 #include <IMP/types.h>
18 #include <boost/random/variate_generator.hpp>
19 #include <boost/random/normal_distribution.hpp>
20 #include <boost/range.hpp>
21 #include <cereal/access.hpp>
22 #include "internal/vector.h"
23 
24 #include <limits>
25 #include <cmath>
26 #include <boost/random/normal_distribution.hpp>
27 
28 #if IMP_HAS_CHECKS >= IMP_INTERNAL
29 #define IMP_ALGEBRA_VECTOR_CHECK check_vector()
30 #define IMP_ALGEBRA_VECTOR_CHECK_OTHER(o) o.check_vector()
31 #else
32 #define IMP_ALGEBRA_VECTOR_CHECK
33 #define IMP_ALGEBRA_VECTOR_CHECK_OTHER(o)
34 #endif
35 
36 /* Should only need to check for "compatible" (same dimension) vectors
37  when using variable-dimension vectors; otherwise, it is checked at
38  compile time */
39 #if IMP_HAS_CHECKS >= IMP_USAGE
40 #define IMP_ALGEBRA_VECTOR_CHECK_INDEX(i) check_index(i)
41 #define IMP_ALGEBRA_VECTOR_CHECK_COMPATIBLE(o) \
42  if (D == -1) { check_compatible_vector(o); }
43 #else
44 #define IMP_ALGEBRA_VECTOR_CHECK_INDEX(i)
45 #define IMP_ALGEBRA_VECTOR_CHECK_COMPATIBLE(o)
46 #endif
47 
48 IMPALGEBRA_BEGIN_NAMESPACE
49 //! A Cartesian vector in D-dimensions.
50 /** Store a vector of Cartesian coordinates. It supports all expected
51  mathematical operators, including using * for the dot product.
52  \see Vector3D
53  \see Vector2D
54 
55  \geometry
56  */
57 template <int D>
58 class VectorBaseD : public GeometricPrimitiveD<D> {
59  friend class cereal::access;
60 
61  template<class Archive>
62  void serialize(Archive &ar) {
63  ar(data_);
64  }
65 
66  void check_vector() const {
67  IMP_INTERNAL_CHECK(!data_.get_is_null(),
68  "Attempt to use uninitialized vector.");
69  }
70  template <int OD>
71  void check_compatible_vector(const VectorBaseD<OD> &o) const {
73  IMP_USAGE_CHECK(o.get_dimension() == get_dimension(),
74  "Dimensions don't match: " << get_dimension() << " vs "
75  << o.get_dimension());
76  }
77  void check_index(unsigned int i) const {
78 #if IMP_HAS_CHECKS < IMP_INTERNAL
79  IMP_UNUSED(i);
80 #endif
81  IMP_INTERNAL_CHECK(i < data_.get_dimension(),
82  "Invalid component of vector requested: "
83  << i << " of " << get_dimension());
84  }
85 
86  public:
87  /** Will accept a list of floats from Python. */
88  template <class Range>
89  explicit VectorBaseD(const Range &r) {
90  if (D != -1 && static_cast<int>(boost::distance(r)) != D) {
91  IMP_THROW("Expected " << D << " but got " << boost::distance(r),
93  }
95  for(double f : r) {
96  IMP_UNUSED(f);
97  IMP_USAGE_CHECK(!IMP::is_nan(f), "NaN passed to constructor");
98  }
99  }
100  data_.set_coordinates(boost::begin(r), boost::end(r));
101  }
102 
103 #ifndef SWIG
104  template <class R>
105  VectorBaseD<D> &operator=(const R &r) {
106  if (D != -1 && static_cast<int>(boost::distance(r)) != D) {
107  IMP_THROW("Expected " << D << " but got " << boost::distance(r),
109  }
111  for(double f : r) {
112  IMP_USAGE_CHECK(!IMP::is_nan(f), "NaN passed in equals");
113  }
114  }
115  data_.set_coordinates(boost::begin(r), boost::end(r));
116  }
117  //! Return the ith Cartesian coordinate.
118  /** In 3D use [0] to get the x coordinate etc. */
119  inline double operator[](unsigned int i) const {
120  IMP_ALGEBRA_VECTOR_CHECK_INDEX(i);
121  IMP_ALGEBRA_VECTOR_CHECK;
122  return data_.get_data()[i];
123  }
124  //! Return the ith Cartesian coordinate.
125  /** In 3D use [0] to get the x coordinate etc. */
126  inline double &operator[](unsigned int i) {
127  IMP_ALGEBRA_VECTOR_CHECK_INDEX(i);
128  return data_.get_data()[i];
129  }
130 
131 #endif
132  //! Default constructor
134 
135  double get_scalar_product(const VectorBaseD &o) const {
136  IMP_ALGEBRA_VECTOR_CHECK;
137  IMP_ALGEBRA_VECTOR_CHECK_OTHER(o);
138  IMP_ALGEBRA_VECTOR_CHECK_COMPATIBLE(o);
139  double ret = 0;
140  for (unsigned int i = 0; i < get_dimension(); ++i) {
141  ret += operator[](i) * o.operator[](i);
142  }
143  return ret;
144  }
145 
146  double get_squared_magnitude() const {
147  // Could be equivalently written as get_scalar_product(*this)
148  // This should be faster though since checks on 'o' can be skipped
149  IMP_ALGEBRA_VECTOR_CHECK;
150  double ret = 0;
151  const double *data = get_data();
152  for (unsigned int i = 0; i < get_dimension(); ++i) {
153  ret += data[i] * data[i];
154  }
155  return ret;
156  }
157 
158  double get_magnitude() const { return std::sqrt(get_squared_magnitude()); }
159 
160  //! Return the distance between this and another vector
161  /** This is essentially identical to (v1 - v2).get_magnitude() but
162  may be slightly more efficient as it avoids creating a temporary
163  vector object. */
164  double get_distance(const VectorBaseD<D> &o) const {
165  IMP_ALGEBRA_VECTOR_CHECK_OTHER(o);
166  IMP_ALGEBRA_VECTOR_CHECK_COMPATIBLE(o);
167  const double *data = get_data(), *odata = o.get_data();
168  double ret = 0;
169  for (unsigned int i = 0; i < get_dimension(); ++i) {
170  ret += (odata[i] - data[i]) * (odata[i] - data[i]);
171  }
172  return std::sqrt(ret);
173  }
174 
175 #ifndef IMP_DOXYGEN
176  double operator*(const VectorBaseD<D> &o) const {
177  return get_scalar_product(o);
178  }
179 
180  VectorBaseD &operator+=(const VectorBaseD &o) {
181  IMP_ALGEBRA_VECTOR_CHECK;
182  IMP_ALGEBRA_VECTOR_CHECK_OTHER(o);
183  IMP_ALGEBRA_VECTOR_CHECK_COMPATIBLE(o);
184  for (unsigned int i = 0; i < get_dimension(); ++i) {
185  operator[](i) += o[i];
186  }
187  return *this;
188  }
189 
190  VectorBaseD &operator-=(const VectorBaseD &o) {
191  IMP_ALGEBRA_VECTOR_CHECK;
192  IMP_ALGEBRA_VECTOR_CHECK_OTHER(o);
193  IMP_ALGEBRA_VECTOR_CHECK_COMPATIBLE(o);
194  for (unsigned int i = 0; i < get_dimension(); ++i) {
195  operator[](i) -= o[i];
196  }
197  return *this;
198  }
199 
200  VectorBaseD &operator/=(double f) {
201  IMP_ALGEBRA_VECTOR_CHECK;
202  for (unsigned int i = 0; i < get_dimension(); ++i) {
203  operator[](i) /= f;
204  }
205  return *this;
206  }
207 
208  VectorBaseD &operator*=(double f) {
209  IMP_ALGEBRA_VECTOR_CHECK;
210  for (unsigned int i = 0; i < get_dimension(); ++i) {
211  operator[](i) *= f;
212  }
213  return *this;
214  }
215 
216  void show(std::ostream &out, std::string delim, bool parens = true) const {
217  IMP_ALGEBRA_VECTOR_CHECK;
218  if (parens) out << "(";
219  for (unsigned int i = 0; i < get_dimension(); ++i) {
220  out << operator[](i);
221  if (i != get_dimension() - 1) {
222  out << delim;
223  }
224  }
225  if (parens) out << ")";
226  }
227  IMP_SHOWABLE_INLINE(VectorBaseD, show(out, ", "););
228 #endif
229 
230 #ifndef SWIG
231  typedef double *iterator;
232  typedef const double *const_iterator;
233  iterator begin() { return data_.get_data(); }
234  iterator end() { return data_.get_data() + get_dimension(); }
235  const_iterator begin() const { return data_.get_data(); }
236  const_iterator end() const { return data_.get_data() + get_dimension(); }
237 
238  typedef double value_type;
239  typedef std::random_access_iterator_tag iterator_category;
240  typedef std::ptrdiff_t difference_type;
241  typedef double *pointer;
242  typedef double &reference;
243  typedef const double &const_reference;
244 
245  static const int DIMENSION = D;
246 #endif
247 
248 #ifndef SWIG
249  // For some reason, this method breaks IMP::atom::get_rmsd() in Python, so
250  // hide it from SWIG
251  Floats get_coordinates() const {
252  return Floats(begin(), end());
253  }
254 
255  //! Return a pointer to the data stored.
256  /** Useful for conversion to other types. */
257  const double *get_data() const { return data_.get_data(); }
258 #endif
259  unsigned int get_dimension() const { return data_.get_dimension(); }
260 
261  private:
262  internal::VectorData<double, D, false> data_;
263 };
264 
265 //! Returns a unit vector pointing at the same direction as this vector.
266 /**
267  @note If the magnitude of this vector is smaller than 1e-12
268  (an arbitrarily selected small number), returns a unit
269  vector pointing at a random direction.
270  */
271 template <class VT>
272 inline VT get_unit_vector(VT vt) {
273  static const double tiny_double =
274  256.0 * std::numeric_limits<double>::epsilon();
275  double mag = vt.get_magnitude();
276  if (mag > tiny_double) {
277  VT ret_value= vt/mag;
278  IMP_USAGE_CHECK(std::abs(ret_value.get_magnitude() - 1.0) < 256.0 * tiny_double,
279  "returned vector is not unit vector");
280  return ret_value;
281  } else {
282  // avoid division by zero - return random unit v
283  // NOTE: (1) avoids vector_generators / SphereD to prevent recursiveness
284  // (2) D might be -1, so use get_dimension()
285  static boost::variate_generator<RandomNumberGenerator,
286  boost::normal_distribution<> >
288  ::boost::normal_distribution<>(0, 1.0));
289  for (unsigned int i = 0; i < vt.get_dimension(); ++i) {
290  vt[i] = generator();
291  }
292  return get_unit_vector(vt);
293  }
294 }
295 
296 //! Returns the magnitude of vt and turns it to a unit vector in place.
297 /**
298  @note If the magnitude of this vector is smaller than 1e-12
299  (an arbitrarily selected small number), vt is turned into a unit
300  vector pointing at a random direction.
301  */
302 template <class VT>
304  const double tiny_double = 1e-12;
305  double mag = vt.get_magnitude();
306  if (mag > tiny_double) {
307  vt /= mag;
308  } else {
309  // avoid division by zero - return random unit v
310  // using get_unit_vector()
311  vt= get_unit_vector(vt);
312  }
313  return mag;
314 }
315 
316 IMPALGEBRA_END_NAMESPACE
317 
318 #endif /* IMPALGEBRA_VECTOR_BASE_D_H */
Base class for geometric types.
Basic types used by IMP.
#define IMP_USAGE_CHECK_VARIABLE(variable)
Definition: check_macros.h:190
#define IMP_SHOWABLE_INLINE(Name, how_to_show)
Declare the methods needed by an object that can be printed.
#define IMP_IF_CHECK(level)
Execute the code block if a certain level checks are on.
Definition: check_macros.h:104
VectorBaseD()
Default constructor.
Definition: VectorBaseD.h:133
IMP::Vector< Float > Floats
Standard way to pass a bunch of Float values.
Definition: types.h:46
VectorBaseD(const Range &r)
Definition: VectorBaseD.h:89
VT get_unit_vector(VT vt)
Returns a unit vector pointing at the same direction as this vector.
Definition: VectorBaseD.h:272
Exception definitions and assertions.
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
Definition: check_macros.h:139
double get_magnitude_and_normalize_in_place(VT &vt)
Returns the magnitude of vt and turns it to a unit vector in place.
Definition: VectorBaseD.h:303
Base class for geometric types.
#define IMP_UNUSED(variable)
double get_distance(const VectorBaseD< D > &o) const
Return the distance between this and another vector.
Definition: VectorBaseD.h:164
Various general useful functions for IMP.
std::ostream & show(Hierarchy h, std::ostream &out=std::cout)
Print the hierarchy using a given decorator to display each node.
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
Definition: check_macros.h:50
Helper macros for throwing and handling exceptions.
#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
Random number generators used by IMP.
const double * get_data() const
Return a pointer to the data stored.
Definition: VectorBaseD.h:257
A Cartesian vector in D-dimensions.
Definition: VectorBaseD.h:58
An exception for an invalid value being passed to IMP.
Definition: exception.h:136
double operator[](unsigned int i) const
Return the ith Cartesian coordinate.
Definition: VectorBaseD.h:119
RandomNumberGenerator random_number_generator
A shared non-GPU random number generator.
double & operator[](unsigned int i)
Return the ith Cartesian coordinate.
Definition: VectorBaseD.h:126