Vector.h

Go to the documentation of this file.

/*  YAP - Yet another PWA toolkit
    Copyright 2015, Technische Universitaet Muenchen,
    Authors: Daniel Greenwald, Johannes Rauch

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/


#ifndef yap_Vector_h
#define yap_Vector_h

#include "fwd/Vector.h"

#include "Matrix.h"

#include <algorithm>
#include <array>
#include <cmath>
#include <iterator>
#include <numeric>
#include <ostream>
#include <string>
#include <type_traits>

namespace yap {

template <typename T, size_t N>
class VectorIterator : public std::iterator<std::random_access_iterator_tag, T, typename std::array<T, N>::difference_type>
{
public:

    VectorIterator(typename std::array<T, N>::iterator it)
        : Iterator_(it)
    {}

    bool operator!=(VectorIterator b) const
    { return this->Iterator_ != b.Iterator_; }

    T& operator*()
    { return *(this->Iterator_); }

    T* operator->()
    { return std::array<T, N>::operator->(this->Iterator_); }

    VectorIterator& operator++()
    { ++(this->Iterator_); return *this; }

    VectorIterator operator++(int)
    { VectorIterator it(*this); (this->Iterator_)++; return it; }

    VectorIterator& operator--()
    { --(this->Iterator_); return *this; }

    VectorIterator operator--(int)
    { VectorIterator it(*this); --(this->Iterator_); return it; }

    VectorIterator& operator+=(typename VectorIterator::difference_type n)
    { this->Iterator_ += n; return *this; }

    VectorIterator& operator-=(typename VectorIterator::difference_type n)
    { this->Iterator_ -= n; return *this; }

    friend const typename VectorIterator::difference_type operator-(VectorIterator lhs, VectorIterator rhs)
    { return lhs.Iterator_ - rhs.Iterator_; }

    VectorIterator operator[](typename VectorIterator::difference_type n) const
    { VectorIterator it(Iterator_[n]); return it; }

    friend const bool operator<(VectorIterator lhs, VectorIterator rhs)
    { return lhs.Iterator_ < rhs.Iterator_; }

    friend bool operator>(VectorIterator lhs, VectorIterator rhs)
    { return lhs.Iterator_ > rhs.Iterator_; }

private:

    typename std::array<T, N>::iterator Iterator_;
};

template <typename T, size_t N>
const VectorIterator<T, N> operator+(VectorIterator<T, N> a, typename VectorIterator<T, N>::difference_type n)
{ return (VectorIterator<T, N>(a) += n); }

template <typename T, size_t N>
const VectorIterator<T, N> operator+(typename VectorIterator<T, N>::difference_type n, VectorIterator<T, N> a)
{ return a + n; }

template <typename T, size_t N>
const VectorIterator<T, N> operator-(VectorIterator<T, N> a, typename VectorIterator<T, N>::difference_type n)
{ return (VectorIterator<T, N>(a) -= n); }

template <typename T, size_t N>
bool operator>=(VectorIterator<T, N> a, VectorIterator<T, N> b)
{ return !(a < b); }

template <typename T, size_t N>
bool operator<=(VectorIterator<T, N> a, VectorIterator<T, N> b)
{ return !(b > a); }

template <typename T, size_t N, typename E>
class VectorExpression
{
public:

    constexpr T operator[](size_t i) const
    { return static_cast<const E&>(*this)[i]; }

    constexpr size_t size() const
    { return N; }

    operator E& ()
    { return static_cast<E&>(*this); }

    constexpr operator const E& () const
    { return static_cast<const E&>(*this); }

};

template <typename T, size_t N>
class Vector : public VectorExpression<T, N, Vector<T, N> >
{
public:
    constexpr Vector(const std::array<T, N>& v) noexcept : Elements_(v) {}

    template <typename E>
    Vector(const VectorExpression<T, N, E>& V)
    { for (size_t i = 0; i < V.size(); ++i) Elements_[i] = V[i]; }

    Vector(T element = 0)
    { Elements_.fill(element); }

    T& operator[](size_t i)
    { return Elements_[i]; }

    const T operator[](size_t i) const
    { return Elements_[i]; }

    T& front()
    { return Elements_.front(); }

    const T& front() const
    { return Elements_.front(); }

    VectorIterator<T, N> begin()
    { return VectorIterator<T, N>(Elements_.begin()); }

    VectorIterator<const T, N> begin() const
    { return VectorIterator<const T, N>(Elements_.begin()); }

    VectorIterator<T, N> end()
    { return VectorIterator<T, N>(Elements_.end()); }

    VectorIterator<const T, N> end() const
    { return VectorIterator<const T, N>(Elements_.end()); }

    virtual auto operator*(const Vector<T, N>& B) const
    -> decltype(T(0) * T(0))
    { return std::inner_product(Elements_.begin(), Elements_.end(), B.Elements_.begin(), T(0) * T(0)); }

    friend bool operator==(const Vector<T, N>& lhs, const Vector<T, N>& rhs)
    { return lhs.Elements_ == rhs.Elements_; }

    friend const T angle(const Vector<T, N>& A, const Vector<T, N>& B)
    {
        T arg = A * B / abs(A) / abs(B);

        // correct for arg just outside boundary (by numerical precision)
        if (std::isfinite(arg) and fabs(arg) > 1)
            arg = (arg > 0) ? (T)1 : (T) - 1;

        return acos(arg);
    }

    friend constexpr T norm(const Vector<T, N>& A)
    { return A * A; }


private:

    std::array<T, N> Elements_;
};

template <typename T, size_t N, typename E1, typename E2>
class VectorAddition : public VectorExpression<T, N, VectorAddition<T, N, E1, E2> >
{
public:
    constexpr VectorAddition(const VectorExpression<T, N, E1>& a, const VectorExpression<T, N, E2>& b)
        : VectorExpression<T, N, VectorAddition<T, N, E1, E2> >(), A_(a), B_(b)
    {}

    constexpr T operator[](size_t i) const
    { return A_[i] + B_[i]; }

private:

    const E1& A_;

    const E2& B_;

};

template <typename T, size_t N, typename E1, typename E2>
constexpr VectorAddition<T, N, E1, E2> operator+(const VectorExpression<T, N, E1>& a, const VectorExpression<T, N, E2>& b)
{ return VectorAddition<T, N, E1, E2>(a, b); }

template <typename T, size_t N, typename E1, typename E2>
class VectorSubtraction : public VectorExpression<T, N, VectorSubtraction<T, N, E1, E2> >
{
public:
    constexpr VectorSubtraction(const VectorExpression<T, N, E1>& a, const VectorExpression<T, N, E2>& b)
        : VectorExpression<T, N, VectorSubtraction<T, N, E1, E2> >(), A_(a), B_(b)
    {}

    constexpr T operator[](size_t i) const
    { return A_[i] - B_[i]; }

private:

    const E1& A_;

    const E2& B_;

};

template <typename T, size_t N, typename E1, typename E2>
constexpr VectorSubtraction<T, N, E1, E2> operator-(const VectorExpression<T, N, E1>& a, const VectorExpression<T, N, E2>& b)
{ return VectorSubtraction<T, N, E1, E2>(a, b); }

template <typename T, size_t N>
std::string to_string(const Vector<T, N>& V)
{
    if (N == 0)
        return "(empty)";
    std::string s = "(";
    std::for_each(V.begin(), V.end(), [&](const T & t) {s += std::to_string(t) + ", ";});
    s.erase(s.size() - 2, 2);
    s += ")";
    return s;
}

template <typename T, size_t N>
std::ostream& operator<<(std::ostream& os, const Vector<T, N>& V)
{ os << to_string(V); return os; }

template <typename T, size_t N>
Vector<T, N>& operator+=(Vector<T, N>& A, const Vector<T, N>& B)
{ std::transform(A.begin(), A.end(), B.begin(), A.begin(), [](const T & a, const T & b) {return a + b;}); return A; }

template <typename T, size_t N>
constexpr Vector<T, N> operator-(const Vector<T, N>& V)
{ return T(-1) * V; }

template <typename T, size_t N>
Vector<T, N>& operator-=(Vector<T, N>& A, const Vector<T, N>& B)
{ std::transform(A.begin(), A.end(), B.begin(), A.begin(), [](const T & a, const T & b) {return a - b;}); return A; }

template <typename T, size_t N>
Vector<T, N>& operator*=(Vector<T, N>& A, const T& B)
{ std::transform(A.begin(), A.end(), A.begin(), [&](const T & v) {return B * v;}); return A; }

template <typename T, size_t N>
const Vector<T, N> operator*(const Vector<T, N>& A, const T& c)
{ auto v = A; return v *= c; return v; }

template <typename T, size_t N>
constexpr Vector<T, N> operator*(const T& c, const Vector<T, N>& A)
{ return A * c; }

template <typename T, size_t N>
Vector<T, N>& operator/=(Vector<T, N>& A, const T& B)
{ std::transform(A.begin(), A.end(), A.begin(), [&](const T & v) {return v / B;}); return A; }

template <typename T, size_t N>
const Vector<T, N> operator/(const Vector<T, N>& A, const T& c)
{ auto v = A; v /= c; return v; }

template <typename T, size_t N>
constexpr T abs(const Vector<T, N>& A)
{ return sqrt(norm(A)); }

template <typename T, size_t N>
const Vector<T, N> unit(const Vector<T, N>& V)
{ T a = abs(V); return (a == 0) ? V : (T(1) / a) * V; }

template <typename T, size_t R, size_t C>
const Vector<T, R> operator*(const Matrix<T, R, C>& M, const Vector<T, C>& V)
{
    Vector<T, R> v;
    for (size_t r = 0; r < R; ++r)
        for (size_t c = 0; c < C; ++c)
            v[r] += M[r][c] * V[c];
    return v;
}

template <typename T1, typename T2, size_t R, size_t C>
auto operator*(const Matrix<T1, R, C>& M, const Vector<T2, C>& V)
-> const Vector<typename std::remove_cv<decltype(operator*(M[0][0], V[0]))>::type, R>
{
    Vector<typename std::remove_cv<decltype(operator*(M[0][0], V[0]))>::type, R> v;
    for (size_t r = 0; r < R; ++r)
        for (size_t c = 0; c < C; ++c)
            v[r] += M[r][c] * V[c];
    return v;
}

template <typename T, size_t N>
const SquareMatrix<T, N> outer(const Vector<T, N>& A, const Vector<T, N>& B)
{
    SquareMatrix<T, N> m;
    for (size_t r = 0; r < N; ++r)
        for (size_t c = 0; c < N; ++c)
            m[r][c] = A[r] * B[c];
    return m;
}

}
#endif