scsl  1.1.1
Shimmering Clarity Standard Library
Vector.h
Go to the documentation of this file.
1 
23 #ifndef SCMATH_GEOM_VECTORS_H
24 #define SCMATH_GEOM_VECTORS_H
25 
26 
27 #include <array>
28 #include <cassert>
29 #include <cmath>
30 #include <initializer_list>
31 #include <ostream>
32 #include <iostream>
33 
34 #include <scmp/Math.h>
35 
36 
37 // This implementation is essentially a C++ translation of a Python library
38 // I wrote for Coursera's "Linear Algebra for Machine Learning" course. Many
39 // of the test vectors come from quiz questions in the class.
40 
41 
42 namespace scmp {
43 namespace geom {
44 
45 
56 template<typename T, size_t N>
57 class Vector {
58 public:
61  {
62  T unitLength = (T) 1.0 / std::sqrt(N);
63  for (size_t i = 0; i < N; i++) {
64  this->arr[i] = unitLength;
65  }
66 
67  scmp::DefaultEpsilon(this->epsilon);
68  }
69 
70 
77  Vector(std::initializer_list<T> ilst)
78  {
79  assert(ilst.size() == N);
80 
81  scmp::DefaultEpsilon(this->epsilon);
82  std::copy(ilst.begin(), ilst.end(), this->arr.begin());
83  }
84 
85 
92  T At(size_t index) const
93  {
94  if (index > this->arr.size()) {
95  throw std::out_of_range("index " +
96  std::to_string(index) + " > " +
97  std::to_string(this->arr.size()));
98  }
99  return this->arr.at(index);
100  }
101 
102 
111  void Set(size_t index, T value)
112  {
113  if (index > this->arr.size()) {
114  throw std::out_of_range("index " +
115  std::to_string(index) + " > " +
116  std::to_string(this->arr.size()));
117  }
118 
119  this->arr[index] = value;
120  }
121 
122 
123 
124 
125 
129  T Magnitude() const
130  {
131  T result = 0;
132 
133  for (size_t i = 0; i < N; i++) {
134  result += (this->arr[i] * this->arr[i]);
135  }
136  return std::sqrt(result);
137  }
138 
139 
147  void
148  SetEpsilon(T eps)
149  {
150  this->epsilon = eps;
151  }
152 
153 
157  bool
158  IsZero() const
159  {
160  for (size_t i = 0; i < N; i++) {
161  if (!scmp::WithinTolerance(this->arr[i], (T) 0.0, this->epsilon)) {
162  return false;
163  }
164  }
165  return true;
166  }
167 
168 
172  Vector
173  UnitVector() const
174  {
175  return *this / this->Magnitude();
176  }
177 
178 
182  bool
183  IsUnitVector() const
184  {
185  return scmp::WithinTolerance(this->Magnitude(), (T) 1.0, this->epsilon);
186  }
187 
188 
193  T
194  Angle(const Vector<T, N> &other) const
195  {
196  auto unitA = this->UnitVector();
197  auto unitB = other.UnitVector();
198 
199  // Can't compute angles with a zero vector.
200  assert(!this->IsZero());
201  assert(!other.IsZero());
202  return static_cast<T>(std::acos(unitA * unitB));
203  }
204 
205 
210  bool
211  IsParallel(const Vector<T, N> &other) const
212  {
213  if (this->IsZero() || other.IsZero()) {
214  return true;
215  }
216 
217  // If the two unit vectors are equal, the two vectors
218  // lie on the same path.
219  //
220  // Context: this used to use Vector::Angle to check for
221  // a zero angle between the two. However, the vagaries
222  // of floating point math meant that while this worked
223  // fine on Linux amd64 builds, it failed on Linux arm64
224  // and MacOS builds. Parallel float vectors would have
225  // an angle of ~0.0003 radians, while double vectors
226  // would have an angle of +NaN. I suspect this is due to
227  // tiny variations in floating point math, such that a dot
228  // product of unit vectors would be just a hair over 1,
229  // e.g. 1.000000001 - which would still fall outside the
230  // domain of acos.
231  auto unitA = this->UnitVector();
232  auto unitB = other.UnitVector();
233 
234  return unitA == unitB;
235  }
236 
237 
243  bool
244  IsOrthogonal(const Vector<T, N> &other) const
245  {
246  if (this->IsZero() || other.IsZero()) {
247  return true;
248  }
249 
250  return scmp::WithinTolerance(*this * other, (T) 0.0, this->epsilon);
251  }
252 
253 
259  Vector
260  ProjectParallel(const Vector<T, N> &basis) const
261  {
262  Vector<T, N> unit_basis = basis.UnitVector();
263 
264  return unit_basis * (*this * unit_basis);
265  }
266 
267 
275  Vector
277  {
278  Vector<T, N> spar = this->ProjectParallel(basis);
279  return *this - spar;
280  }
281 
282 
291  Vector
292  Cross(const Vector<T, N> &other) const
293  {
294  assert(N == 3);
295  if (N != 3) {
296  throw std::out_of_range("Cross-product can only called on Vector<T, 3>.");
297  }
298 
299  return Vector<T, N>{
300  (this->arr[1] * other.arr[2]) - (other.arr[1] * this->arr[2]),
301  -((this->arr[0] * other.arr[2]) - (other.arr[0] * this->arr[2])),
302  (this->arr[0] * other.arr[1]) - (other.arr[0] * this->arr[1])
303  };
304  }
305 
306 
312  Vector
313  operator+(const Vector<T, N> &other) const
314  {
315  Vector<T, N> vec;
316 
317  for (size_t i = 0; i < N; i++) {
318  vec.arr[i] = this->arr[i] + other.arr[i];
319  }
320 
321  return vec;
322  }
323 
324 
330  Vector
331  operator-(const Vector<T, N> &other) const
332  {
333  Vector<T, N> vec;
334 
335  for (size_t i = 0; i < N; i++) {
336  vec.arr[i] = this->arr[i] - other.arr[i];
337  }
338 
339  return vec;
340  }
341 
342 
347  Vector
348  operator*(const T k) const
349  {
350  Vector<T, N> vec;
351 
352  for (size_t i = 0; i < N; i++) {
353  vec.arr[i] = this->arr[i] * k;
354  }
355 
356  return vec;
357  }
358 
359 
364  Vector
365  operator/(const T k) const
366  {
367  Vector<T, N> vec;
368 
369  for (size_t i = 0; i < N; i++) {
370  vec.arr[i] = this->arr[i] / k;
371  }
372 
373  return vec;
374  }
375 
376 
381  T
382  operator*(const Vector<T, N> &other) const
383  {
384  T result = 0;
385 
386  for (size_t i = 0; i < N; i++) {
387  result += (this->arr[i] * other.arr[i]);
388  }
389 
390  return result;
391  }
392 
393 
399  bool
400  operator==(const Vector<T, N> &other) const
401  {
402  for (size_t i = 0; i < N; i++) {
403  if (!scmp::WithinTolerance(this->arr[i], other.arr[i], this->epsilon)) {
404  return false;
405  }
406  }
407  return true;
408  }
409 
410 
416  bool
417  operator!=(const Vector<T, N> &other) const
418  {
419  return !(*this == other);
420  }
421 
422 
435  const T &
436  operator[](size_t i) const
437  {
438  return this->arr[i];
439  }
440 
441 
447  friend std::ostream &
448  operator<<(std::ostream &outs, const Vector<T, N> &vec)
449  {
450  outs << "<";
451  for (size_t i = 0; i < N; i++) {
452  outs << vec.arr[i];
453  if (i < (N - 1)) {
454  outs << ", ";
455  }
456  }
457  outs << ">";
458  return outs;
459  }
460 
461 private:
462  static const size_t dim = N;
463  T epsilon;
464  std::array<T, N> arr;
465 };
466 
470 
475 
479 
483 
487 
491 
495 
499 
500 
501 } // namespace geom
502 } // namespace scmp
503 
504 
505 #endif // SCMATH_GEOM_VECTORS_H
Common math functions.
Vectors represent a direction and Magnitude.
Definition: Vector.h:57
Vector ProjectOrthogonal(const Vector< T, N > &basis)
Project this vector perpendicularly onto some basis vector.
Definition: Vector.h:276
Vector operator/(const T k) const
Scalar division.
Definition: Vector.h:365
T Magnitude() const
Compute the length of the vector.
Definition: Vector.h:129
Vector Cross(const Vector< T, N > &other) const
Compute the cross product of two vectors.
Definition: Vector.h:292
bool IsUnitVector() const
Determine if this is a unit vector.
Definition: Vector.h:183
friend std::ostream & operator<<(std::ostream &outs, const Vector< T, N > &vec)
Write a vector a stream in the form "<i, j, ...>".
Definition: Vector.h:448
void Set(size_t index, T value)
Set a new value for the vector.
Definition: Vector.h:111
bool IsZero() const
Determine whether this is a zero vector.
Definition: Vector.h:158
bool operator==(const Vector< T, N > &other) const
Vector equivalence.
Definition: Vector.h:400
void SetEpsilon(T eps)
Set equivalence tolerance.
Definition: Vector.h:148
Vector operator*(const T k) const
Scalar multiplication.
Definition: Vector.h:348
T At(size_t index) const
Return the element At index i.
Definition: Vector.h:92
Vector UnitVector() const
Obtain the unit vector for this vector.
Definition: Vector.h:173
bool IsParallel(const Vector< T, N > &other) const
Determine whether two vectors are parallel.
Definition: Vector.h:211
Vector operator-(const Vector< T, N > &other) const
Vector subtraction.
Definition: Vector.h:331
T operator*(const Vector< T, N > &other) const
Compute the Dot product between two vectors.
Definition: Vector.h:382
T Angle(const Vector< T, N > &other) const
Compute the Angle between two vectors.
Definition: Vector.h:194
const T & operator[](size_t i) const
Array indexing into vector.
Definition: Vector.h:436
Vector()
Construct a unit vector of a given type and size.
Definition: Vector.h:60
bool operator!=(const Vector< T, N > &other) const
Vector non-equivalence.
Definition: Vector.h:417
Vector operator+(const Vector< T, N > &other) const
Vector addition.
Definition: Vector.h:313
Vector ProjectParallel(const Vector< T, N > &basis) const
Project this vector onto some basis vector.
Definition: Vector.h:260
Vector(std::initializer_list< T > ilst)
Construct a Vector with initial values.
Definition: Vector.h:77
bool IsOrthogonal(const Vector< T, N > &other) const
Determine if two vectors are orthogonal or perpendicular to each other.
Definition: Vector.h:244
Vector< float, 3 > Vector3F
Type alias for a three-dimensional float vector.
Definition: Vector.h:482
Vector< double, 2 > Vector2D
Type alias for a two-dimensional double vector.
Definition: Vector.h:490
Vector< double, 3 > Vector3D
Type alias for a three-dimensional double vector.
Definition: Vector.h:494
Vector< double, 4 > Vector4D
Type alias for a four-dimensional double vector.
Definition: Vector.h:498
Vector< float, 4 > Vector4F
Type alias for a four-dimensional float vector.
Definition: Vector.h:486
Vector< float, 2 > Vector2F
Type alias for a two-dimensional float vector.
Definition: Vector.h:478
Shimmering Clarity Math & Physics toolkit.
Definition: estimation.h:31
void DefaultEpsilon(double &epsilon)
Get the default epsilon value.