investigate reduction within parallel_for

[paciorek]

for parallelized declFuns, we have need for code like this:

  nc <- nClass(
        Cpublic = list(
            go = nFunction(
                fun = function(x = 'numericVector') {
                    logProb_y <- x
                    sum <- 0
                    parallel_for(i, 1:length(x), {
                        logProb_y[i] <- 2*x[i]
                        sum <- sum + logProb_y[i]
                    })
                    return(sum)
                },
                returnType = 'numericScalar'
            )
        )
    )
    Cnc <- nCompile(nc)

Currently that gives this error in our compilation processing.

Error in TensorReduction(<environment>, NULL) : 
  could not find function "TensorReduction"

So we'll want to investigate if this is possible with TBB. Perhaps only with TBB reduction functionality given the need to write to shared sum.

[paciorek]

Claude indicates this is possible using TBB parallel_reduce. Something like this:

#include <tbb/parallel_reduce.h>
#include <tbb/blocked_range.h>
#include <vector>
#include <iostream>

int main() {
    const size_t N = 1000;
    std::vector<double> v(N);

    // Lambda-based parallel_reduce: transform + sum in one pass
    double total = tbb::parallel_reduce(
        tbb::blocked_range<size_t>(0, N),

        // Identity value (starting sum)
        0.0,

        // Per-chunk work: assign new values AND accumulate local sum
        [&v](const tbb::blocked_range<size_t>& r, double local_sum) {
            for (size_t i = r.begin(); i != r.end(); ++i) {
                v[i] = static_cast<double>(i * i);  // transform
                local_sum += v[i];                   // reduce
            }
            return local_sum;
        },

        // Combine results from different chunks
        [](double a, double b) { return a + b; }
    );

    std::cout << "Sum: " << total << std::endl;
}

[paciorek]

We'd presumably want code like this, across some combination of a calc_1_par and primitive calc_par (in place of calc_one):

            parallel_for_reduce(i, 1:instr$lens[1], {
                myidx <- instr$values[[1]][1]+i
                model$logProb_y[myidx] <- dnorm(model$y[myidx], mean = model$mu, sd = model$lifted_sqrt_oPtau_cP, log = 1)
                logProb <- logProb + model$logProb_y[myidx]
            })