Skip to content

Commit

Permalink
Analysis Refactoring - Part 2 (#1795)
Browse files Browse the repository at this point in the history
Co-authored-by: Tristan Youngs <tristan.youngs@stfc.ac.uk>
  • Loading branch information
2 people authored and rhinoella committed Apr 2, 2024
1 parent 59d7a56 commit 6692aab
Show file tree
Hide file tree
Showing 34 changed files with 1,082 additions and 643 deletions.
3 changes: 3 additions & 0 deletions src/analyser/CMakeLists.txt
@@ -1,7 +1,10 @@
add_library(
analyser
dataExporter.h
dataNormaliser1D.cpp
dataNormaliser1D.h
dataNormaliser2D.cpp
dataNormaliser2D.h
dataNormaliser3D.cpp
dataNormaliser3D.h
dataNormaliserBase.h
Expand Down
31 changes: 31 additions & 0 deletions src/analyser/dataExporter.h
@@ -0,0 +1,31 @@
// SPDX-License-Identifier: GPL-3.0-or-later
// Copyright (c) 2024 Team Dissolve and contributors

#pragma once

#include "base/processPool.h"

template <typename DataND, typename DataNDExportFileFormat> class DataExporter
{
public:
// Try to export the specified data, if a valid filename has been provided
static bool exportData(const DataND &targetData, DataNDExportFileFormat &fileAndFormat, const ProcessPool &procPool)
{
if (fileAndFormat.hasFilename())
{
if (procPool.isMaster())
{
if (fileAndFormat.exportData(targetData))
procPool.decideTrue();
else
{
procPool.decideFalse();
return false;
}
}
else if (!procPool.decision())
return false;
}
return true;
}
};
8 changes: 4 additions & 4 deletions src/analyser/dataNormaliser1D.cpp
Expand Up @@ -2,8 +2,6 @@
// Copyright (c) 2024 Team Dissolve and contributors

#include "analyser/dataNormaliser1D.h"
#include "expression/expression.h"
#include "expression/variable.h"
#include "math/data1D.h"
#include "math/integrator.h"

Expand All @@ -14,8 +12,10 @@ void DataNormaliser1D::normalise(NormalisationFunction1D normalisationFunction)
const auto &xs = targetData_.xAxis();
auto &values = targetData_.values();

const auto xDelta = xs.size() > 1 ? xs[1] - xs[0] : 1.0;

for (auto i = 0; i < xs.size(); ++i)
values.at(i) = normalisationFunction(xs[i], values.at(i));
values.at(i) = normalisationFunction(xs[i], xDelta, values.at(i));
}

void DataNormaliser1D::normaliseByGrid() { Messenger::warn("Grid normalisation not implemented for 1D data."); }
Expand Down Expand Up @@ -57,4 +57,4 @@ void DataNormaliser1D::normaliseTo(double value, bool absolute)
auto sum = absolute ? Integrator::absSum(targetData_) : Integrator::sum(targetData_);
targetData_ /= sum;
targetData_ *= value;
}
}
3 changes: 1 addition & 2 deletions src/analyser/dataNormaliser1D.h
Expand Up @@ -5,9 +5,8 @@

#include "analyser/dataNormaliserBase.h"
#include "math/data1D.h"
#include <string_view>

using NormalisationFunction1D = std::function<double(const double &, const double &)>;
using NormalisationFunction1D = std::function<double(const double &, const double &, const double &)>;
class DataNormaliser1D : public DataNormaliserBase<Data1D, NormalisationFunction1D>
{
public:
Expand Down
72 changes: 72 additions & 0 deletions src/analyser/dataNormaliser2D.cpp
@@ -0,0 +1,72 @@
// SPDX-License-Identifier: GPL-3.0-or-later
// Copyright (c) 2024 Team Dissolve and contributors

#include "analyser/dataNormaliser2D.h"
#include "math/data2D.h"
#include "math/integrator.h"

DataNormaliser2D::DataNormaliser2D(Data2D &targetData) : DataNormaliserBase<Data2D, NormalisationFunction2D>(targetData) {}

void DataNormaliser2D::normalise(NormalisationFunction2D normalisationFunction)
{
const auto &xs = targetData_.xAxis();
const auto &ys = targetData_.yAxis();
auto &values = targetData_.values();

const auto xDelta = xs.size() > 1 ? xs[1] - xs[0] : 1.0;
const auto yDelta = ys.size() > 1 ? ys[1] - ys[0] : 1.0;

for (auto i = 0; i < xs.size(); ++i)
{
for (auto j = 0; j < ys.size(); ++j)
{
values[{i, j}] = normalisationFunction(xs[i], xDelta, ys[j], yDelta, values[{i, j}]);
}
}
}

void DataNormaliser2D::normaliseByGrid() { Messenger::warn("Grid normalisation not implemented for 2D data."); }

void DataNormaliser2D::normaliseBySphericalShell()
{
// We expect x values to be centre-bin values, and regularly spaced
const auto &xAxis = targetData_.xAxis();
const auto &yAxis = targetData_.yAxis();
auto &values = targetData_.values();

if (xAxis.size() < 2)
return;

// Derive first left-bin boundary from the delta between points 0 and 1
auto leftBin = xAxis[0] - (xAxis[1] - xAxis[0]) * 0.5;
auto r1Cubed = pow(leftBin, 3);
for (auto n = 0; n < xAxis.size(); ++n)
{
auto r2Cubed = 0.0, rightBin = 0.0;
for (auto m = 0; m < yAxis.size(); ++m)
{
// Get new right-bin from existing left bin boundary and current bin centre
auto rightBin = leftBin + 2 * (xAxis[n] - leftBin);
auto r2Cubed = pow(rightBin, 3);

// Calculate divisor for normalisation
auto divisor = (4.0 / 3.0) * PI * (r2Cubed - r1Cubed);

// Peform normalisation step
values[{n, m}] /= divisor;
if (targetData_.valuesHaveErrors())
targetData_.error(n, m) /= divisor;
}

// Overwrite old values
r1Cubed = r2Cubed;
leftBin = rightBin;
}
}

void DataNormaliser2D::normaliseTo(double value, bool absolute)
{
auto sum = absolute ? Integrator::absSum(targetData_) : Integrator::sum(targetData_);
targetData_ /= sum;
targetData_ *= value;
}
24 changes: 24 additions & 0 deletions src/analyser/dataNormaliser2D.h
@@ -0,0 +1,24 @@
// SPDX-License-Identifier: GPL-3.0-or-later
// Copyright (c) 2024 Team Dissolve and contributors

#pragma once

#include "analyser/dataNormaliserBase.h"
#include "math/data2D.h"

using NormalisationFunction2D =
std::function<double(const double &, const double &, const double &, const double &, const double &)>;
class DataNormaliser2D : public DataNormaliserBase<Data2D, NormalisationFunction2D>
{
public:
DataNormaliser2D(Data2D &targetData);

/*
* Normalisation Functions
*/
public:
void normalise(NormalisationFunction2D normalisationFunction) override;
void normaliseByGrid() override;
void normaliseBySphericalShell() override;
void normaliseTo(double value = 1.0, bool absolute = true) override;
};
4 changes: 2 additions & 2 deletions src/analyser/dataNormaliser3D.h
Expand Up @@ -5,9 +5,9 @@

#include "analyser/dataNormaliserBase.h"
#include "math/data3D.h"
#include <string_view>

using NormalisationFunction3D = std::function<double(const double &, const double &, const double &, const double &)>;
using NormalisationFunction3D = std::function<double(const double &, const double &, const double &, const double &,
const double &, const double &, const double &)>;
class DataNormaliser3D : public DataNormaliserBase<Data3D, NormalisationFunction3D>
{
public:
Expand Down
139 changes: 9 additions & 130 deletions src/modules/angle/angle.cpp
Expand Up @@ -8,140 +8,19 @@
#include "keywords/speciesSiteVector.h"
#include "keywords/vec3Double.h"
#include "module/context.h"
#include "procedure/nodes/calculateAngle.h"
#include "procedure/nodes/calculateDistance.h"
#include "procedure/nodes/collect1D.h"
#include "procedure/nodes/collect2D.h"
#include "procedure/nodes/collect3D.h"
#include "procedure/nodes/operateExpression.h"
#include "procedure/nodes/operateNormalise.h"
#include "procedure/nodes/operateNumberDensityNormalise.h"
#include "procedure/nodes/operateSitePopulationNormalise.h"
#include "procedure/nodes/operateSphericalShellNormalise.h"
#include "procedure/nodes/process1D.h"
#include "procedure/nodes/process2D.h"
#include "procedure/nodes/process3D.h"
#include "procedure/nodes/select.h"

AngleModule::AngleModule() : Module(ModuleTypes::Angle), analyser_(ProcedureNode::AnalysisContext)
AngleModule::AngleModule() : Module(ModuleTypes::Angle)
{
// Select: Site 'A'
selectA_ = analyser_.createRootNode<SelectProcedureNode>("A");
auto &forEachA = selectA_->branch()->get();

// -- Select: Site 'B'
selectB_ = forEachA.create<SelectProcedureNode>("B");
auto &forEachB = selectB_->branch()->get();

// -- -- Calculate: 'rAB'
auto calcAB = forEachB.create<CalculateDistanceProcedureNode>("rAB", selectA_, selectB_);

// -- -- Collect1D: 'RDF(AB)'
collectAB_ =
forEachB.create<Collect1DProcedureNode>({}, calcAB, ProcedureNode::AnalysisContext, rangeAB_.x, rangeAB_.y, rangeAB_.z);

// -- -- Select: Site 'C'
selectC_ = forEachB.create<SelectProcedureNode>("C");
auto &forEachC = selectC_->branch()->get();

// -- -- -- Calculate: 'rBC'
auto calcBC = forEachC.create<CalculateDistanceProcedureNode>({}, selectB_, selectC_);

// -- -- -- Calculate: 'aABC'
calculateAngle_ = forEachC.create<CalculateAngleProcedureNode>({}, selectA_, selectB_, selectC_);
calculateAngle_->keywords().set("Symmetric", symmetric_);

// -- -- -- Collect1D: 'RDF(BC)'
collectBC_ =
forEachC.create<Collect1DProcedureNode>({}, calcBC, ProcedureNode::AnalysisContext, rangeBC_.x, rangeBC_.y, rangeBC_.z);

// -- -- -- Collect1D: 'ANGLE(ABC)'
collectABC_ = forEachC.create<Collect1DProcedureNode>({}, calculateAngle_, ProcedureNode::AnalysisContext, angleRange_.x,
angleRange_.y, angleRange_.z);

// -- -- -- Collect2D: 'DAngle (A-B)-C'
collectDAngleAB_ =
forEachC.create<Collect2DProcedureNode>({}, calcAB, calculateAngle_, ProcedureNode::AnalysisContext, rangeAB_.x,
rangeAB_.y, rangeAB_.z, angleRange_.x, angleRange_.y, angleRange_.z);

// -- -- -- Collect2D: 'DAngle A-(B-C)'
collectDAngleBC_ =
forEachC.create<Collect2DProcedureNode>({}, calcBC, calculateAngle_, ProcedureNode::AnalysisContext, rangeBC_.x,
rangeBC_.y, rangeBC_.z, angleRange_.x, angleRange_.y, angleRange_.z);

// -- -- -- Collect3D: 'rAB vs rBC vs aABC'
collectDDA_ = forEachC.create<Collect3DProcedureNode>({}, calcAB, calcBC, calculateAngle_, ProcedureNode::AnalysisContext,
rangeAB_, rangeBC_, angleRange_);

// Process1D: 'RDF(AB)'
processAB_ = analyser_.createRootNode<Process1DProcedureNode>("RDF(AB)", collectAB_);
processAB_->keywords().set("LabelValue", std::string("g\\sub{AB}(r)"));
processAB_->keywords().set("LabelX", std::string("r, \\symbol{Angstrom}"));
auto &rdfABNormalisation = processAB_->branch()->get();
rdfABNormalisation.create<OperateSitePopulationNormaliseProcedureNode>({},
ConstNodeVector<SelectProcedureNode>({selectA_}));
rdfABNormalisation.create<OperateNumberDensityNormaliseProcedureNode>({}, ConstNodeVector<SelectProcedureNode>({selectB_}));
rdfABNormalisation.create<OperateSphericalShellNormaliseProcedureNode>({});

// Process1D: 'RDF(BC)'
processBC_ = analyser_.createRootNode<Process1DProcedureNode>("RDF(BC)", collectBC_);
processBC_->keywords().set("LabelValue", std::string("g\\sub{BC}(r)"));
processBC_->keywords().set("LabelX", std::string("r, \\symbol{Angstrom}"));
auto &rdfBCNormalisation = processBC_->branch()->get();
rdfBCNormalisation.create<OperateSitePopulationNormaliseProcedureNode>(
{}, ConstNodeVector<SelectProcedureNode>({selectB_, selectA_}));
rdfBCNormalisation.create<OperateNumberDensityNormaliseProcedureNode>({}, ConstNodeVector<SelectProcedureNode>({selectC_}));
rdfBCNormalisation.create<OperateSphericalShellNormaliseProcedureNode>({});

// Process1D: 'ANGLE(ABC)'
processAngle_ = analyser_.createRootNode<Process1DProcedureNode>("Angle(ABC)", collectABC_);
processAngle_->keywords().set("LabelValue", std::string("Normalised Frequency"));
processAngle_->keywords().set("LabelX", std::string("\\symbol{theta}, \\symbol{degrees}"));
auto &angleNormalisation = processAngle_->branch()->get();
angleNormalisation.create<OperateExpressionProcedureNode>({}, "value/sin(toRad(x))");
angleNormalisation.create<OperateNormaliseProcedureNode>({}, 1.0);

// Process2D: 'DAngle (A-B)-C'
processDAngleAB_ = analyser_.createRootNode<Process2DProcedureNode>("DAngle((A-B)-C)", collectDAngleAB_);
processDAngleAB_->keywords().set("LabelValue", std::string("g\\sub{AB}(r)"));
processDAngleAB_->keywords().set("LabelX", std::string("r, \\symbol{Angstrom}"));
processDAngleAB_->keywords().set("LabelY", std::string("\\symbol{theta}, \\symbol{degrees}"));
auto &dAngleABNormalisation = processDAngleAB_->branch()->get();
dAngleABNormalisationExpression_ = dAngleABNormalisation.create<OperateExpressionProcedureNode>(
{}, fmt::format("{} * value/sin(toRad(y))/sin(toRad(yDelta))", symmetric_ ? 1.0 : 2.0));
dAngleABNormalisation.create<OperateSitePopulationNormaliseProcedureNode>(
{}, ConstNodeVector<SelectProcedureNode>({selectA_, selectC_}));
dAngleABNormalisation.create<OperateNumberDensityNormaliseProcedureNode>({},
ConstNodeVector<SelectProcedureNode>({selectB_}));
dAngleABNormalisation.create<OperateSphericalShellNormaliseProcedureNode>({});

// Process2D: 'DAngle A-(B-C)'
processDAngleBC_ = analyser_.createRootNode<Process2DProcedureNode>("DAngle(A-(B-C))", collectDAngleBC_);
processDAngleBC_->keywords().set("LabelValue", std::string("g\\sub{BC}(r)"));
processDAngleBC_->keywords().set("LabelX", std::string("r, \\symbol{Angstrom}"));
processDAngleBC_->keywords().set("LabelY", std::string("\\symbol{theta}, \\symbol{degrees}"));
auto &dAngleBCNormalisation = processDAngleBC_->branch()->get();
dAngleBCNormalisationExpression_ = dAngleBCNormalisation.create<OperateExpressionProcedureNode>(
{}, fmt::format("{} * value/sin(toRad(y))/sin(toRad(yDelta))", symmetric_ ? 1.0 : 2.0));
dAngleBCNormalisation.create<OperateSitePopulationNormaliseProcedureNode>(
{}, ConstNodeVector<SelectProcedureNode>({selectB_, selectA_}));
dAngleBCNormalisation.create<OperateNumberDensityNormaliseProcedureNode>({},
ConstNodeVector<SelectProcedureNode>({selectC_}));
dAngleBCNormalisation.create<OperateSphericalShellNormaliseProcedureNode>({});

/*
* Keywords
*/

keywords_.addTarget<ConfigurationKeyword>("Configuration", "Set target configuration for the module", targetConfiguration_);

keywords_.setOrganisation("Options", "Sites", "Specify sites defining the angle interaction A-B-C.");
keywords_.add<SpeciesSiteVectorKeyword>("SiteA", "Specify site(s) which represent 'A' in the interaction A-B-C",
selectA_->speciesSites(), selectA_->axesRequired());
keywords_.add<SpeciesSiteVectorKeyword>("SiteB", "Specify site(s) which represent 'B' in the interaction A-B-C",
selectB_->speciesSites(), selectB_->axesRequired());
keywords_.add<SpeciesSiteVectorKeyword>("SiteC", "Specify site(s) which represent 'C' in the interaction A-B-C",
selectC_->speciesSites(), selectC_->axesRequired());
keywords_.add<SpeciesSiteVectorKeyword>("SiteA", "Specify site(s) which represent 'A' in the interaction A-B-C", a_);
keywords_.add<SpeciesSiteVectorKeyword>("SiteB", "Specify site(s) which represent 'B' in the interaction A-B-C", b_);
keywords_.add<SpeciesSiteVectorKeyword>("SiteC", "Specify site(s) which represent 'C' in the interaction A-B-C", c_);

keywords_.setOrganisation("Options", "Ranges", "Ranges over which to bin quantities from the calculation.");
keywords_.add<Vec3DoubleKeyword>("RangeAB", "Range (min, max, binwidth) of A-B distance binning", rangeAB_,
Expand All @@ -167,16 +46,16 @@ AngleModule::AngleModule() : Module(ModuleTypes::Angle), analyser_(ProcedureNode

keywords_.setOrganisation("Export");
keywords_.add<FileAndFormatKeyword>("ExportAB", "File format and file name under which to save calculated A-B RDF data",
processAB_->exportFileAndFormat(), "EndExportAB");
exportFileAndFormatAB_, "EndExportAB");
keywords_.add<FileAndFormatKeyword>("ExportBC", "File format and file name under which to save calculated B-C RDF data",
processBC_->exportFileAndFormat(), "EndExportBC");
exportFileAndFormatBC_, "EndExportBC");
keywords_.add<FileAndFormatKeyword>("ExportAngle",
"File format and file name under which to save calculated A-B-C angle histogram",
processAngle_->exportFileAndFormat(), "EndExportAngle");
exportFileAndFormatAngle_, "EndExportAngle");
keywords_.add<FileAndFormatKeyword>("ExportDAngleAB",
"File format and file name under which to save calculated (A-B)-C distance-angle map",
processDAngleAB_->exportFileAndFormat(), "EndExportDAngleAB");
exportFileAndFormatDAngleAB_, "EndExportDAngleAB");
keywords_.add<FileAndFormatKeyword>("ExportDAngleBC",
"File format and file name under which to save calculated A-(B-C) distance-angle map",
processDAngleBC_->exportFileAndFormat(), "EndExportDAngleBC");
exportFileAndFormatDAngleBC_, "EndExportDAngleBC");
}

0 comments on commit 6692aab

Please sign in to comment.