Difference operators in statistics action

src/multio/action/encode/GribEncoder.cc

@@ -57,7 +57,7 @@ const std::map<const std::string, const std::int64_t> ops_to_code{
     {"instant", 0000}, {"average", 1000}, {"accumulate", 2000}, {"maximum", 3000}, {"minimum", 4000}, {"stddev", 5000}};
 
 const std::map<const std::string, const std::int64_t> type_of_statistical_processing{
-    {"average", 0}, {"accumulate", 1}, {"maximum", 2}, {"minimum", 3}, {"stddev", 6}};
+    {"average", 0}, {"accumulate", 1}, {"maximum", 2}, {"minimum", 3}, {"difference", 4}, {"stddev", 6}, {"inverse-difference", 8}};
 
 const std::map<const std::string, const std::string> category_to_levtype{
     {"ocean-grid-coordinate", "oceanSurface"}, {"ocean-2d", "oceanSurface"}, {"ocean-3d", "oceanModelLevel"}};
@@ -751,18 +751,14 @@ void setDateAndStatisticalFields(GribEncoder& g, const message::Metadata& in,
         }
 
         if (operation && (*operation != "instant")) {
-            static const std::map<const std::string, const std::int64_t> TYPE_OF_STATISTICAL_PROCESSING{
-                {"average", 0}, {"accumulate", 1}, {"maximum", 2}, {"minimum", 3}, {"stddev", 6}};
-            if (auto searchStat = TYPE_OF_STATISTICAL_PROCESSING.find(*operation);
-                searchStat != TYPE_OF_STATISTICAL_PROCESSING.end()) {
-                g.setValue("typeOfStatisticalProcessing", searchStat->second);
-            }
-            else {
+            const auto searchStat = type_of_statistical_processing.find(*operation);
+            if (searchStat == std::end(type_of_statistical_processing)) {
                 std::ostringstream oss;
                 oss << "setDateAndStatisticalFields - Cannot map value \"" << *operation
                     << "\"for key \"operation\" (statistical output) to a valid grib2 type of statistical processing.";
                 throw eckit::UserError(oss.str(), Here());
             }
+            g.setValue("typeOfStatisticalProcessing", searchStat->second);
         }
 
 

src/multio/action/statistics/CMakeLists.txt

@@ -25,6 +25,8 @@ list( APPEND _statistics_sources
     operations/OperationWithDeaccumulatedData.h
     operations/Accumulate.h
     operations/Average.h
+    operations/Difference.h
+    operations/InverseDifference.h
     operations/FluxAverage.h
     operations/Instant.h
     operations/Minimum.h

src/multio/action/statistics/Operations.h

@@ -20,6 +20,8 @@
 
 #include "multio/action/statistics/operations/Accumulate.h"
 #include "multio/action/statistics/operations/Average.h"
+#include "multio/action/statistics/operations/Difference.h"
+#include "multio/action/statistics/operations/InverseDifference.h"
 #include "multio/action/statistics/operations/FluxAverage.h"
 #include "multio/action/statistics/operations/Instant.h"
 #include "multio/action/statistics/operations/Maximum.h"
@@ -40,6 +42,12 @@ std::unique_ptr<Operation> make_operation(const std::string& opname, long sz, st
     if (opname == "average") {
         return std::make_unique<Average<Precision>>(opname, sz, win, cfg);
     }
+    if (opname == "difference") {
+        return std::make_unique<Difference<Precision>>(opname, sz, win, cfg);
+    }
+    if (opname == "inverse-difference") {
+        return std::make_unique<InverseDifference<Precision>>(opname, sz, win, cfg);
+    }
     if (opname == "flux-average") {
         return std::make_unique<FluxAverage<Precision>>(opname, sz, win, cfg);
     }
@@ -80,6 +88,14 @@ std::unique_ptr<Operation> load_operation(const std::string& opname, std::shared
         found = true;
         ret = std::make_unique<Average<Precision>>(opname, win, IOmanager, opt);
     }
+    if (opname == "difference") {
+        found = true;
+        ret = std::make_unique<Difference<Precision>>(opname, win, IOmanager, opt);
+    }
+    if (opname == "inverse-difference") {
+        found = true;
+        ret = std::make_unique<InverseDifference<Precision>>(opname, win, IOmanager, opt);
+    }
     if (opname == "flux-average") {
         found = true;
         ret = std::make_unique<FluxAverage<Precision>>(opname, win, IOmanager, opt);

src/multio/action/statistics/operations/Difference.h

@@ -0,0 +1,61 @@
+
+#pragma once
+
+#include "multio/LibMultio.h"
+#include "multio/action/statistics/operations/OperationWithDeaccumulatedData.h"
+
+namespace multio::action {
+
+template <typename T, typename = std::enable_if_t<std::is_floating_point_v<T>>>
+class Difference final : public OperationWithDeaccumulatedData<T> {
+public:
+    using OperationWithDeaccumulatedData<T>::name_;
+    using OperationWithDeaccumulatedData<T>::logHeader_;
+    using OperationWithDeaccumulatedData<T>::initValues_;
+    using OperationWithDeaccumulatedData<T>::values_;
+    using OperationWithDeaccumulatedData<T>::win_;
+    using OperationWithDeaccumulatedData<T>::checkSize;
+    using OperationWithDeaccumulatedData<T>::checkTimeInterval;
+
+
+    Difference(const std::string& name, long sz, const OperationWindow& win, const StatisticsConfiguration& cfg) :
+        OperationWithDeaccumulatedData<T>{name, "difference", sz, true, win, cfg} {}
+
+    Difference(const std::string& name, const OperationWindow& win, std::shared_ptr<StatisticsIO>& IOmanager,
+               const StatisticsOptions& opt) :
+        OperationWithDeaccumulatedData<T>{name, "difference", true, win, IOmanager, opt} {};
+
+    void compute(eckit::Buffer& buf, const StatisticsConfiguration& cfg) override {
+        checkTimeInterval(cfg);
+        LOG_DEBUG_LIB(LibMultio) << logHeader_ << ".compute().count=" << win_.count() << std::endl;
+        auto val = static_cast<T*>(buf.data());
+        cfg.bitmapPresent() ? computeWithMissing(val, cfg) : computeWithoutMissing(val, cfg);
+        return;
+    }
+
+    void updateData(const void* data, long sz, const StatisticsConfiguration& cfg) override {
+        checkSize(sz, cfg);
+        LOG_DEBUG_LIB(LibMultio) << logHeader_ << ".update().count=" << win_.count() << std::endl;
+        auto val = static_cast<const T*>(data);
+        std::copy(val, val + (sz / sizeof(T)), values_.begin());
+        return;
+    }
+
+private:
+    void computeWithoutMissing(T* val, const StatisticsConfiguration& cfg) {
+        std::transform(values_.begin(), values_.end(), initValues_.begin(), val,
+                       [](T v1, T v2) { return static_cast<T>(v1 - v2); });
+        return;
+    }
+
+    void computeWithMissing(T* val, const StatisticsConfiguration& cfg) {
+        double m = cfg.missingValue();
+        std::transform(values_.begin(), values_.end(), initValues_.begin(), val,
+                       [m](T v1, T v2) { return static_cast<T>(m == v1 || m == v2 ? m : v1 - v2); });
+        return;
+    }
+
+    void print(std::ostream& os) const override { os << logHeader_; }
+};
+
+}  // namespace multio::action

src/multio/action/statistics/operations/InverseDifference.h

@@ -0,0 +1,62 @@
+
+#pragma once
+
+#include "multio/LibMultio.h"
+#include "multio/action/statistics/operations/OperationWithDeaccumulatedData.h"
+
+namespace multio::action {
+
+template <typename T, typename = std::enable_if_t<std::is_floating_point_v<T>>>
+class InverseDifference final : public OperationWithDeaccumulatedData<T> {
+public:
+    using OperationWithDeaccumulatedData<T>::name_;
+    using OperationWithDeaccumulatedData<T>::logHeader_;
+    using OperationWithDeaccumulatedData<T>::initValues_;
+    using OperationWithDeaccumulatedData<T>::values_;
+    using OperationWithDeaccumulatedData<T>::win_;
+    using OperationWithDeaccumulatedData<T>::checkSize;
+    using OperationWithDeaccumulatedData<T>::checkTimeInterval;
+
+
+    InverseDifference(const std::string& name, long sz, const OperationWindow& win, const StatisticsConfiguration& cfg) :
+        OperationWithDeaccumulatedData<T>{name, "inverse-difference", sz, true, win, cfg} {}
+
+    InverseDifference(const std::string& name, const OperationWindow& win, std::shared_ptr<StatisticsIO>& IOmanager,
+                      const StatisticsOptions& opt) :
+        OperationWithDeaccumulatedData<T>{name, "inverse-difference", true, win, IOmanager, opt} {};
+
+    void compute(eckit::Buffer& buf, const StatisticsConfiguration& cfg) override {
+        checkTimeInterval(cfg);
+        LOG_DEBUG_LIB(LibMultio) << logHeader_ << ".compute().count=" << win_.count() << std::endl;
+        auto val = static_cast<T*>(buf.data());
+        cfg.bitmapPresent() ? computeWithMissing(val, cfg) : computeWithoutMissing(val, cfg);
+        return;
+    }
+
+    void updateData(const void* data, long sz, const StatisticsConfiguration& cfg) override {
+        checkSize(sz, cfg);
+        LOG_DEBUG_LIB(LibMultio) << logHeader_ << ".update().count=" << win_.count() << std::endl;
+        auto val = static_cast<const T*>(data);
+        std::copy(val, val + (sz / sizeof(T)), values_.begin());
+        return;
+    }
+
+private:
+
+    void computeWithoutMissing(T* val, const StatisticsConfiguration& cfg) {
+        std::transform(values_.begin(), values_.end(), initValues_.begin(), val,
+                       [](T v1, T v2) { return static_cast<T>(v2 - v1); });
+        return;
+    }
+
+    void computeWithMissing(T* val, const StatisticsConfiguration& cfg) {
+        double m = cfg.missingValue();
+        std::transform(values_.begin(), values_.end(), initValues_.begin(), val,
+                       [m](T v1, T v2) { return static_cast<T>(m == v1 || m == v2 ? m : v2 - v1); });
+        return;
+    }
+
+    void print(std::ostream& os) const override { os << logHeader_; }
+};
+
+}  // namespace multio::action

tests/multio/action/statistics/CMakeLists.txt

@@ -126,4 +126,20 @@ add_subdirectory(restart)
 # add moving mask tests
 add_subdirectory(moving-mask)
 
+ecbuild_add_test(
+    TARGET ${PREFIX}_difference
+    SOURCES test_difference.cc
+    NO_AS_NEEDED
+    LIBS multio-action-debug-sink multio-action-statistics
+    ENVIRONMENT "MULTIO_SERVER_CONFIG_PATH=${CMAKE_CURRENT_SOURCE_DIR}/config"
+)
+
+ecbuild_add_test(
+    TARGET ${PREFIX}_inverse_difference
+    SOURCES test_inverse_difference.cc
+    NO_AS_NEEDED
+    LIBS multio-action-debug-sink multio-action-statistics
+    ENVIRONMENT "MULTIO_SERVER_CONFIG_PATH=${CMAKE_CURRENT_SOURCE_DIR}/config"
+)
+
 endif()

tests/multio/action/statistics/testStatisticsUtils.h

@@ -0,0 +1,95 @@
+
+# pragma once
+
+#include <random>
+
+#include "eckit/config/YAMLConfiguration.h"
+#include "eckit/testing/Test.h"
+
+#include "multio/action/Plan.h"
+#include "multio/config/MultioConfiguration.h"
+#include "multio/config/PathConfiguration.h"
+
+namespace multio::test {
+
+
+namespace {
+    config::ConfigAndPaths configFromActions(const std::string& actions) {
+        config::ConfigAndPaths configAndPaths;
+        configAndPaths.paths = config::defaultConfigPaths();
+        configAndPaths.parsedConfig = eckit::LocalConfiguration{eckit::YAMLConfiguration(actions)};
+        return configAndPaths;
+    }
+}
+
+
+class MultioTestEnvironment {
+public:
+    MultioTestEnvironment() = delete;
+
+    MultioTestEnvironment(config::MultioConfiguration multioConfig) :
+        multioConfig_{std::move(multioConfig)},
+        plan_{config::ComponentConfiguration(multioConfig_.parsedConfig(), multioConfig_)} {}
+
+    MultioTestEnvironment(const std::string& actions) : MultioTestEnvironment(configFromActions(actions)) {}
+
+    config::MultioConfiguration& multioConfig() {
+        return multioConfig_;
+    }
+
+    multio::action::Plan& plan() {
+        return plan_;
+    }
+
+    std::queue<message::Message>& debugSink() {
+        return multioConfig_.debugSink();
+    }
+
+private:
+    config::MultioConfiguration multioConfig_;
+    multio::action::Plan plan_;
+};
+
+
+template<typename T>
+std::vector<T> randomVector(size_t size, T min, T max, std::uint32_t seed = 42) {
+    std::mt19937 gen(seed);
+
+    std::vector<T> v(size);
+
+    if constexpr(std::is_integral<T>::value) {
+        std::uniform_int_distribution<T> dis(min, max);
+        std::transform(v.begin(), v.end(), v.begin(), [&dis, &gen]() { return dis(gen); });
+    }
+    else if constexpr(std::is_floating_point<T>::value) {
+        std::uniform_real_distribution<T> dis(min, max);
+        std::transform(v.begin(), v.end(), v.begin(), [&dis, &gen](T val) { return dis(gen); });
+    }
+
+    return v;
+}
+
+
+multio::message::Message createStatisticsMessage(long step, std::vector<double> payload) {
+    eckit::Buffer payloadBuf{payload.data(), payload.size() * sizeof(double)};
+    multio::message::Metadata metadata{{
+        { "paramId", 0 },
+        { "level", 0 },
+        { "levtype", "none" },
+        { "gridType", "none" },
+        { "precision", "double" },
+        { "startDate", 0 },
+        { "startTime", 0 },
+        { "step", step },
+        { "bitmapPresent", false },
+        { "missingValue", -999.0 }
+    }};
+
+    return multio::message::Message{
+        {multio::message::Message::Tag::Field, {}, {}, std::move(metadata)},
+        std::move(payloadBuf)
+    };
+}
+
+
+}  // namespace multio::test