Блог

У меня есть несколько примеров в моем коде, где у меня есть условие, основанное на коэффициентах массивов 1xN, и мне нужно задать целые столбцы массивов MxN в зависимости от этих условий. В моем случае N равно Eigen::Dynamic , а M колеблется от 2 до 4, но в каждом случае является константой времени компиляции.

Вот простая функция, иллюстрирующая, что я имею в виду, с a массивами 1xN и b являющимися массивами 1xN, которые формируют условие, c являющимися массивом 2xN с дополнительными данными и res являющимися выходным параметром, столбцы которого всегда задаются как единое целое:

 #include <iostream>
#include <Eigen/Dense>

using namespace Eigen;

template<Index nRows>
using ArrayNXd = Array<double, nRows, Dynamic>;

using Array1Xd = ArrayNXd<1>;
using Array2Xd = ArrayNXd<2>;
using Array3Xd = ArrayNXd<3>;

void asFunction(
    Array3Xdamp; res,
    const Array1Xdamp; a, const Array1Xdamp; b, const Array2Xdamp; c
){
    for (Index col{0}; col<a.cols();   col){
        if ( a[col] > b[col] )
            res.col(col) = Array3d{
                 a[col]   b[col],
                (a[col]   b[col]) * c(0, col),
                (a[col] - b[col]) * c(1, col)
            };
        else
            res.col(col) = Array3d{
                 a[col] - b[col],
                 a[col]   b[col],
                (a[col]   b[col]) * (a[col] - b[col])
            };
    }
}


int main(){
    Array1Xd a (3), b(3);
    Array2Xd c (2, 3);
    
    a << 1, 2, 3;
    b << 0, 1, 2;
    c <<
        0, 1, 2,
        1, 2, 3;

    Array3Xd res (3,3);
    
    asFunction(res, a, b, c);

    std::cout << "as function:n" << res << "n";

    return 0;
}
 

Функции, подобные этой, используются в разделе, посвященном производительности* моего кода, и я чувствую, что оставляю производительность на столе, потому что использование циклов с Eigen типами обычно не является оптимальным решением.

*да, я составил его профиль.

Я написал ту же функцию , что и a NullaryExpr , но это было немного медленнее. Я думаю, что это имеет смысл, учитывая дополнительные оценки условий и ветвления для каждой строки:

 #include <iostream>
#include <Eigen/Dense>

using namespace Eigen;

template<Index nRows>
using ArrayNXd = Array<double, nRows, Dynamic>;

using Array1Xd = ArrayNXd<1>;
using Array2Xd = ArrayNXd<2>;
using Array3Xd = ArrayNXd<3>;

class MyFunctor
{
public:
    using Scalar = double;

    static constexpr Index
        RowsAtCompileTime { 3 },
        MaxRowsAtCompileTime { 3 },
        ColsAtCompileTime { Dynamic },
        MaxColsAtCompileTime { Dynamic };

    using DenseType = Array<
        Scalar  ,    RowsAtCompileTime,    ColsAtCompileTime,
        ColMajor, MaxRowsAtCompileTime, MaxColsAtCompileTime
    >;

private:
    typename Array1Xd::Nested m_a;
    typename Array1Xd::Nested m_b;
    typename Array2Xd::Nested m_c;

public:
    MyFunctor(
        const Array1Xdamp; a,
        const Array1Xdamp; b,
        const Array2Xdamp; c
    ) : m_a {a}, m_b {b}, m_c{c}
    {}

    bool cond(Index col) const {
        return m_a[col] > m_b[col];
    }

    Scalar func1(Index col) const {
        return m_a[col]   m_b[col];
    }

    Scalar func2(Index col) const {
        return m_a[col] - m_b[col];
    }

    Scalar func3(Index row, Index col) const {
        switch(row){
            case 0: return func1(col);
            case 1: return func1(col) * m_c(0, col);
            case 2: return func2(col) * m_c(1, col);
            default: __builtin_unreachable();
        }
    }

    Scalar func4(Index row, Index col) const {
        switch (row){
            case 0: return func2(col);
            case 1: return func1(col);
            case 2: return func1(col) / func2(col);
            default: __builtin_unreachable();
        }
    }

    Scalar operator() (Index row, Index col) const {
        if ( cond(col) )
            return func3(row, col);
        else
            return func4(row, col);
    }
};

using MyReturnType = Eigen::CwiseNullaryOp<
    MyFunctor, typename MyFunctor::DenseType
>;


MyReturnType asFunctor(
    const Array1Xdamp; a,
    const Array1Xdamp; b,
    const Array2Xdamp; c
){
    using DenseType = typename MyFunctor::DenseType;
    return DenseType::NullaryExpr(
        3, a.cols(),
        MyFunctor(a, b, c)
    );
}


int main(){
    Array1Xd a (3), b(3);
    Array2Xd c (2, 3);
    
    a << 1, 2, 3;
    b << 0, 1, 2;
    c <<
        0, 1, 2,
        1, 2, 3;

    std::cout << "as functor:n" << asFunctor(a,b,c) << "n";

    return 0;
}
 

Мой вопрос: существует ли более эффективный способ реализации логики, аналогичной приведенной выше (оценить скалярное условие для каждого столбца матрицы, возвращаемые значения для всего столбца на основе условия) с использованием eigen библиотеки?

Примечание: использование выражения было бы несколько предпочтительнее, потому что мне не нужно беспокоиться о выделении памяти, внешних параметрах и т. Д., А код может быть написан с учетом скаляров, что делает его гораздо более понятным.

Редактировать: Примечание 2: Я тоже пытался использовать <Condition>.template replicate<nRows,1>().select(..., ...) , но это было медленнее и труднее читать.

Вы можете использовать метод выбора Eigen, но он работает только для скаляров, поэтому вам нужно пройти по одному измерению.

 const auto condition = a > b;
res.row(0) = condition.select(a   b /*true*/, a - b /*false*/);
res.row(1) = condition.select((a   b) * c.row(0), a   b);
res.row(2) = condition.select((a - b) * c.row(1), (a   b) * (a - b));
 

Обратите внимание, что вы, вероятно, будете быстрее, если перенесете все свои массивы. Затем итерация идет столбец за столбцом, которые векторизованы намного лучше, так как Eigen является основным столбцом.

поэтому я посмотрел только на этот фрагмент кода

     for (Index col{0}; col<a.cols();   col){
        if ( a[col] > b[col] )
            res.col(col) = Array3d{
                 a[col]   b[col],
                (a[col]   b[col]) * c(0, col),
                (a[col] - b[col]) * c(1, col)
            };
        else
            res.col(col) = Array3d{
                 a[col] - b[col],
                 a[col]   b[col],
                (a[col]   b[col]) * (a[col] - b[col])
            };
    }
 

Я подозреваю, но не могу доказать, что эти a[col] и b[col] получают доступ каждый раз, когда вы им звоните. Возможно, вам захочется попробовать создать короткие временные интервалы для значений, которые вы используете повторно. Например:
поэтому я посмотрел только на этот фрагмент кода

     for (Index col{0}; col<a.cols();   col){
        auto acol=a[col];
        auto bcol=b[col];
        auto apb=acol bcol;
        auto amb=acol-bcol;
        if ( acol > bcol )
            res.col(col) = Array3d{
                 apb,
                (apb) * c(0, col),
                (amb) * c(1, col)
            };
        else
            res.col(col) = Array3d{
                 amb,
                 apb,
                (apb) * (amb)
            };
    }
 

и да, я знаю, что это не совсем то, чего ты хотел. может быть, это поможет тебе