My original Javascript function looks like this:

function fitIT(IT, Xs, ys){
    try{
        let ysHat = ys.map(y => [y]);
        let XsHat = Xs.map(function(X){
            return zip(IT.funcs, IT.terms).map(function([f, term]){
                return f(zip(X, term).map(function([x, t]){
                    return Math.pow(x, t);
                }).reduce(prod, 1));
            });
        });
    
        let XsHatT = math.transpose(XsHat);
    
        let q1 = math.inv(math.multiply(XsHatT, XsHat));
        let q2 = math.multiply(XsHatT, ysHat);
        
        let ws = math.multiply(q1, q2);

        IT.coeffs = ws.map(w => w[0]);
    }
    catch(err){ //não foi possível encontrar inversa, coeficientes ficam 1
        //console.log(err);
    }
    
    let MAE = zip(Xs, ys).map(function([X, y]){
        return Math.abs(evalIT(IT, X) - y);
    }).reduce(sum, 0) / Xs.length;

    let score = 1/(1 + MAE);

    IT.score = isFinite(score) ? score : 0.0;

    return IT;
}

I am trying to re-write it in C++ using Eigen:

void fit(const Eigen::Matrix<TestProject::Scalar, Eigen::Dynamic, 1>& X, const Eigen::Map<const Eigen::Matrix<TestProject::Scalar, Eigen::Dynamic, 1>>& y) {
    try {
        Eigen::VectorXd yHat = Eigen::Map<const Eigen::VectorXd>(y.data(), static_cast<Eigen::Index>(y.size()));
        Eigen::Matrix<TestProject::Scalar, 1, Eigen::Dynamic> XHat(Functions.size());
        std::transform(
            std::execution::par,
            Polynomials.begin(), Polynomials.end(),
            Functions.begin(),
            XHat.data(),
            [&](const Vector<Scalar>& p, auto f) {
               Scalar transformedX = std::transform_reduce(
                       std::execution::par,
                       X.begin(),
                       X.end(),
                       p.begin(),
                       1.0,
                       std::multiplies(),
                       [](double x, double p) { return std::pow(x, p); }
               );
               return f(transformedX);
            }
        );
        Eigen::MatrixXd XsHatT = XHat.transpose();

        Eigen::MatrixXd q1 = (XsHatT * XHat).inverse();
        Eigen::MatrixXd q2 = XsHatT * yHat;

        Eigen::VectorXd ws = q1 * q2;

        Coefficients.resize(ws.size());
        for (int i = 0; i < ws.size(); ++i) {
            Coefficients[i] = ws[i];
        }
    } catch (const std::exception& e) {
        std::cerr << "Matrix inversion failed: " << e.what() << std::endl;
        Coefficients.assign(Functions.size(), 1.0);
    }

    Scalar MAE = std::transform_reduce(
            X.data(), X.data() + X.size(), y.begin(), 0.0, std::plus<>(),
            [&](Scalar X, Scalar y) {
                Scalar yHat = eval({X});
                return std::abs(yHat - y);
            }
    ) / static_cast<double>(X.size());

    Scalar score = 1 / (1 + MAE);
    Fitness = std::isfinite(score) ? score : 0.0;
}

On a different class I’ve read the data from file and stored them in an Eigen matrix. Here I’m trying to split them and call the function:

Eigen::Ref<const TestProject::Dataset::Matrix> dataMatrix = dataset->Values();
Eigen::Index targetIndex = dataMatrix.cols() - 1;
Eigen::Map<const Eigen::Matrix<TestProject::Scalar, Eigen::Dynamic, 1>> ysHat(dataMatrix.col(targetIndex).data(), dataMatrix.rows());
TestProject::Dataset::Matrix XsHat(dataMatrix.rows(), dataMatrix.cols() - 1);
XsHat << dataMatrix.leftCols(targetIndex), dataMatrix.rightCols(dataMatrix.cols() - targetIndex - 1);

newInd.fit(XsHat, ysHat);

This throws an error: void Eigen::PlainObjectBase<Derived>::resizeLike(const Eigen::EigenBase<OtherDerived>&) [with OtherDerived = Eigen::Array<double, -1, -1>; Derived = Eigen::Matrix<double, -1, 1>]: Assertion 'other.rows() == 1 || other.cols() == 1' failed.

If I print ysHat I get a 500×1 vector, and XsHat a 500×5 matrix.