eigen_util.h

Go to the documentation of this file.

/* 
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Daniel Abele
*
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
*     http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#ifndef MIO_MATH_EIGEN_UTIL_H
#define MIO_MATH_EIGEN_UTIL_H
 
#include "memilio/math/eigen.h" // IWYU pragma: keep
 
#include <utility>
#include <iterator>
 
namespace mio
{
 
template <typename T>
struct Seq {
    Seq(T start_, T n_, T stride_ = 1)
        : start(start_)
        , n(n_)
        , stride(stride_)
    {
        assert(start >= 0);
        assert(n >= 0);
        assert(stride >= 1);
    }
    T start, n, stride = 1;
};
 
template <class M>
concept IsDynamicVector = (std::remove_reference_t<M>::RowsAtCompileTime == Eigen::Dynamic &&
                           std::remove_reference_t<M>::ColsAtCompileTime == 1) ||
                          (std::remove_reference_t<M>::RowsAtCompileTime == 1 &&
                           std::remove_reference_t<M>::ColsAtCompileTime == Eigen::Dynamic);
 
template <class M>
concept IsDynamicMatrix = std::remove_reference_t<M>::RowsAtCompileTime == Eigen::Dynamic &&
                          std::remove_reference_t<M>::ColsAtCompileTime == Eigen::Dynamic;
 
template <class M>
Eigen::Index major_size(M&& m)
{
    return std::remove_reference_t<M>::IsRowMajor ? m.rows() : m.cols();
}
 
template <class M>
Eigen::Index minor_size(M&& m)
{
    return std::remove_reference_t<M>::IsRowMajor ? m.cols() : m.rows();
}
 
template <class M>
struct CVPlainMatrix {
    using Type = typename M::PlainMatrix;
};
template <class M>
struct CVPlainMatrix<Eigen::Ref<const M>> {
    using Type = const M;
};
template <class M>
using CVPlainMatrixT = typename CVPlainMatrix<M>::Type;
 
template <IsDynamicVector V>
auto slice(V&& v, Seq<Eigen::Index> elems)
{
    return Eigen::Map<CVPlainMatrixT<std::decay_t<V>>, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>(
        v.data() + elems.start, elems.n, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>{1, elems.stride});
}
 
template <IsDynamicMatrix M>
auto slice(M&& m, Seq<Eigen::Index> rows, Seq<Eigen::Index> cols)
{
    assert(rows.start + rows.stride * rows.n <= m.rows());
    assert(cols.start + cols.stride * cols.n <= m.cols());
 
    auto majSpec   = std::remove_reference_t<M>::IsRowMajor ? rows : cols;
    auto minSpec   = std::remove_reference_t<M>::IsRowMajor ? cols : rows;
    auto majStride = majSpec.stride * minor_size(m);
    auto minStride = minSpec.stride;
    auto data      = m.data() + majSpec.start * minor_size(m) + minSpec.start;
 
    return Eigen::Map<CVPlainMatrixT<std::remove_reference_t<M>>, Eigen::Unaligned,
                      Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>(data, rows.n, cols.n, {majStride, minStride});
}
 
template <typename M>
auto reshape(M&& m, Eigen::Index rows, Eigen::Index cols)
{
    assert(rows * cols == m.rows() * m.cols());
    assert(rows >= 1);
    assert(cols >= 1);
 
    return Eigen::Map<CVPlainMatrixT<std::remove_reference_t<M>>>(m.data(), rows, cols);
}
 
template <class M>
concept IsMatrixExpression = std::is_base_of_v<Eigen::EigenBase<M>, M>;
 
template <class A, class B>
auto max(const Eigen::MatrixBase<A>& a, B&& b)
{
    using std::max;
    return a.binaryExpr(std::forward<B>(b), [](auto a_i, auto b_i) {
        return max(a_i, b_i);
    });
}
 
template <class Rng, class F>
auto map(const Rng& v, F f)
{
    using Result = std::invoke_result_t<F, const typename Rng::value_type&>;
    using Scalar = std::decay_t<Result>;
    using Array  = Eigen::Array<Scalar, Eigen::Dynamic, 1>;
    return Array::NullaryExpr(v.size(), [f, &v](Eigen::Index i) -> Scalar {
        return f(v[size_t(i)]);
    });
}
 
namespace details
{
//true if elements returned by matrix(i, j) are references where matrix is of type M;
//false if the elements are temporaries, e.g. for expressions like Eigen::MatrixXd::Constant(r, c, v).
template <class M>
concept IsElementReference =
    std::is_reference_v<decltype(std::declval<M>()(std::declval<Eigen::Index>(), std::declval<Eigen::Index>()))>;
} // namespace details
 
template <class M, bool IsConst = false>
class RowMajorIterator
{
public:
    using MatrixRef = std::conditional_t<IsConst, const M&, M&>;
    using MatrixPtr = std::conditional_t<IsConst, const M*, M*>;
    static_assert(IsConst || details::IsElementReference<M>, "Iterator must be const if matrix is not in memory.");
 
    using iterator_category = std::random_access_iterator_tag;
    using value_type        = typename M::Scalar;
    using reference         = std::conditional_t<IsConst, const value_type&, value_type&>;
    using difference_type   = Eigen::Index;
 
    // Operator-> returns a pointer to an element if the element can be referenced.
    // If the matrix returns temporaries, a proxy is returned that holds a local copy
    // and forwards the address of that copy. Adress is pointer-to-const because modifying
    // the local inside the proxy is almost certainly not what is intended.
    struct Proxy {
        value_type value;
        const value_type* operator->() const
        {
            return &value;
        }
    };
    using pointer = std::conditional_t<details::IsElementReference<MatrixRef>,
                                       std::conditional_t<IsConst, const value_type*, value_type*>, Proxy>;
 
    RowMajorIterator(MatrixRef m, Eigen::Index i)
        : m_matrix(&m)
        , m_i(i)
    {
    }
 
    explicit RowMajorIterator()
        : m_matrix(nullptr)
        , m_i(0)
    {
    }
 
    RowMajorIterator& operator++()
    {
        ++m_i;
        return *this;
    }
    RowMajorIterator operator++(int)
    {
        auto cpy = *this;
        ++m_i;
        return cpy;
    }
    RowMajorIterator operator+(difference_type n) const
    {
        return {*m_matrix, m_i + n};
    }
    friend RowMajorIterator operator+(difference_type n, const RowMajorIterator& iter)
    {
        return {*iter.m_matrix, iter.m_i + n};
    }
    RowMajorIterator& operator+=(difference_type n)
    {
        m_i += n;
        return *this;
    }
    RowMajorIterator& operator--()
    {
        --m_i;
        return *this;
    }
    RowMajorIterator operator--(int)
    {
        auto cpy = *this;
        --m_i;
        return cpy;
    }
    RowMajorIterator operator-(difference_type n) const
    {
        return {*m_matrix, m_i - n};
    }
    RowMajorIterator& operator-=(difference_type n)
    {
        m_i -= n;
        return *this;
    }
    difference_type operator-(RowMajorIterator r) const
    {
        return m_i - r.m_i;
    }
 
    decltype(auto) operator*() const
    {
        return (*m_matrix)(m_i / m_matrix->cols(), m_i % m_matrix->cols());
    }
 
    pointer operator->() const
        requires details::IsElementReference<MatrixRef>
    {
        return &(**this);
    }
 
    pointer operator->() const
        requires(!details::IsElementReference<MatrixRef>)
    {
        return Proxy{**this};
    }
    bool operator<(const RowMajorIterator& other) const
    {
        return m_i < other.m_i;
    }
    bool operator>(const RowMajorIterator& other) const
    {
        return other < *this;
    }
    bool operator<=(const RowMajorIterator& other) const
    {
        return !(other < *this);
    }
    bool operator>=(const RowMajorIterator& other) const
    {
        return !(*this < other);
    }
    bool operator==(const RowMajorIterator& other) const
    {
        return !(*this < other) && !(other < *this);
    }
    bool operator!=(const RowMajorIterator& other) const
    {
        return !(*this == other);
    }
 
private:
    MatrixPtr m_matrix;
    Eigen::Index m_i;
};
 
template <IsMatrixExpression M>
    requires details::IsElementReference<M>
RowMajorIterator<M, false> begin(M& m)
{
    return {m, 0};
}
 
template <IsMatrixExpression M>
RowMajorIterator<M, true> begin(const M& m)
{
    return {m, 0};
}
 
template <IsMatrixExpression M>
RowMajorIterator<M, true> cbegin(const M& m)
{
    return {m, 0};
}
 
template <IsMatrixExpression M>
    requires details::IsElementReference<M>
RowMajorIterator<M, false> end(M& m)
{
    return {m, m.size()};
}
 
template <IsMatrixExpression M>
RowMajorIterator<M, true> end(const M& m)
{
    return {m, m.size()};
}
 
template <IsMatrixExpression M>
RowMajorIterator<M, true> cend(const M& m)
{
    return {m, m.size()};
}
 
} // namespace mio
 
#endif // MIO_MATH_EIGEN_UTIL_H