feat: Data history storage (for averaging) (#2206)

Co-authored-by: Tristan Youngs <trisyoungs@googlemail.com>

Author: trisyoungs

src/classes/partialSet.cpp

@@ -499,6 +499,13 @@ void PartialSet::operator*=(const double factor)
     unboundTotal_ *= factor;
 }
 
+PartialSet PartialSet::operator*(const double factor) const
+{
+    auto result = (*this);
+    result *= factor;
+    return result;
+}
+
 /*
  * Searchers
  */

src/classes/partialSet.h

@@ -135,6 +135,7 @@ class PartialSet
     void operator+=(const PartialSet &source);
     void operator-=(const double delta);
     void operator*=(const double factor);
+    PartialSet operator*(const double factor) const;
 
     /*
      * Searchers

src/math/CMakeLists.txt

@@ -39,6 +39,7 @@ add_library(
   histogram2D.h
   histogram3D.cpp
   histogram3D.h
+  history.h
   integerHistogram1D.cpp
   integerHistogram1D.h
   integrator.cpp

src/math/history.h

@@ -0,0 +1,45 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2025 Team Dissolve and contributors
+
+#pragma once
+
+#include "base/serialiser.h"
+#include <memory>
+#include <vector>
+
+// Data History
+template <class T> class History
+{
+    private:
+    // Stored historical data
+    std::vector<std::unique_ptr<T>> history_;
+
+    public:
+    // Update history with supplied data and return current average
+    T average(const T &currentData, int averagingLength)
+    {
+        // Push the current data onto the history stack
+        history_.emplace_back(std::make_unique<T>(currentData));
+
+        // Prune old data to get to the averagingLength
+        while (history_.size() > averagingLength)
+            history_.erase(history_.begin());
+
+        // Perform averaging of the datasets that we have
+        T averaged;
+        auto weight = 1.0 / history_.size();
+        for (auto &data : history_)
+            averaged += *data * weight;
+
+        return averaged;
+    };
+    // Express data as a serialisable value
+    SerialisedValue serialise()
+        requires(std::is_base_of_v<Serialisable<>, T>)
+    {
+        SerialisedValue result;
+        result["size"] = history_.size();
+        Serialisable<>::fromVectorToTable(history_, "data", result);
+        return result;
+    }
+};

src/nodes/gr/gr.cpp

@@ -11,8 +11,6 @@ GRNode::GRNode(Graph *parentGraph) : Node(parentGraph)
     addOption<std::optional<Number>>("Range", "Maximum r to calculate g(r) out to", requestedRange_);
     addOption<std::optional<Number>>("Averaging", "Number of historical partial sets to combine into final partials",
                                      averagingLength_);
-    addOption<Averaging::AveragingScheme>("AveragingScheme", "Weighting scheme to use when averaging partials",
-                                          averagingScheme_);
     addOption<Function1DWrapper>("IntraBroadening", "Type of broadening to apply to intramolecular g(r)", intraBroadening_);
     addOption<std::optional<Number>>("Smoothing", "Specifies the degree of smoothing to apply to calculated g(r)", nSmooths_);
     addOption<bool>("Save", "Whether to save partials and total functions to disk", save_);

src/nodes/gr/gr.h

@@ -8,13 +8,10 @@
 #include "classes/partialSet.h"
 #include "classes/species.h"
 #include "items/list.h"
-#include "math/averaging.h"
-#include "math/data1D.h"
 #include "math/function1D.h"
-#include "nodes/graph.h"
+#include "math/history.h"
 #include "nodes/node.h"
 #include "nodes/number.h"
-#include "nodes/parameter.h"
 #include <vector>
 
 class GRNode : public Node
@@ -47,12 +44,12 @@ class GRNode : public Node
     Configuration *targetConfiguration_{nullptr};
     // Raw simulation g(r)
     std::optional<PartialSet> rawGR_;
+    // Historical raw g(r)
+    History<PartialSet> rawGRHistory_;
     // Unweighted g(r)
     std::optional<PartialSet> unweightedGR_;
     // Number of historical partial sets to combine into final partials
     std::optional<Number> averagingLength_{5};
-    // Weighting scheme to use when averaging partials
-    Averaging::AveragingScheme averagingScheme_{Averaging::LinearAveraging};
     // Bin width (spacing in r) to use
     Number binWidth_{0.001};
     // Perform internal check of calculated partials against a set calculated by a simple unoptimised double-loop

src/nodes/gr/process.cpp

@@ -25,8 +25,7 @@ NodeConstants::ProcessResult GRNode::process()
         message("Partials will be calculated out to {} Angstroms.\n", requestedRange_.value().asDouble());
     message("Bin-width to use is {} Angstroms.\n", binWidth_.asDouble());
     if (averagingLength_)
-        message("Partials will be averaged over {} sets (scheme = {}).\n", averagingLength_.value().asDouble(),
-                Averaging::averagingSchemes().keyword(averagingScheme_));
+        message("Partials will be averaged over {} sets.\n", averagingLength_.value().asDouble());
     else
         message("No averaging of partials will be performed.\n");
     if (intraBroadening_.form() == Functions1D::Form::None)
@@ -80,17 +79,9 @@ NodeConstants::ProcessResult GRNode::process()
     bool alreadyUpToDate;
     calculateRawGR(grRange, alreadyUpToDate);
 
-    // Perform averagingLength_ of unweighted partials if requested, and if we're not already up-to-date
-    /*
+    // Perform averaging of unweighted partials if requested, and if we're not already up-to-date
     if ((averagingLength_.value_or(1) > 1) && (!alreadyUpToDate))
-    {
-        // Store the current fingerprint, since we must ensure we retain it in the averaged T.
-        std::string currentFingerprint{rawGR_.fingerprint()};
-
-        Averaging::average<PartialSet>(dissolve().processingModuleData(), std::format("{}//OriginalGR",
-    targetConfiguration_->niceName()), name(), averagingLength_.value().asDouble(), averagingScheme_);
-    }
-    */
+        (*rawGR_) = rawGRHistory_.average(*rawGR_, averagingLength_.value().asInteger());
 
     /*
     // Perform internal test of original g(r)?