Reduction bug (?) with sparse tensors

[ryanmrichard]

Not sure if this is a feature or a bug, but as far as I can tell the reduction functor is not used to reduce the tiles to scalars if one of the tiles is zero. For reductions without the zero-product property, this produces an incorrect result. MWE below.

class MyReducer {
public:
    // Required typedefs
    using result_type = double;
    using first_argument_type  = TA::TensorD;
    using second_argument_type = TA::TensorD;

    using result_ref       = result_type&;
    using const_result_ref = const result_type&;
    using const_arg1_ref   = const first_argument_type&;
    using const_arg2_ref   = const second_argument_type&;

    result_type operator()() const { return 0.0; }

    const_result_ref operator()(const_result_ref r) const { return r; }

    void operator()(result_ref result, const_result_ref arg) const {
        result = result + arg;
    }

    void operator()(result_ref result, const_arg1_ref first,
                    const_arg2_ref second) const {
        // Sums up the difference between the tiles using the fact that there's only one element
        result += first(0) - second(0);
    }
};

I would expect the following to print 1, but it prints 0:

TA::TiledRange trange{{0, 1}};
//LHS is vector with a single element 1
TA::TSpArrayD lhs(world, trange, {1});
//RHS is a sparse vector with no dense tiles
auto rhs = TA::make_array<TA::TSpArrayD>(world, trange,
    [](TA::TensorD& tile, const TA::Range& r){
        TA::TensorD temp(r);
        temp(0) = 0.0;
        tile = temp;
        return tile.norm();
    });
    MyReducer r;
    // Does 1 - 0 and gets 0
    auto v =lhs("i").reduce(rhs("i"), r);
    std::cout << v.get() << std::endl;

If you change temp(0) = 0.0; to temp(0) = 2.0 the code correctly prints -1.

[ryanmrichard]

Interestingly it appears to work if the this tensor for the reduction (lhs in my MWE) is the empty sparse tensor.

[evaleev]

this is not a bug, this is a missing feature ... the binary reduction op is run over a logical-and of the shapes to the two arguments. Indeed, note that in MyReducer there are no operator() overloads where the first or the second argument is "missing".

[justusc]

@evaleev is correct. I didn't encounter any cases where it wasn't necessary to have a logical-or reduction, or could not be accomplished by other means. In this case, you can:

double v = (lhs("i") - rhs("i")).sum();

The feature could be added with a small amount of work in the binary reduction expression function. I can point out the bit of code, if anyone wants to take a stab at it.

[ryanmrichard]

@justusc the subtraction was an example. What I'm literally trying to do is to implement Numpy's allclose, i.e., std::fabs(lhs) <= atol + rtol * std::fabs(rhs) followed by logical-and (lhs and rhs are elements of the tensors atol and rtol define "close") are followed by logical-and. The problem is the comparison is skipped for zero tiles, which is fine only if both tiles are supposed to be zero.

[evaleev]

@justusc thanks, I think expr.h:665-671 + minor mods elsewhere ... with c++17 it will be a relatively straightforward spaghetty of if constexpr

[justusc]

@evaleev yes, that is what I was going to point out. On a side note, c++17 would have made a lot of TA code much simpler.