Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions src/classes/partialSet.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
*/
Expand Down
1 change: 1 addition & 0 deletions src/classes/partialSet.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions src/math/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@ add_library(
histogram2D.h
histogram3D.cpp
histogram3D.h
history.h
integerHistogram1D.cpp
integerHistogram1D.h
integrator.cpp
Expand Down
45 changes: 45 additions & 0 deletions src/math/history.h
Original file line number Diff line number Diff line change
@@ -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;
Comment on lines +24 to +32

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
// 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;
// How many items to average
auto length = history.size() < averagingLength ? history.size() : averagingLength;
// Perform averaging of the datasets that we have
T averaged;
auto weight = 1.0 / length;
for (auto &data : std::span(history_.rbegin(), ristory_.rbegin()+length))
averaged += *data * weight;

This is an alternate implementation that doesn't lose history. This would allow us to look at multiple averaging lengths simultaneously (the current implementation essentially locks you into the shortest length). The disadvantage is that it can grow without bound, though that could be fixed by replacing the vector with a ring buffer.

This isn't a necessary change - just a suggestion.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Interesting thought. While I like the idea, my concern is that this could add a lot of bloat as many of the objects being stored here are pretty chunky (I'm looking at you, PartialSet). As such, I would prefer to stick to the ultra-simple, shortest-length version for the time being.


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;
}
};
2 changes: 0 additions & 2 deletions src/nodes/gr/gr.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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_);
Expand Down
9 changes: 3 additions & 6 deletions src/nodes/gr/gr.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
15 changes: 3 additions & 12 deletions src/nodes/gr/process.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)?
Expand Down
Loading