In [1]:

    

#include <eigen3/Eigen/Core>
#include <iostream>
#include <tuple>


    

In [2]:

    

using EigenArray2d = Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
using EigenArray1d = Eigen::Array<double, Eigen::Dynamic, 1>;


    

In [3]:

    

constexpr auto ni = 3;
constexpr auto nj = 2;

auto vec_ni = EigenArray1d(ni);
vec_ni << 1.0, 2.0, 3.0;

auto vec_nj = EigenArray1d(nj);
vec_nj << 40.0, 50.0;

auto mat_ni_nj = EigenArray2d(ni, nj);
mat_ni_nj <<
    600.0, 700.0,
    800.0, 900.0,
    1000.0, 1100.0;

std::cout << "vec_ni:\n" << vec_ni << "\n\n";
std::cout << "vec_nj:\n" << vec_nj << "\n\n";
std::cout << "mat_ni_nj:\n" << mat_ni_nj << "\n\n";


    

    




vec_ni:
1
2
3

vec_nj:
40
50

mat_ni_nj:
 600  700
 800  900
1000 1100

In [4]:

    

std::cout << "vec_nj * mat_ni_nj:\n"
    << mat_ni_nj.rowwise() * vec_nj.transpose() << "\n\n";


    

    




vec_nj * mat_ni_nj:
24000 35000
32000 45000
40000 55000

In [5]:

    

std::cout << "vec_ni * mat_ni_nj.T:\n"
    << (mat_ni_nj.transpose().rowwise() * vec_ni.transpose()) << "\n\n";


    

    




vec_ni * mat_ni_nj.T:
 600 1600 3000
 700 1800 3300

In [6]:

    

std::cout << "np.sum(vec_nj * mat_ni_nj, axis=1):\n"
    << (mat_ni_nj.rowwise() * vec_nj.transpose()).rowwise().sum() << "\n\n";


    

    




np.sum(vec_nj * mat_ni_nj, axis=1):
59000
77000
95000

In [7]:

    

auto mat_nj_ni = (mat_ni_nj.transpose().rowwise() * vec_ni.transpose()).eval();


    

In [8]:

    

std::cout << "mat_nj_ni / vec_ni:\n"
    << mat_nj_ni.rowwise() / vec_ni.transpose() << "\n\n";


    

    




mat_nj_ni / vec_ni:
 600  800 1000
 700  900 1100

In [9]:

    

// Replacing `mat_nj_ni` with it's expression

std::cout << "(vec_ni * mat_ni_nj.T) / vec_ni:\n"
    << (mat_ni_nj.transpose().rowwise() * vec_ni.transpose()).rowwise() / vec_ni.transpose() << "\n\n";


    

    




(vec_ni * mat_ni_nj.T) / vec_ni:
 600  800 1000
 700  900 1100

In [10]:

    

// Replacing rightmost `vec_ni` by the sum expression

std::cout << "(vec_ni * mat_ni_nj.T) / np.sum(vec_nj * mat_ni_nj, axis=1):\n"
    << (mat_ni_nj.transpose().rowwise() * vec_ni.transpose()).rowwise() / (mat_ni_nj.rowwise() * vec_nj.transpose()).rowwise().sum().transpose() << "\n\n";


    

    




(vec_ni * mat_ni_nj.T) / np.sum(vec_nj * mat_ni_nj, axis=1):
0.0101695 0.0207792 0.0315789
0.0118644 0.0233766 0.0347368