[embray]

Two bugs, both the symptom and the cause. The cause was insufficient NaN comparison in ImageMode::getBackGuess--on debug builds a bare assert(!isnan(mean)) caught mean=NaN, but on release builds this fell through. Because of -ffast-math the NaN was simply allowed to propagate through the statistics.

In KappaSigmaBinning::computeBins this NaN propagation also meant that a NaN could poison the whole cell, resulting in ndata==0. This was also asserted against on debug builds, but allowed to propagate on release builds.

Uses the fastmath_isnan and fastmath_isinf from SEUtils since these were introduced exactly for this reason, that we still need to sweep for non-finite values even of -ffast-math builds which are created by default on release.

Update: This turned out to be a more serious issue affecting this code on nearly flat background cells with close adjacent values, resulting in some caes in catastrophic cancellation. I asked Claude to benchmark different alternative algorithms, namely shifted data and Welford's, against similar data that reproduced the problem, both for speed and accuracy. The benchmark code is attached for reference and results from my machine. It looks like the shifted-data approach is just as fast as the old code, and still highly accurate. So I implemented it myself.
variance-algorithms.cpp
variance-algorithms.txt

There was also a comment in the code that suggested it was faster that boost::accumulator. It turns out that actually has two algorithms for calculating the variance: accumlators::lazy_variance and accumulators::variance. That claim was probably against the "lazy" variant, whereas the non-lazy variant is essentially Weldford's.

In any case the claim still seems to hold that this is a bit faster.

---

This fix was also applied on my ASDF branch at #631, though it really should be a stand-alone fix. It has nothing to do specifically with ASDF. I first found it on ASDF coadd files from Roman, but the same bugs occur if you convert the file to FITS.

[embray]

It turns out @braingram is still tripping at this exception which I had argued "shouldn't be reachable anymore" and yet apparently is on arm64. I will have to trace the code path again and see if there are any other ways to trip on NaNs, esp. with -ffast-math. I don't think the problem is likely to due the variance calculation anymore but anything is possible I suppose. In the meantime will try to get my hands on an arm64 mac to test on.

[embray]

This situation is definitely worse than I thought, and it is both compiler and architecture specific to some extent. When release builds of SX++ are created they use -ffast-math by default, which means -ffinite-math-only. In our discussion yesterday there was some vague mention that we tend to just hope there are no NaNs, but that's clearly not a valid assumption, as evidenced by #626, which does not carry any explanation or background, but is clear from the patch that they were showing up in PSFx files.

The attempted fix of type-punning floats so ints and using bit twiddling to check for non-finite values looks quite reasonable on its face, and was found to work! In fact, it works for me too, where I was compilign with my system gcc v13.3. However, modern versions of clang recognize this pattern and actually compile it out entirely which is rather aggressive to me. But -ffinite-math-only really does mean "I promise that none of my data or calculations will ever produce NaNs or Infs". This is a bold promise to make especially if you are dealing then with external data sources over which the software has no control. And in fact, produced undefined behavior.

Take this simple example:

#include <stdio.h>
#include <stdlib.h>

int main(int argc, char **argv) {
    if (argc < 2) {
        fprintf(stderr, "usage: %s <float>\n", argv[0]);
        return 1;
    }
    float x = strtof(argv[1], NULL);     /* from outside: the compiler cannot fold it */
    int lt = x < 0.0f, eq = x == 0.0f, gt = x > 0.0f;
    printf("%s (x=%f):\n  x<0  => %d\n  x==0 => %d\n  x>0  => %d\n  sum  =  %d\n",
           argv[1], x, lt, eq, gt, lt + eq + gt);
    return 0;
}

Here are the results I get for gcc vs clang with different compiler flags. This is just on the same x86_64 machine:

test on argv[1]=NaN)	clang -O2 -ffast-math	gcc -O2 -ffastmath
x < 0	1	0
x == 0	1	1
x != 0	0	0
x > 0	0	0
sum	2	1

The contradiction is self-evident: a value cannot be both less than 0 and equal zero as in the gcc case, not even if it's not a number. Either way, the disagreement between compilers is more significant.

If there are branches on comparisons like that it can result in entirely different application behavior based on compiler alone, much less machine architecture and how these comparisons are compiled in assembly.

A concrete, real outcome of this difference showed up for me running SX++ compiled under gcc vs clang on the same input file. The gcc build took only a few seconds to run, while the clang build took over half an hour. I traced the source of this behavior to Lutz::labelImage:

  for (int y = 0; y < lines; y++) {
    LutzStatus ps = LutzStatus::COMPLETE;
    LutzStatus cs = LutzStatus::NONOBJECT;

    if (y % chunk_height == 0) {
      chunk = image.getChunk(0, y, image.getWidth(), std::min(chunk_height, lines - y));
    }

    int dy = y % chunk_height;
    for (int x=0; x < width; x++) {
      DetectionImage::PixelType value = (x == width - 1) ? 0.0 : chunk->getValue(x, dy);

      LutzMarker last_marker = marker[x];
      marker[x] = LutzMarker::ZERO;

      bool in_object = value > 0.0;
      if (in_object) {
        ...
      }
  ...
}

The crux is the value > 0.0 comparison. As we just saw if value=NaN on GCC the result is 0, while on clang it's 1! This turns out to be rather catastrophic if you have a border of all NaN values. It means each "connected" group of NaNs is merged into the previous one (PixelGroup::merge_pixel_list) which does pixel_list.insert(...), where std::vector<PixelCoordinate> pixel_list. This process repeats once per scan line as long as the pixel group remains open:

          } else {
            // Add group to current group
            auto prev_group = inc_group_map.at(x);
            inc_group_map.erase(x);

            group_stack.back().merge_pixel_list(prev_group);
          }

That case in particular didn't use std::move either resulting in repeating copies of the same PixelGroup being created, and growing with each line resulting in quadratic behavior. This is not noticeable for a typical source, but as runs of NaNs are being treated here as belonging to a connected group. If the border of the image outside the detector footprint is all NaN this can be quite enormous. In fact this is a bug in Lutz outside the -ffast-math and NaN handling: Any oversized connected source could in principle trigger this bug, I think. Though I admit I don't fully understand the algorithm enough to tell if this was intentional or not.

This also resulted in memory thrashing since this large pixel list buffer quickly grows > 32 MB (M_MMAP_THRESHOLD on my system) so malloc could not reuse it. Each copy is a new mmap/munmap, which is actually what led me to this analysis in the first place. Running the slow clang build under ptrace showed it spending it most its time on repeated large mmap/munmap calls that I eventually traced here.

The whole Lutz business is an aside from the main problem I'm analyzing here, but is a catastrophic result of the aforementioned UB surrounding NaNs.

ISTM the only truly reliable workaround is to disable this with -fno-finite-math-only, though I can't say for sure what other negative impacts that would have in areas that lack nan/inf checks.