Vector.h

Go to the documentation of this file.

 
 
#ifndef SCMATH_GEOM_VECTORS_H
#define SCMATH_GEOM_VECTORS_H
 
 
#include <array>
#include <cassert>
#include <cmath>
#include <initializer_list>
#include <ostream>
#include <iostream>
 
#include <scmp/Math.h>
 
 
// This implementation is essentially a C++ translation of a Python library
// I wrote for Coursera's "Linear Algebra for Machine Learning" course. Many
// of the test vectors come from quiz questions in the class.
 
 
namespace scmp {
namespace geom {
 
 
template<typename T, size_t N>
class Vector {
public:
    Vector()
    {
        T           unitLength = (T) 1.0 / std::sqrt(N);
        for (size_t i          = 0; i < N; i++) {
            this->arr[i] = unitLength;
        }
 
        scmp::DefaultEpsilon(this->epsilon);
    }
 
 
    Vector(std::initializer_list<T> ilst)
    {
        assert(ilst.size() == N);
 
        scmp::DefaultEpsilon(this->epsilon);
        std::copy(ilst.begin(), ilst.end(), this->arr.begin());
    }
 
 
    T At(size_t index) const
    {
        if (index > this->arr.size()) {
            throw std::out_of_range("index " +
                        std::to_string(index) + " > " +
                        std::to_string(this->arr.size()));
        }
        return this->arr.at(index);
    }
 
 
    void Set(size_t index, T value)
    {
        if (index > this->arr.size()) {
            throw std::out_of_range("index " +
                        std::to_string(index) + " > " +
                        std::to_string(this->arr.size()));
        }
 
        this->arr[index] = value;
    }
 
 
 
 
 
    T Magnitude() const
    {
        T result = 0;
 
        for (size_t i = 0; i < N; i++) {
            result += (this->arr[i] * this->arr[i]);
        }
        return std::sqrt(result);
    }
 
 
    void
    SetEpsilon(T eps)
    {
        this->epsilon = eps;
    }
 
 
    bool
    IsZero() const
    {
        for (size_t i = 0; i < N; i++) {
            if (!scmp::WithinTolerance(this->arr[i], (T) 0.0, this->epsilon)) {
                return false;
            }
        }
        return true;
    }
 
 
    Vector
    UnitVector() const
    {
        return *this / this->Magnitude();
    }
 
 
    bool
    IsUnitVector() const
    {
        return scmp::WithinTolerance(this->Magnitude(), (T) 1.0, this->epsilon);
    }
 
 
    T
    Angle(const Vector<T, N> &other) const
    {
        auto unitA = this->UnitVector();
        auto unitB = other.UnitVector();
 
        // Can't compute angles with a zero vector.
        assert(!this->IsZero());
        assert(!other.IsZero());
        return static_cast<T>(std::acos(unitA * unitB));
    }
 
 
    bool
    IsParallel(const Vector<T, N> &other) const
    {
        if (this->IsZero() || other.IsZero()) {
            return true;
        }
 
        // If the two unit vectors are equal, the two vectors
        // lie on the same path.
        //
        // Context: this used to use Vector::Angle to check for
        // a zero angle between the two. However, the vagaries
        // of floating point math meant that while this worked
        // fine on Linux amd64 builds, it failed on Linux arm64
        // and MacOS builds. Parallel float vectors would have
        // an angle of ~0.0003 radians, while double vectors
        // would have an angle of +NaN. I suspect this is due to
        // tiny variations in floating point math, such that a dot
        // product of unit vectors would be just a hair over 1,
        // e.g. 1.000000001 - which would still fall outside the
        // domain of acos.
        auto unitA = this->UnitVector();
        auto unitB = other.UnitVector();
 
        return unitA == unitB;
    }
 
 
    bool
    IsOrthogonal(const Vector<T, N> &other) const
    {
        if (this->IsZero() || other.IsZero()) {
            return true;
        }
 
        return scmp::WithinTolerance(*this * other, (T) 0.0, this->epsilon);
    }
 
 
    Vector
    ProjectParallel(const Vector<T, N> &basis) const
    {
        Vector<T, N> unit_basis = basis.UnitVector();
 
        return unit_basis * (*this * unit_basis);
    }
 
 
    Vector
    ProjectOrthogonal(const Vector<T, N> &basis)
    {
        Vector<T, N> spar = this->ProjectParallel(basis);
        return *this - spar;
    }
 
 
    Vector
    Cross(const Vector<T, N> &other) const
    {
        assert(N == 3);
        if (N != 3) {
            throw std::out_of_range("Cross-product can only called on Vector<T, 3>.");
        }
 
        return Vector<T, N>{
            (this->arr[1] * other.arr[2]) - (other.arr[1] * this->arr[2]),
            -((this->arr[0] * other.arr[2]) - (other.arr[0] * this->arr[2])),
            (this->arr[0] * other.arr[1]) - (other.arr[0] * this->arr[1])
        };
    }
 
 
    Vector
    operator+(const Vector<T, N> &other) const
    {
        Vector<T, N> vec;
 
        for (size_t i = 0; i < N; i++) {
            vec.arr[i] = this->arr[i] + other.arr[i];
        }
 
        return vec;
    }
 
 
    Vector
    operator-(const Vector<T, N> &other) const
    {
        Vector<T, N> vec;
 
        for (size_t i = 0; i < N; i++) {
            vec.arr[i] = this->arr[i] - other.arr[i];
        }
 
        return vec;
    }
 
 
    Vector
    operator*(const T k) const
    {
        Vector<T, N> vec;
 
        for (size_t i = 0; i < N; i++) {
            vec.arr[i] = this->arr[i] * k;
        }
 
        return vec;
    }
 
 
    Vector
    operator/(const T k) const
    {
        Vector<T, N> vec;
 
        for (size_t i = 0; i < N; i++) {
            vec.arr[i] = this->arr[i] / k;
        }
 
        return vec;
    }
 
 
    T
    operator*(const Vector<T, N> &other) const
    {
        T result = 0;
 
        for (size_t i = 0; i < N; i++) {
            result += (this->arr[i] * other.arr[i]);
        }
 
        return result;
    }
 
 
    bool
    operator==(const Vector<T, N> &other) const
    {
        for (size_t i = 0; i < N; i++) {
            if (!scmp::WithinTolerance(this->arr[i], other.arr[i], this->epsilon)) {
                return false;
            }
        }
        return true;
    }
 
 
    bool
    operator!=(const Vector<T, N> &other) const
    {
        return !(*this == other);
    }
 
 
    const T &
    operator[](size_t i) const
    {
        return this->arr[i];
    }
 
 
    friend std::ostream &
    operator<<(std::ostream &outs, const Vector<T, N> &vec)
    {
        outs << "<";
        for (size_t i = 0; i < N; i++) {
            outs << vec.arr[i];
            if (i < (N - 1)) {
                outs << ", ";
            }
        }
        outs << ">";
        return outs;
    }
 
private:
    static const size_t dim = N;
    T                   epsilon;
    std::array<T, N>    arr;
};
 
 
 
typedef Vector<float, 2> Vector2F;
 
typedef Vector<float, 3> Vector3F;
 
typedef Vector<float, 4> Vector4F;
 
typedef Vector<double, 2> Vector2D;
 
typedef Vector<double, 3> Vector3D;
 
typedef Vector<double, 4> Vector4D;
 
 
} // namespace geom
} // namespace scmp
 
 
#endif // SCMATH_GEOM_VECTORS_H