New operator MDR, for refactoring floating point arrays using MGARD's…

… new MDR extension. (#3826)

- Only works with BP5 for now since BP5 deserializer had to be modified for reading with accuracy
 - Use mgard cmake flag MGARD_ENABLE_MDR to set ADIOS2_HAVE_MGARD_MDR flag and build refactor operator conditionally
 - Added Accuracy struct to adios types (error, norm, and rel|abs flag)
 - Set/GetAccuracy() functions to Variable, VariableBase and Operator classes.
 - BP5 deserializer now prepares the operator with accuracy before calling Decompress().
 - bpls new option --error "double,double|inf,rel|abs" added.
 - Add tests for MDR operator and Accuracy default behavior (for all operators)
 - Add rounding error constants adios2::ops::mdr::DOUBLE_ROUNDING_ERROR_LIMIT and FLOAT_ROUNDING_ERROR_LIMIT,
      beacuse GetAccuracy() is limited by rounding errors in MDR

CMakeLists.txt

@@ -243,7 +243,7 @@ endif()
 
 set(ADIOS2_CONFIG_OPTS
     DataMan DataSpaces HDF5 HDF5_VOL MHS SST Fortran MPI Python Blosc2 BZip2
-    LIBPRESSIO MGARD PNG SZ ZFP DAOS IME O_DIRECT Sodium Catalyst SysVShMem UCX
+    LIBPRESSIO MGARD MGARD_MDR PNG SZ ZFP DAOS IME O_DIRECT Sodium Catalyst SysVShMem UCX
     ZeroMQ Profiling Endian_Reverse Derived_Variable AWSSDK GPU_Support CUDA Kokkos
     Kokkos_CUDA Kokkos_HIP Kokkos_SYCL
 )

bindings/CXX11/adios2/cxx11/Variable.cpp

@@ -78,6 +78,13 @@ namespace adios2
     }                                                                                              \
                                                                                                    \
     template <>                                                                                    \
+    void Variable<T>::SetAccuracy(const adios2::Accuracy &a)                                       \
+    {                                                                                              \
+        helper::CheckForNullptr(m_Variable, "in call to Variable<T>::SetAccuracy");                \
+        m_Variable->SetAccuracy(a);                                                                \
+    }                                                                                              \
+                                                                                                   \
+    template <>                                                                                    \
     size_t Variable<T>::SelectionSize() const                                                      \
     {                                                                                              \
         helper::CheckForNullptr(m_Variable, "in call to Variable<T>::SelectionSize");              \
@@ -219,6 +226,13 @@ namespace adios2
     }                                                                                              \
                                                                                                    \
     template <>                                                                                    \
+    adios2::Accuracy Variable<T>::GetAccuracy()                                                    \
+    {                                                                                              \
+        helper::CheckForNullptr(m_Variable, "in call to Variable<T>::GetAccuracy");                \
+        return m_Variable->GetAccuracy();                                                          \
+    }                                                                                              \
+                                                                                                   \
+    template <>                                                                                    \
     std::vector<std::vector<typename Variable<T>::Info>> Variable<T>::AllStepsBlocksInfo()         \
     {                                                                                              \
         return DoAllStepsBlocksInfo();                                                             \

bindings/CXX11/adios2/cxx11/Variable.h

@@ -209,6 +209,12 @@ class Variable
      */
     void SetStepSelection(const adios2::Box<size_t> &stepSelection);
 
+    /**
+     * Sets the requested accuracy for the next read operation.
+     * The actual accuracy after the read is provided in GetAccuracy()
+     */
+    void SetAccuracy(const adios2::Accuracy &a);
+
     /**
      * Returns the number of elements required for pre-allocation based on
      * current count and stepsCount
@@ -335,6 +341,13 @@ class Variable
      */
     T Max(const size_t step = adios2::DefaultSizeT) const;
 
+    /**
+     * Get the provided accuracy for the last read operation.
+     * Most operations provide data as it was written, meaning that
+     * error is reported as 0.0
+     */
+    adios2::Accuracy GetAccuracy();
+
     /** Contains block information for a particular Variable<T> */
     struct Info
     {

cmake/DetectOptions.cmake

@@ -175,13 +175,17 @@ if(TARGET LibPressio::libpressio)
 endif()
 
 # MGARD
+set(ADIOS2_HAVE_MGARD_MDR FALSE)
 if(ADIOS2_USE_MGARD STREQUAL AUTO)
   find_package(mgard CONFIG)
 elseif(ADIOS2_USE_MGARD)
   find_package(mgard REQUIRED CONFIG)
 endif()
 if(mgard_FOUND)
   set(ADIOS2_HAVE_MGARD TRUE)
+  if(MGARD_ENABLE_MDR)
+     set(ADIOS2_HAVE_MGARD_MDR TRUE)
+ endif()
 endif()
 
 # PNG

examples/basics/globalArrayND/globalArrayNDWrite.cpp

@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
     MPI_Comm_rank(MPI_COMM_WORLD, &rank);
     MPI_Comm_size(MPI_COMM_WORLD, &nproc);
 #endif
-    const int NSTEPS = 5;
+    const int NSTEPS = 1;
 
 #if ADIOS2_USE_MPI
     adios2::ADIOS adios(MPI_COMM_WORLD);
@@ -58,7 +58,7 @@ int main(int argc, char *argv[])
 #endif
 
     // Application variables for output
-    const unsigned int Nx = 10;
+    const unsigned int Nx = 100000;
     // Global 2D array, size of nproc x Nx, with 1D decomposition
     // Each process writes one "row" of the 2D matrix.
     std::vector<double> row(Nx);
@@ -84,6 +84,9 @@ int main(int argc, char *argv[])
 
         io.DefineAttribute<int>("nsteps", NSTEPS);
 
+        adios2::Operator op = adios.DefineOperator("mdr", "mdr");
+        varGlobalArray.AddOperation(op, {{"accuracy", std::to_string(0.1)}});
+
         // Open file. "w" means we overwrite any existing file on disk,
         // but Advance() will append steps to the same file.
         adios2::Engine writer = io.Open("globalArray.bp", adios2::Mode::Write);

source/adios2/CMakeLists.txt

@@ -323,6 +323,9 @@ endif()
 if(ADIOS2_HAVE_MGARD)
   target_sources(adios2_core PRIVATE operator/compress/CompressMGARDPlus.cpp)
   target_sources(adios2_core PRIVATE operator/compress/CompressMGARD.cpp)
+  if(ADIOS2_HAVE_MGARD_MDR)
+      target_sources(adios2_core PRIVATE operator/refactor/RefactorMDR.cpp)
+  endif()
   target_link_libraries(adios2_core PRIVATE mgard::mgard)
 endif()
 

source/adios2/common/ADIOSTypes.h

@@ -310,6 +310,22 @@ using Box = std::pair<T, T>;
 template <typename T, typename Enable = void>
 struct TypeInfo;
 
+/** Data accuracy **/
+
+/* Error. Accuracy can be requested for reading data.
+   norm: 0.0 = L2, inf() = Linf
+   relative: true = relative error, false = absolute error
+ */
+struct Accuracy
+{
+    double error;
+    double norm;
+    bool relative;
+};
+
+constexpr double L2_norm = 0.0;
+constexpr double Linf_norm = std::numeric_limits<double>::infinity();
+
 /**
  *  Return the actual size in bytes of elements of the given type.  Returns -1
  * for strings.
@@ -405,6 +421,19 @@ constexpr char s[] = "s";
 }
 #endif
 
+#ifdef ADIOS2_HAVE_MGARD_MDR
+
+constexpr char MDR[] = "mdr";
+
+namespace mdr
+{
+constexpr double DOUBLE_ROUNDING_ERROR_LIMIT = 5.0e-16;
+constexpr double FLOAT_ROUNDING_ERROR_LIMIT = 3.0e-7;
+namespace key
+{}
+}
+#endif
+
 #ifdef ADIOS2_HAVE_LIBPRESSIO
 constexpr char LossyLIBPRESSIO[] = "libpressio";
 namespace libpressio

source/adios2/core/Info.cpp

@@ -72,6 +72,9 @@ static const char *const operators[] = {
     "MGARD",
     "MGARDPlus",
 #endif
+#ifdef ADIOS2_HAVE_MGARD_MDR
+    "MDR",
+#endif
 #ifdef ADIOS2_HAVE_SZ
     "SZ",
 #endif

source/adios2/core/Operator.cpp

@@ -11,6 +11,8 @@
 #include "Operator.h"
 #include "adios2/helper/adiosFunctions.h"
 
+#include <iostream>
+
 namespace adios2
 {
 namespace core
@@ -30,6 +32,9 @@ void Operator::SetParameter(const std::string key, const std::string value) noex
 
 Params &Operator::GetParameters() noexcept { return m_Parameters; }
 
+void Operator::SetAccuracy(const adios2::Accuracy &a) noexcept { m_AccuracyRequested = a; }
+adios2::Accuracy Operator::GetAccuracy() const noexcept { return m_AccuracyProvided; }
+
 #define declare_type(T)                                                                            \
                                                                                                    \
     void Operator::RunCallback1(const T *arg0, const std::string &arg1, const std::string &arg2,   \
@@ -96,6 +101,12 @@ Dims Operator::ConvertDims(const Dims &dimensions, const DataType type, const si
 
 size_t Operator::GetHeaderSize() const { return 0; }
 
+size_t Operator::GetEstimatedSize(const size_t ElemCount, const size_t ElemSize, const size_t ndims,
+                                  const size_t *dims) const
+{
+    return ElemCount * ElemSize + 128;
+};
+
 // PRIVATE
 void Operator::CheckCallbackType(const std::string type) const
 {

source/adios2/core/Operator.h

@@ -37,6 +37,7 @@ class Operator
         COMPRESS_SZ = 6,
         COMPRESS_ZFP = 7,
         COMPRESS_MGARDPLUS = 8,
+        REFACTOR_MDR = 41,
         CALLBACK_SIGNATURE1 = 51,
         CALLBACK_SIGNATURE2 = 52,
         PLUGIN_INTERFACE = 53,
@@ -56,6 +57,9 @@ class Operator
 
     Params &GetParameters() noexcept;
 
+    void SetAccuracy(const adios2::Accuracy &a) noexcept;
+    adios2::Accuracy GetAccuracy() const noexcept;
+
 #define declare_type(T)                                                                            \
     virtual void RunCallback1(const T *, const std::string &, const std::string &,                 \
                               const std::string &, const size_t, const Dims &, const Dims &,       \
@@ -68,6 +72,10 @@ class Operator
 
     virtual size_t GetHeaderSize() const;
 
+    /** Give an upper bound estimate how big the transformed data could be */
+    virtual size_t GetEstimatedSize(const size_t ElemCount, const size_t ElemSize,
+                                    const size_t ndims, const size_t *dims) const;
+
     /**
      * @param dataIn
      * @param blockStart
@@ -94,6 +102,11 @@ class Operator
     /** Parameters associated with a particular Operator */
     Params m_Parameters;
 
+    /** user requested accuracy */
+    Accuracy m_AccuracyRequested = {0.0, 0.0, false};
+    /** provided accuracy */
+    Accuracy m_AccuracyProvided = {0.0, 0.0, false};
+
     /**
      * Used by lossy compressors with a limitation on complex data types or
      * dimentions Returns a adios2::Dims object that meets the requirement of a

source/adios2/core/VariableBase.cpp

@@ -204,6 +204,14 @@ void VariableBase::SetMemorySelection(const Box<Dims> &memorySelection)
     m_MemoryCount = memorySelection.second;
 }
 
+void VariableBase::SetAccuracy(const adios2::Accuracy &a) noexcept
+{
+    m_AccuracyRequested = a;
+    m_AccuracyProvided = {0.0, a.norm, a.relative};
+}
+adios2::Accuracy VariableBase::GetAccuracy() const noexcept { return m_AccuracyProvided; }
+adios2::Accuracy VariableBase::GetAccuracyRequested() const noexcept { return m_AccuracyRequested; }
+
 size_t VariableBase::GetAvailableStepsStart() const { return m_AvailableStepsStart; }
 
 size_t VariableBase::GetAvailableStepsCount() const { return m_AvailableStepsCount; }

source/adios2/core/VariableBase.h

@@ -98,6 +98,11 @@ class VariableBase
 
     Engine *m_Engine = nullptr;
 
+    /** user requested accuracy */
+    Accuracy m_AccuracyRequested = {0.0, 0.0, false};
+    /** provided accuracy */
+    Accuracy m_AccuracyProvided = {0.0, 0.0, false};
+
     /** Index to Step and blocks' (inside a step) characteristics position in a
      * serial metadata buffer
      * <pre>
@@ -166,6 +171,22 @@ class VariableBase
      */
     void SetMemorySelection(const Box<Dims> &memorySelection);
 
+    /**
+     * Sets the requested accuracy for the next read operation.
+     * The actual accuracy after the read is provided in GetAccuracy()
+     */
+    void SetAccuracy(const adios2::Accuracy &a) noexcept;
+
+    /**
+     * Get the provided accuracy for the last read operation.
+     * Most operations provide data as it was written, meaning that
+     * error is reported as 0.0
+     */
+    adios2::Accuracy GetAccuracy() const noexcept;
+
+    /** Return the requested accuracy set by user with SetAccuracy */
+    adios2::Accuracy GetAccuracyRequested() const noexcept;
+
     size_t GetAvailableStepsStart() const;
 
     size_t GetAvailableStepsCount() const;

source/adios2/operator/OperatorFactory.cpp

@@ -31,6 +31,10 @@
 #include "adios2/operator/compress/CompressMGARDPlus.h"
 #endif
 
+#ifdef ADIOS2_HAVE_MGARD_MDR
+#include "adios2/operator/refactor/RefactorMDR.h"
+#endif
+
 #ifdef ADIOS2_HAVE_PNG
 #include "adios2/operator/compress/CompressPNG.h"
 #endif
@@ -74,6 +78,8 @@ std::string OperatorTypeToString(const Operator::OperatorType type)
         return "sz";
     case Operator::COMPRESS_ZFP:
         return "zfp";
+    case Operator::REFACTOR_MDR:
+        return "mdr";
     case Operator::PLUGIN_INTERFACE:
         return "plugin";
     default:
@@ -139,6 +145,12 @@ std::shared_ptr<Operator> MakeOperator(const std::string &type, const Params &pa
     {
 #ifdef ADIOS2_HAVE_ZFP
         ret = std::make_shared<compress::CompressZFP>(parameters);
+#endif
+    }
+    else if (typeLowerCase == "mdr")
+    {
+#ifdef ADIOS2_HAVE_MGARD_MDR
+        ret = std::make_shared<refactor::RefactorMDR>(parameters);
 #endif
     }
     else if (typeLowerCase == "plugin")

source/adios2/operator/OperatorFactory.h

@@ -17,6 +17,8 @@ namespace adios2
 namespace core
 {
 
+std::string OperatorTypeToString(const Operator::OperatorType type);
+
 std::shared_ptr<Operator> MakeOperator(const std::string &type, const Params &parameters);
 
 size_t Decompress(const char *bufferIn, const size_t sizeIn, char *dataOut, MemorySpace memSpace,