Different values for Eigen::VectorXd (autocasted to np.ndarray in Python) after at least two function calls

[MartinRJDagleish]

Hey,

I am currently trying to write a few tests and to speed-up a few Python functions. As I do not need a full convolution and I tend to cut away a large portion, I wanted to test, if I could use my OpenMP "knowledge" to write a fast parallel version that would instead do less computation and return the array with the correct size already.

This has multiple problems:

1. Depending on the order (and amount) of functions calls, the values seem to change?
2. Weird artifacts happen sometimes (part of array; around 5-10 values; being NaN)
3. Using const does not help
4. views do not change or help in this (tried that)
5. The problem was even worse, when I had the same name for the Eigen::VectorXd in both the serial and the parallel function version.

I tried to check if different objects in memory were created, which "seemed" to happen (see test_script.py).

My C++ code extracted:

Eigen::VectorXd
ser_convolve_cut(const nb::ndarray<double, nb::ndim<1>> &irf_vals,
                 const nb::ndarray<double, nb::ndim<1>> &model_eval) {
  Eigen::VectorXd ser_convlv_res(model_eval.size());

  size_t model_eval_size = model_eval.size();
  for (size_t i = 0; i < model_eval_size; ++i) { // outer loop for result
    for (size_t k = 0; k < i; k++) { // inner loop for convolution
      ser_convlv_res(i) += model_eval(i - k) * irf_vals(k);
    }
  }

  return ser_convlv_res;
}

Eigen::VectorXd par_convolve_cut(
    const nb::ndarray<double, nb::ndim<1>, nb::c_contig, nb::device::cpu>
        &irf_vals,
    const nb::ndarray<double, nb::ndim<1>, nb::c_contig, nb::device::cpu>
        &model_eval) {
  Eigen::VectorXd par_convlv_res(model_eval.size());
 //  auto irf_view = irf_vals.view();
 // auto model_view = model_eval.view();

  size_t model_eval_size = model_eval.size();
  size_t i, k;
#pragma omp parallel for schedule(dynamic) shared(par_convlv_res)              \
    firstprivate(irf_vals, model_eval) private(i, k)
  for (i = 0; i < model_eval_size; ++i) { // outer loop for result
    //   size_t start_k = i >= model_eval_size ? i - model_eval_size + 1 : 0;
    // size_t end_k = i < model_eval_size ? i : model_eval_size - 1;
    for (k = 0; k < i; k++) { // inner loop for convolution
      par_convlv_res(i) += model_eval(i - k) * irf_vals(k);
    }
  }

  return par_convlv_res;
}

...

NB_MODULE(nanobind_test, m) {
  m.def("ser_convolve_cut", &ser_convolve_cut,
        "Convolve two 1D arrays serially");

  m.def("par_convolve_cut", &par_convolve_cut,
        "Convolve two 1D arrays using OpenMP parallel for loop");
}

test_script.py:

import numpy as np
import matplotlib.pyplot as plt
import nanobind_test as nb

irf = np.array(...) # values
fitted_func = np.array(...) # values

# THE ORDER MATTERS A LOT HERE! 
par_res = nb.par_convolve_cut(irf, fitted_func)
ser_res = nb.ser_convolve_cut(irf, fitted_func)
ref_res = np.convolve(irf, fitted_func)

# print(ser_res.data)
# ser_res = nb.ser_convolve_cut(irf, fitted_func)
# print(ser_res.data)
# ser_res = nb.ser_convolve_cut(irf, fitted_func)
# print(ser_res.data)
# ser_res = nb.ser_convolve_cut(irf, fitted_func)
# print(ser_res.data)
# par_res = nb.par_convolve_cut(irf, fitted_func)
# ser_res = nb.ser_convolve_cut(irf, fitted_func)
print("serial has nan:",np.isnan(ser_res).any())
# print("parallel has nan:",np.isnan(par_res).any())

plt.plot(ref_res)
plt.plot(ser_res)
plt.plot(par_res)
plt.legend(["REF", "SER", "PAR"])
plt.show()

Plot for reference:
Order: par, ser, ref

Order: ser, par, ref

Order: ref, ser, par

Order: ref, ser, par, ser, par (more computation, make the results grow; although the previous should be overwritten)

[wjakob]

Why use nd-arrays for some arguments and Eigen arrays for others? For the serial version, you aren't specifying any contiguity constraints—have you considered that the arrays might not be dense in memory? For the parallel version, have you tested correctness of the parallelization outside of nanobind?

[MartinRJDagleish]

I used nd-arrays for the the arguments and Eigen arrays for the result, because that was my initial understanding after reading the docs.

I thought that the nd-array type was mainly for the interop. between numpy (Python) and th nd-array (C++), when you accept parameters and then when you want to "do some computation in C++" you should use Eigen, because those can be "autoconverted" back to numpy arrays.

I had trouble finding the code on how to initialize an C++ array with nd-array, so I used Eigen for that.

I am now retrying, if the contiguity does help or not, but from my tests I thought that were not difference there.

The parallel version should be correct and does produce the correct result, but as mentioned it is dependent on the call order.