calcBandWidth scans only +l translations -> asymmetric convolution error for mu>0

[guycohen]

Summary

OperatorTree::calcBandWidth only inspects l>=0 (via getNode(depth, l)) and ignores the symmetric -l. This makes the computed bandwidth one‑sided. isOutsideBand then prunes operator terms on the opposite side, which shows up as a large error for positive shifts (e.g., mu>0), while mu<0 looks fine.

Minimal reproducer

#include "MRCPP/MWFunctions"
#include "MRCPP/MWOperators"
#include "MRCPP/Printer"
#include <array>
#include <cmath>
#include <iostream>

int main() {
    mrcpp::Printer::init(0);

    const double t_max = 10.0;
    const double alpha0 = 1.0;
    const double alpha = 1.0;
    const double prec = 1.0e-5;

    for (double mu : {5.0, -5.0}) {
        std::cout << "mu = " << mu << "\n";
        mrcpp::BoundingBox<1> bbox(0, {-1}, {2}, {t_max}, false);
        mrcpp::MultiResolutionAnalysis<1> mra(bbox, 7);

        auto delta = [alpha0](const mrcpp::Coord<1> &r) {
            return std::exp(-alpha0 * r[0] * r[0]);
        };

        mrcpp::FunctionTree<1> delta_tree(mra);
        mrcpp::project<1, double>(prec, delta_tree, delta);

        mrcpp::GaussExp<1> kernel;
        kernel.append(mrcpp::GaussFunc<1>(alpha, 1.0, {mu}));

        mrcpp::ConvolutionOperator<1> op(mra, kernel, prec);
        mrcpp::FunctionTree<1> out(mra);
        mrcpp::apply(prec, out, op, delta_tree);

        double best_t = 0.0, best_val = -1.0;
        for (int i = 0; i <= 400; ++i) {
            double t = -10.0 + 0.05 * i;
            double v = out.evalf({t});
            if (v > best_val) { best_val = v; best_t = t; }
        }
        std::cout << "peak at t=" << best_t << "\n";
    }
}

Observed (before patch)

mu = 5.000000e+00
peak at t=5.500000000000e+00
mu = -5.000000000000e+00
peak at t=-5.000000000000e+00

Observed (after patch)

This asymmetry disappears if bandwidth is computed from both +l and -l:
mu = 5.000000e+00
peak at t=5.000000000000e+00
mu = -5.000000000000e+00
peak at t=-5.000000000000e+00

Patch

  diff --git a/src/trees/OperatorTree.cpp b/src/trees/OperatorTree.cpp
  index XXXXXXX..YYYYYYY 100644
  --- a/src/trees/OperatorTree.cpp
  +++ b/src/trees/OperatorTree.cpp
  @@ -119,11 +119,13 @@ void OperatorTree::calcBandWidth(double prec) {
           int l = 0;
           bool done = false;
           while (not done) {
               done = true;
  -            MWNode<2> &node = getNode(depth, l);
  +            MWNode<2> &node_pos = getNode(depth, l);
  +            MWNode<2> &node_neg = (l == 0) ? node_pos : getNode(depth, -l);
               double thrs = std::max(MachinePrec, prec / (8.0 * (1 << depth)));
               for (int k = 0; k < 4; k++) {
  -                if (node.getComponentNorm(k) > thrs) {
  +                double norm_pos = node_pos.getComponentNorm(k);
  +                double norm_neg = node_neg.getComponentNorm(k);
  +                if (std::max(norm_pos, norm_neg) > thrs) {
                       this->bandWidth->setWidth(depth, k, l);
                       done = false;
                   }
               }