diff --git a/Framework/Algorithms/CMakeLists.txt b/Framework/Algorithms/CMakeLists.txt index 3a44091c90bf..d2159bed58bb 100644 --- a/Framework/Algorithms/CMakeLists.txt +++ b/Framework/Algorithms/CMakeLists.txt @@ -87,6 +87,7 @@ set(SRC_FILES src/CreateGroupingWorkspace.cpp src/CreateLogPropertyTable.cpp src/CreateLogTimeCorrection.cpp + src/CreateMonteCarloWorkspace.cpp src/CreatePSDBleedMask.cpp src/CreatePeaksWorkspace.cpp src/CreateSampleWorkspace.cpp @@ -234,16 +235,16 @@ set(SRC_FILES src/Plus.cpp src/PointByPointVCorrection.cpp src/PoissonErrors.cpp - src/PolarizationCorrections/HeliumAnalyserEfficiency.cpp - src/PolarizationCorrections/PolarizationCorrectionsHelpers.cpp - src/PolarizationCorrections/SpinStateValidator.cpp - src/PolarizationCorrections/PolarizerEfficiency.cpp src/PolarizationCorrectionFredrikze.cpp src/PolarizationCorrectionWildes.cpp - src/PolarizationEfficiencyCor.cpp src/PolarizationCorrections/DepolarizedAnalyserTransmission.cpp src/PolarizationCorrections/FlipperEfficiency.cpp + src/PolarizationCorrections/HeliumAnalyserEfficiency.cpp + src/PolarizationCorrections/PolarizationCorrectionsHelpers.cpp src/PolarizationCorrections/PolarizationEfficienciesWildes.cpp + src/PolarizationCorrections/PolarizerEfficiency.cpp + src/PolarizationCorrections/SpinStateValidator.cpp + src/PolarizationEfficiencyCor.cpp src/PolynomialCorrection.cpp src/Power.cpp src/PowerLawCorrection.cpp @@ -437,6 +438,7 @@ set(INC_FILES inc/MantidAlgorithms/CreateGroupingWorkspace.h inc/MantidAlgorithms/CreateLogPropertyTable.h inc/MantidAlgorithms/CreateLogTimeCorrection.h + inc/MantidAlgorithms/CreateMonteCarloWorkspace.h inc/MantidAlgorithms/CreatePSDBleedMask.h inc/MantidAlgorithms/CreatePeaksWorkspace.h inc/MantidAlgorithms/CreateSampleWorkspace.h @@ -588,15 +590,15 @@ set(INC_FILES inc/MantidAlgorithms/Plus.h inc/MantidAlgorithms/PointByPointVCorrection.h inc/MantidAlgorithms/PoissonErrors.h + inc/MantidAlgorithms/PolarizationCorrectionFredrikze.h + inc/MantidAlgorithms/PolarizationCorrectionWildes.h inc/MantidAlgorithms/PolarizationCorrections/DepolarizedAnalyserTransmission.h + inc/MantidAlgorithms/PolarizationCorrections/FlipperEfficiency.h inc/MantidAlgorithms/PolarizationCorrections/HeliumAnalyserEfficiency.h inc/MantidAlgorithms/PolarizationCorrections/PolarizationCorrectionsHelpers.h - inc/MantidAlgorithms/PolarizationCorrections/SpinStateValidator.h - inc/MantidAlgorithms/PolarizationCorrections/FlipperEfficiency.h - inc/MantidAlgorithms/PolarizationCorrections/PolarizerEfficiency.h inc/MantidAlgorithms/PolarizationCorrections/PolarizationEfficienciesWildes.h - inc/MantidAlgorithms/PolarizationCorrectionFredrikze.h - inc/MantidAlgorithms/PolarizationCorrectionWildes.h + inc/MantidAlgorithms/PolarizationCorrections/PolarizerEfficiency.h + inc/MantidAlgorithms/PolarizationCorrections/SpinStateValidator.h inc/MantidAlgorithms/PolarizationEfficiencyCor.h inc/MantidAlgorithms/PolynomialCorrection.h inc/MantidAlgorithms/Power.h @@ -794,6 +796,7 @@ set(TEST_FILES CreateGroupingWorkspaceTest.h CreateLogPropertyTableTest.h CreateLogTimeCorrectionTest.h + CreateMonteCarloWorkspaceTest.h CreatePSDBleedMaskTest.h CreatePeaksWorkspaceTest.h CreateSampleWorkspaceTest.h @@ -937,15 +940,15 @@ set(TEST_FILES PlusTest.h PointByPointVCorrectionTest.h PoissonErrorsTest.h + PolarizationCorrectionFredrikzeTest.h + PolarizationCorrectionWildesTest.h PolarizationCorrections/DepolarizedAnalyserTransmissionTest.h + PolarizationCorrections/FlipperEfficiencyTest.h PolarizationCorrections/HeliumAnalyserEfficiencyTest.h PolarizationCorrections/PolarizationCorrectionsHelpersTest.h - PolarizationCorrections/SpinStateValidatorTest.h - PolarizationCorrections/FlipperEfficiencyTest.h - PolarizationCorrections/PolarizerEfficiencyTest.h PolarizationCorrections/PolarizationEfficienciesWildesTest.h - PolarizationCorrectionFredrikzeTest.h - PolarizationCorrectionWildesTest.h + PolarizationCorrections/PolarizerEfficiencyTest.h + PolarizationCorrections/SpinStateValidatorTest.h PolarizationEfficiencyCorTest.h PolynomialCorrectionTest.h PowerLawCorrectionTest.h diff --git a/Framework/Algorithms/inc/MantidAlgorithms/CreateMonteCarloWorkspace.h b/Framework/Algorithms/inc/MantidAlgorithms/CreateMonteCarloWorkspace.h new file mode 100644 index 000000000000..535a39ab2b7a --- /dev/null +++ b/Framework/Algorithms/inc/MantidAlgorithms/CreateMonteCarloWorkspace.h @@ -0,0 +1,40 @@ +// Mantid Repository : https://github.com/mantidproject/mantid +// +// Copyright © 2025 ISIS Rutherford Appleton Laboratory UKRI, +// NScD Oak Ridge National Laboratory, European Spallation Source, +// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS +// SPDX - License - Identifier: GPL - 3.0 + +#pragma once + +#include +#include + +#include "MantidAPI/Algorithm.h" +#include "MantidAlgorithms/DllConfig.h" + +namespace Mantid { +namespace Algorithms { +using namespace std; + +/** CreateMonteCarloWorkspace : The algorithm generates a simulated workspace by sampling from the probability +distribution of input data, useful for testing of fitting functions and modeling. + */ +class MANTID_ALGORITHMS_DLL CreateMonteCarloWorkspace : public API::Algorithm { +public: + const string name() const override; + int version() const override; + const string category() const override; + const string summary() const override; + + Mantid::HistogramData::HistogramY fillHistogramWithRandomData(const std::vector &cdf, int numIterations, + int seedInput, API::Progress &progress); + std::vector computeNormalizedCDF(const Mantid::HistogramData::HistogramY &yData); + int integrateYData(const Mantid::HistogramData::HistogramY &yData); + +private: + void init() override; + void exec() override; +}; + +} // namespace Algorithms +} // namespace Mantid diff --git a/Framework/Algorithms/src/CreateMonteCarloWorkspace.cpp b/Framework/Algorithms/src/CreateMonteCarloWorkspace.cpp new file mode 100644 index 000000000000..47ad958422e8 --- /dev/null +++ b/Framework/Algorithms/src/CreateMonteCarloWorkspace.cpp @@ -0,0 +1,163 @@ +// Mantid Repository : https://github.com/mantidproject/mantid +// +// Copyright © 2025 ISIS Rutherford Appleton Laboratory UKRI, +// NScD Oak Ridge National Laboratory, European Spallation Source, +// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS +// SPDX - License - Identifier: GPL - 3.0 + +#include +#include +#include +#include +#include // For std::accumulate +#include +#include + +#include "MantidAPI/MatrixWorkspace.h" +#include "MantidAPI/Progress.h" +#include "MantidAPI/Workspace.h" +#include "MantidAPI/WorkspaceFactory.h" +#include "MantidAlgorithms/CreateMonteCarloWorkspace.h" + +#include "MantidKernel/BoundedValidator.h" +#include "MantidKernel/Logger.h" + +namespace { +Mantid::Kernel::Logger g_log("CreateMonteCarloWorkspace"); +} +namespace Mantid { +namespace Algorithms { +using Mantid::Kernel::Direction; +using namespace Mantid::API; +using namespace std; + +// Register the algorithm into the AlgorithmFactory +DECLARE_ALGORITHM(CreateMonteCarloWorkspace) + +//---------------------------------------------------------------------------------------------- + +/// Algorithms name for identification. @see Algorithm::name +const std::string CreateMonteCarloWorkspace::name() const { return "CreateMonteCarloWorkspace"; } + +/// Algorithm's version for identification. @see Algorithm::version +int CreateMonteCarloWorkspace::version() const { return 1; } + +/// Algorithm's category for identification. @see Algorithm::category +const std::string CreateMonteCarloWorkspace::category() const { return "Simulation"; } + +/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary +const std::string CreateMonteCarloWorkspace::summary() const { + return "Creates a randomly simulated workspace by sampling from the probability distribution of input data."; +} + +//---------------------------------------------------------------------------------------------- +/** Initialize the algorithm's properties. + */ +void CreateMonteCarloWorkspace::init() { + auto mustBePositive = std::make_shared>(); + mustBePositive->setLower(0); + + declareProperty(std::make_unique>("InputWorkspace", "", Direction::Input), + "Input Workspace containing data to be simulated"); + declareProperty("Seed", 32, mustBePositive, + "Integer seed that initialises the random-number generator, for reproducibility"); + declareProperty("MonteCarloEvents", 0, mustBePositive, + "Number of Monte Carlo events to simulate. Defaults to integral of input workspace if 0."); + declareProperty(std::make_unique>("OutputWorkspace", "", Direction::Output), + "Name of output workspace."); +} + +//---------------------------------------------------------------------------------------------- + +Mantid::HistogramData::HistogramY CreateMonteCarloWorkspace::fillHistogramWithRandomData(const std::vector &cdf, + int numIterations, + int seedInput, + API::Progress &progress) { + + Mantid::HistogramData::HistogramY outputY(cdf.size(), 0.0); + std::mt19937 gen(seedInput); + std::uniform_real_distribution<> dis(0.0, 1.0); + + int progressInterval = std::max(1, numIterations / 100); // Update progress every 1% + + for (int i = 0; i < numIterations; ++i) { + double randomNum = dis(gen); + auto it = std::lower_bound(cdf.begin(), cdf.end(), randomNum); + size_t index = std::distance(cdf.begin(), it); + + if (index < outputY.size()) { + outputY[index] += 1.0; + } + + if (i % progressInterval == 0) { + progress.report("Generating random data..."); + } + } + return outputY; +} + +/** + * Compute a normalized CDF [0..1] from the given histogram data. + */ +std::vector CreateMonteCarloWorkspace::computeNormalizedCDF(const Mantid::HistogramData::HistogramY &yData) { + std::vector cdf(yData.size()); + std::partial_sum(yData.begin(), yData.end(), cdf.begin()); + double totalCounts = cdf.back(); + + if (totalCounts > 0.0) { + // Normalize the CDF + std::transform(cdf.begin(), cdf.end(), cdf.begin(), [totalCounts](double val) { return val / totalCounts; }); + } else { + g_log.warning("Total counts are zero; normalization skipped."); + } + return cdf; +} + +/** + * Determine how many iterations to use for MC sampling. + * If userMCEvents > 0, use that directly; otherwise use the integral of the input data. + */ +int CreateMonteCarloWorkspace::integrateYData(const Mantid::HistogramData::HistogramY &yData) { + double totalCounts = std::accumulate(yData.begin(), yData.end(), 0.0); + int iterations = static_cast(std::round(totalCounts)); + + if (iterations == 0) { + g_log.warning("Total counts in the input workspace round to 0. No Monte Carlo events will be generated."); + } + return iterations; +} + +//---------------------------------------------------------------------------------------------- +/** Execute the algorithm. + */ +void CreateMonteCarloWorkspace::exec() { + MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace"); + int seedInput = getProperty("Seed"); + int userMCEvents = getProperty("MonteCarloEvents"); + + const auto &originalYData = inputWs->y(0); // Counts in each bin + + int numIterations = (userMCEvents > 0) ? userMCEvents : integrateYData(originalYData); + + API::Progress progress(this, 0.0, 1.0, 101); + progress.report("Computing normalized CDF..."); + std::vector cdf = computeNormalizedCDF(originalYData); + + MatrixWorkspace_sptr outputWs = WorkspaceFactory::Instance().create(inputWs, 1); + outputWs->setSharedX(0, inputWs->sharedX(0)); + + // Fill the bins with random data, following the distribution in the CDF + Mantid::HistogramData::HistogramY outputY = fillHistogramWithRandomData(cdf, numIterations, seedInput, progress); + outputWs->mutableY(0) = outputY; + + // Calculate errors as the square root of the counts + Mantid::HistogramData::HistogramE outputE(outputY.size()); + std::transform(outputY.begin(), outputY.end(), outputE.begin(), [](double count) { return std::sqrt(count); }); + outputWs->mutableE(0) = outputE; + + g_log.warning("Only the first spectrum is being plotted."); + + setProperty("OutputWorkspace", outputWs); +} + +} // namespace Algorithms +} // namespace Mantid diff --git a/Framework/Algorithms/test/CreateMonteCarloWorkspaceTest.h b/Framework/Algorithms/test/CreateMonteCarloWorkspaceTest.h new file mode 100644 index 000000000000..f16ce1b3eeea --- /dev/null +++ b/Framework/Algorithms/test/CreateMonteCarloWorkspaceTest.h @@ -0,0 +1,156 @@ +#include "MantidAPI/MatrixWorkspace.h" +#include "MantidAlgorithms/CreateMonteCarloWorkspace.h" +#include "MantidDataObjects/Workspace2D.h" +#include "MantidFrameworkTestHelpers/WorkspaceCreationHelper.h" +#include "MantidHistogramData/Histogram.h" +#include + +#include +#include + +using namespace Mantid::Algorithms; +using namespace Mantid::API; +using namespace Mantid::HistogramData; + +using Mantid::Algorithms::CreateMonteCarloWorkspace; + +class CreateMonteCarloWorkspaceTest : public CxxTest::TestSuite { +public: + static CreateMonteCarloWorkspaceTest *createSuite() { return new CreateMonteCarloWorkspaceTest(); } + static void destroySuite(CreateMonteCarloWorkspaceTest *suite) { delete suite; } + + void test_Init() { + CreateMonteCarloWorkspace alg; + TS_ASSERT_THROWS_NOTHING(alg.initialize()); + TS_ASSERT(alg.isInitialized()); + } + + void test_integrateYData() { + CreateMonteCarloWorkspace alg; + Mantid::HistogramData::HistogramY yData = {1.0, 2.0, 3.0, 4.0}; + int iterations = alg.integrateYData(yData); + TS_ASSERT_EQUALS(iterations, 10); // Verify yData sums to 10 (1+2+3+4) + } + + void test_computeNormalizedCDF() { + CreateMonteCarloWorkspace alg; + Mantid::HistogramData::HistogramY yData = {1.0, 2.0, 3.0, 4.0}; + std::vector cdf = alg.computeNormalizedCDF(yData); + TS_ASSERT_EQUALS(cdf.size(), yData.size()); + TS_ASSERT_DELTA(cdf.back(), 1.0, 1e-6); // Check last element is normalized to 1.0 + } + + void test_fillHistogramWithRandomData() { + CreateMonteCarloWorkspace alg; + std::vector cdf = {0.1, 0.3, 0.6, 1.0}; + Mantid::API::Progress progress(nullptr, 0.0, 1.0, 1); // Dummy progress + Mantid::HistogramData::HistogramY outputY = alg.fillHistogramWithRandomData(cdf, 100, 32, progress); + + auto sumCounts = std::accumulate(outputY.begin(), outputY.end(), 0.0); + TS_ASSERT_EQUALS(sumCounts, 100); // Ensure total number of counts is correct + } + + void test_exec_with_custom_MCEvents() { + auto inputWS = createInputWorkspace(10, 5.0); // 10 bins, each bin has 5.0 + auto outputWS = runMonteCarloWorkspace(inputWS, 32, 100, "MonteCarloTest_CustomMC"); + + TS_ASSERT(outputWS); + const auto &outputY = outputWS->y(0); + auto sumOutput = std::accumulate(outputY.begin(), outputY.end(), 0.0); + TS_ASSERT_DELTA(sumOutput, 100.0, 1e-6); // Verify sum matches custom MC events + + removeWorkspace("MonteCarloTest_CustomMC"); + } + + void test_exec_without_custom_events() { + // Passing zero => use the input data's sum (10 bins * 5.0 = 50 total) + auto inputWS = createInputWorkspace(10, 5.0); + auto outputWS = runMonteCarloWorkspace(inputWS, 32, 0, "MonteCarloTest_Default"); + + TS_ASSERT(outputWS); + const auto &outputY = outputWS->y(0); + auto sumOutput = std::accumulate(outputY.begin(), outputY.end(), 0.0); + TS_ASSERT_DELTA(sumOutput, 50.0, 1e-6); // Sum matches input data's total counts + + removeWorkspace("MonteCarloTest_Default"); + } + + void test_reproducibility_with_seed() { + // Both run with the same seed and should produce identical Y values + auto inputWS = createInputWorkspace(10, 5.0); + + auto outputWS1 = runMonteCarloWorkspace(inputWS, 42, 0, "MonteCarloTest_WS1"); + auto outputWS2 = runMonteCarloWorkspace(inputWS, 42, 0, "MonteCarloTest_WS2"); + + TS_ASSERT(outputWS1); + TS_ASSERT(outputWS2); + + const auto &outputY1 = outputWS1->y(0); + const auto &outputY2 = outputWS2->y(0); + + TS_ASSERT_EQUALS(outputY1.size(), outputY2.size()); + for (size_t i = 0; i < outputY1.size(); ++i) { + TS_ASSERT_EQUALS(outputY1[i], outputY2[i]); + } + + removeWorkspace("MonteCarloTest_WS1"); + removeWorkspace("MonteCarloTest_WS2"); + } + + void test_error_calculation() { + // We fill the input with perfect squares to easily check sqrt in the result + auto inputWS = WorkspaceCreationHelper::create2DWorkspace(1, 10); + auto &yData = inputWS->mutableY(0); + yData = {1.0, 4.0, 9.0, 16.0, 25.0, 36.0, 49.0, 64.0, 81.0, 100.0}; + + auto outputWS = runMonteCarloWorkspace(inputWS, 32, 0, "MonteCarloTest_Error"); + TS_ASSERT(outputWS); + + const auto &outputY = outputWS->y(0); + const auto &outputE = outputWS->e(0); + + TS_ASSERT_EQUALS(outputY.size(), outputE.size()); + for (size_t i = 0; i < outputY.size(); ++i) { + TS_ASSERT_DELTA(outputE[i], std::sqrt(outputY[i]), 1e-6); + } + + removeWorkspace("MonteCarloTest_Error"); + } + +private: + static MatrixWorkspace_sptr createInputWorkspace(int numBins, double initialValue); + static MatrixWorkspace_sptr runMonteCarloWorkspace(const MatrixWorkspace_sptr &inputWS, int seed, int mcEvents, + const std::string &outputName); + static void removeWorkspace(const std::string &workspaceName); +}; + +// -- Helper method implementations below -- + +MatrixWorkspace_sptr CreateMonteCarloWorkspaceTest::createInputWorkspace(int numBins, double initialValue) { + auto ws = WorkspaceCreationHelper::create2DWorkspace(1, numBins); + auto &yData = ws->mutableY(0); + std::fill(yData.begin(), yData.end(), initialValue); + return ws; +} + +MatrixWorkspace_sptr CreateMonteCarloWorkspaceTest::runMonteCarloWorkspace(const MatrixWorkspace_sptr &inputWS, + int seed, int mcEvents, + const std::string &outputName) { + CreateMonteCarloWorkspace alg; + alg.initialize(); + alg.setProperty("InputWorkspace", inputWS); + alg.setProperty("Seed", seed); + alg.setProperty("MonteCarloEvents", mcEvents); + alg.setPropertyValue("OutputWorkspace", outputName); + + TS_ASSERT_THROWS_NOTHING(alg.execute()); + TS_ASSERT(alg.isExecuted()); + + using Mantid::API::AnalysisDataService; + return AnalysisDataService::Instance().retrieveWS(outputName); +} + +void CreateMonteCarloWorkspaceTest::removeWorkspace(const std::string &workspaceName) { + using Mantid::API::AnalysisDataService; + AnalysisDataService::Instance().remove(workspaceName); +} diff --git a/docs/source/algorithms/CreateMonteCarloWorkspace-v1.rst b/docs/source/algorithms/CreateMonteCarloWorkspace-v1.rst new file mode 100644 index 000000000000..94def9f3f8c5 --- /dev/null +++ b/docs/source/algorithms/CreateMonteCarloWorkspace-v1.rst @@ -0,0 +1,53 @@ +.. algorithm:: + +.. summary:: + +.. relatedalgorithms:: + +.. properties:: + +Description +----------- +The algorithm generates a simulated workspace by sampling from the probability +distribution of input data, useful for testing of fitting functions and modeling. + +Usage +----- + +**Example - CreateMonteCarloWorkspace** + +.. testcode:: Create simulation and compare + + # import mantid algorithms, numpy and matplotlib + from mantid.simpleapi import * + import matplotlib.pyplot as plt + import numpy as np + from mantid.api import FunctionFactory + + func = FunctionWrapper(FunctionFactory.createInitialized("name=BackToBackExponential,I=25000,A=0.06,B=0.015,X0=30000,S=30;name=FlatBackground,A0=50")) + # create input workspace + x = np.linspace(29650.0, 30500.0, 201) + y = func(x) + e = np.sqrt(y) + ws = CreateWorkspace(DataX=x, DataY=y, DataE=e, UnitX="TOF") + # call algorithm + ws_mc = CreateMonteCarloWorkspace(InputWorkspace=ws, Seed=0) + + fig, axes = plt.subplots(subplot_kw={'projection': 'mantid'}) + axes.plot(ws, label='input', wkspIndex=0) + axes.plot(ws_mc, label='CreateMonteCarloWorkspace output', wkspIndex=0, alpha=0.75) + legend = axes.legend(fontsize=8.0).set_draggable(True).legend + fig.show() + + +.. image:: ../../../images/CreateMonteCarloWorkspace_spectrum.png + :alt: Overplot of simulated data over input data + :width: 500px + :height: 400px + :scale: 100% + :align: center + :class: custom-class + +.. categories:: + +.. sourcelink:: diff --git a/docs/source/release/v6.12.0/Inelastic/Algorithms/New_features/38196_CreateMonteCarloWorkspace.rst b/docs/source/release/v6.12.0/Inelastic/Algorithms/New_features/38196_CreateMonteCarloWorkspace.rst new file mode 100644 index 000000000000..47217f56131d --- /dev/null +++ b/docs/source/release/v6.12.0/Inelastic/Algorithms/New_features/38196_CreateMonteCarloWorkspace.rst @@ -0,0 +1 @@ +- New algorithm :ref:`CreateMonteCarloWorkspace ` that creates a randomly simulated workspace by sampling from the probability distribution of input data. diff --git a/images/CreateMonteCarloWorkspace_spectrum.png b/images/CreateMonteCarloWorkspace_spectrum.png new file mode 100644 index 000000000000..7cefd398f01d Binary files /dev/null and b/images/CreateMonteCarloWorkspace_spectrum.png differ