NaN Handling

src/odc/Comparator.cc

@@ -8,7 +8,6 @@
  * does it submit to any jurisdiction.
  */
 
-#include "eckit/config/Resource.h"
 #include "eckit/exception/Exceptions.h"
 #include "eckit/filesystem/PathName.h"
 #include "eckit/log/Log.h"
@@ -41,8 +40,7 @@ class ValuesDifferent : public Exception {
 
 namespace odc {
 
-Comparator::Comparator(bool skipTestingHaveMissing) :
-    skipTestingHaveMissing_(skipTestingHaveMissing), nRow_(0), NaN_isOK_(Resource<bool>("$odc_NAN_IS_OK", false)) {}
+Comparator::Comparator(bool skipTestingHaveMissing) : skipTestingHaveMissing_(skipTestingHaveMissing), nRow_(0) {}
 
 
 void Comparator::compare(const PathName& p1, const PathName& p2) {
@@ -168,12 +166,11 @@ void Comparator::compare(int nCols, const double* data1, const double* data2, co
                     case INTEGER:
                     case BITFIELD:
                     case DOUBLE:
-                        if (!(same(*pdata1, *pdata2) || (NaN_isOK_ && (::isnan(*pdata1) && ::isnan(*pdata2)))))
+                        if (!same(*pdata1, *pdata2))
                             raiseNotEqual(column, *pdata1, *pdata2);
                         break;
                     case REAL:
-                        if (!(same(float(*pdata1), float(*pdata2)) ||
-                              (NaN_isOK_ && (::isnan(*pdata1) && ::isnan(*pdata2)))))
+                        if (!same(float(*pdata1), float(*pdata2)))
                             raiseNotEqual(column, *pdata1, *pdata2);
                         break;
                     case IGNORE:

src/odc/Comparator.h

@@ -77,15 +77,16 @@ class Comparator {
         return relativeError;
     }
 
-    inline static int same(double A, double B) { return err(A, B) < maxRelativeError; }
+    inline static bool same(double A, double B) {
+        return err(A, B) < maxRelativeError || (std::isnan(A) && std::isnan(B));
+    }
 
     void raiseNotEqual(const core::Column&, double, double);
 
 private:
 
     bool skipTestingHaveMissing_;
     long nRow_;
-    bool NaN_isOK_;
 };
 
 template <typename T1, typename T2>

src/odc/codec/Integer.h

@@ -11,6 +11,11 @@
 #ifndef odc_core_codec_Integer_H
 #define odc_core_codec_Integer_H
 
+#include <cmath>
+
+#include "eckit/exception/Exceptions.h"
+
+#include "odc/ODBAPISettings.h"
 #include "odc/core/Codec.h"
 
 /// @note We have some strange behaviour in here. In particular, we support BOTH decoding
@@ -36,7 +41,9 @@ class BaseCodecInteger : public core::DataStreamCodec<ByteOrder> {
 public:  // methods
 
     BaseCodecInteger(api::ColumnType type, const std::string& name, double minmaxmissing = odc::MDI::integerMDI()) :
-        core::DataStreamCodec<ByteOrder>(name, type), castedMissingValue_(static_cast<ValueType>(minmaxmissing)) {
+        core::DataStreamCodec<ByteOrder>(name, type),
+        castedMissingValue_(static_cast<ValueType>(minmaxmissing)),
+        integersAsDoubles_(ODBAPISettings::instance().integersAsDoubles()) {
 
         this->min_          = minmaxmissing;
         this->max_          = minmaxmissing;
@@ -65,6 +72,9 @@ class BaseCodecInteger : public core::DataStreamCodec<ByteOrder> {
 
     void gatherStats(const double& v) override {
         static_assert(sizeof(ValueType) == sizeof(v), "unsafe casting check");
+        if (integersAsDoubles_ && std::isnan(v)) {
+            throw eckit::UserError("NaN is not a valid value in INTEGER/BITFIELD column '" + this->name() + "'");
+        }
         const ValueType& val(reinterpret_cast<const ValueType&>(v));
         core::Codec::gatherStats(val);
     }
@@ -75,6 +85,8 @@ class BaseCodecInteger : public core::DataStreamCodec<ByteOrder> {
     ///         directly where needed is to work around a Cray 8.7 compiler bug, where
     ///         where the punned version gets optimised out
     ValueType castedMissingValue_;
+
+    const bool integersAsDoubles_;
 };
 
 

src/odc/core/Codec.cc

@@ -10,6 +10,8 @@
 
 #include "odc/core/Codec.h"
 
+#include <cmath>
+
 #include "eckit/exception/Exceptions.h"
 
 #include "odc/core/CodecFactory.h"
@@ -98,6 +100,8 @@ void Codec::missingValue(double v) {
 }
 
 void Codec::gatherStats(const double& v) {
+    if (std::isnan(v))
+        return;
     if (v == missingValue_) {
         hasMissing_ = 1;
     }

tests/core/CMakeLists.txt

@@ -25,6 +25,7 @@ list( APPEND _core_odc_tests
     test_text_reader
     test_table_iterator
     test_initial_missing
+    test_comparator
 )
 
 foreach( _test ${_core_odc_tests} )

tests/core/test_codecs_end_to_end.cc

@@ -8,13 +8,19 @@
  * does it submit to any jurisdiction.
  */
 
+#include <cmath>
+#include <limits>
+
 #include "eckit/io/Buffer.h"
 #include "eckit/io/DataHandle.h"
 #include "eckit/io/MemoryHandle.h"
 #include "eckit/testing/Test.h"
 
+#include "odc/MDI.h"
 #include "odc/Reader.h"
 #include "odc/Writer.h"
+#include "odc/api/ColumnType.h"
+#include "odc/api/Odb.h"
 #include "odc/codec/Integer.h"
 #include "odc/codec/String.h"
 #include "odc/core/MetaData.h"
@@ -304,6 +310,121 @@ CASE("Missing values are encoded and decoded correctly") {
 }
 
 
+CASE("NaN before real values preserves correct min/max in column metadata") {
+
+    const double NaN = std::numeric_limits<double>::quiet_NaN();
+
+    eckit::Buffer buf(4096);
+    eckit::MemoryHandle writeDH(buf);
+
+    {
+        odc::Writer<> oda(writeDH);
+        odc::Writer<>::iterator writer = oda.begin();
+
+        writer->setNumberOfColumns(1);
+        writer->setColumn(0, "col0", odc::api::DOUBLE);
+        writer->writeHeader();
+
+        (*writer)[0] = NaN;
+        ++writer;
+        (*writer)[0] = 1.1;
+        ++writer;
+        (*writer)[0] = 3.3;
+        ++writer;
+    }
+
+    eckit::MemoryHandle dh(buf.data(), static_cast<size_t>(writeDH.position()));
+    dh.openForRead();
+    odc::Reader oda(dh);
+    odc::Reader::iterator it = oda.begin();
+
+    EXPECT_EQUAL(it->columns()[0]->min(), 1.1);
+    EXPECT_EQUAL(it->columns()[0]->max(), 3.3);
+}
+
+
+CASE("NaN round-trips as NaN in REAL/DOUBLE columns") {
+
+    const double NaN = std::numeric_limits<double>::quiet_NaN();
+
+    for (auto type : {odc::api::DOUBLE, odc::api::REAL}) {
+
+        eckit::Buffer buf(4096);
+        eckit::MemoryHandle writeDH(buf);
+
+        {
+            odc::Writer<> oda(writeDH);
+            odc::Writer<>::iterator writer = oda.begin();
+
+            writer->setNumberOfColumns(1);
+            writer->setColumn(0, "col0", type);
+            writer->writeHeader();
+
+            (*writer)[0] = 1.0;
+            ++writer;
+            (*writer)[0] = NaN;
+            ++writer;
+            (*writer)[0] = 3.0;
+            ++writer;
+        }
+
+        eckit::MemoryHandle dh(buf.data(), static_cast<size_t>(writeDH.position()));
+        dh.openForRead();
+        odc::Reader oda(dh);
+
+        odc::Reader::iterator it = oda.begin();
+        EXPECT_EQUAL((*it)[0], 1.0);
+        ++it;
+        EXPECT(std::isnan((*it)[0]));
+        ++it;
+        EXPECT_EQUAL((*it)[0], 3.0);
+        ++it;
+        EXPECT_EQUAL(it, oda.end());
+    }
+}
+
+CASE("NaN is rejected in INTEGER/BITFIELD columns") {
+
+    const double NaN = std::numeric_limits<double>::quiet_NaN();
+
+    for (auto type : {odc::api::INTEGER, odc::api::BITFIELD}) {
+
+        eckit::Buffer buf(4096);
+        eckit::MemoryHandle writeDH(buf);
+
+        odc::Writer<> oda(writeDH);
+        odc::Writer<>::iterator writer = oda.begin();
+
+        writer->setNumberOfColumns(1);
+        writer->setColumn(0, "col0", type);
+        writer->writeHeader();
+
+        (*writer)[0] = 1.0;
+        ++writer;
+        (*writer)[0] = NaN;
+        EXPECT_THROWS_AS(++writer, eckit::UserError);
+    }
+}
+
+
+CASE("api::encode rejects NaN in INTEGER/BITFIELD columns") {
+
+    const double NaN = std::numeric_limits<double>::quiet_NaN();
+    std::vector<double> values{1.0, NaN, 3.0};
+
+    for (auto type : {odc::api::INTEGER, odc::api::BITFIELD}) {
+
+        eckit::Buffer buf(4096);
+        eckit::MemoryHandle writeDH(buf);
+        writeDH.openForWrite(0);
+
+        std::vector<odc::api::ColumnInfo> columns{{"col0", type, sizeof(double), {}}};
+        std::vector<odc::api::ConstStridedData> strides{{values.data(), values.size(), sizeof(double), sizeof(double)}};
+        EXPECT_THROWS_AS(odc::api::encode(writeDH, columns, strides), eckit::UserError);
+    }
+}
+
+
 // ------------------------------------------------------------------------------------------------------
 
 int main(int argc, char* argv[]) {

tests/core/test_comparator.cc

@@ -0,0 +1,45 @@
+/*
+ * (C) Copyright 1996-2012 ECMWF.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation nor
+ * does it submit to any jurisdiction.
+ */
+
+#include <limits>
+
+#include "eckit/testing/Test.h"
+
+#include "odc/Comparator.h"
+
+using namespace eckit::testing;
+
+// ------------------------------------------------------------------------------------------------------
+
+CASE("Comparator treats two NaNs as equal") {
+    const double NaN = std::numeric_limits<double>::quiet_NaN();
+
+    EXPECT(odc::Comparator::same(NaN, NaN));
+    EXPECT(!odc::Comparator::same(NaN, 1.0));
+    EXPECT(!odc::Comparator::same(1.0, NaN));
+}
+
+CASE("Comparator treats identical values as equal") {
+    EXPECT(odc::Comparator::same(0.0, 0.0));
+    EXPECT(odc::Comparator::same(1.0, 1.0));
+    EXPECT(odc::Comparator::same(-1.0, -1.0));
+}
+
+CASE("Comparator treats different values as unequal") {
+    EXPECT(!odc::Comparator::same(0.0, 1.0));
+    EXPECT(!odc::Comparator::same(1.0, 2.0));
+    EXPECT(!odc::Comparator::same(1.0, -1.0));
+}
+
+// ------------------------------------------------------------------------------------------------------
+
+int main(int argc, char* argv[]) {
+    return run_tests(argc, argv);
+}