Adjust Algorithm and tests based on review

Updated Script and Image Update Framework/Algorithms/src/CreateMonteCarloWorkspace.cpp Co-authored-by: RichardWaiteSTFC <[email protected]> Update Style of code and image change copyright year to 2025 Update Framework/Algorithms/src/CreateMonteCarloWorkspace.cpp Co-authored-by: thomashampson <[email protected]> Update Framework/Algorithms/src/CreateMonteCarloWorkspace.cpp Co-authored-by: thomashampson <[email protected]> Update Framework/Algorithms/src/CreateMonteCarloWorkspace.cpp Co-authored-by: thomashampson <[email protected]> Rename Variables Rename Variables
mantidproject · Jan 30, 2025 · 9910e15 · 9910e15
1 parent c5e67ae
commit 9910e15
Show file tree

Hide file tree

Showing 6 changed files with 103 additions and 142 deletions.
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/CreateMonteCarloWorkspace.h b/Framework/Algorithms/inc/MantidAlgorithms/CreateMonteCarloWorkspace.h
@@ -1,6 +1,6 @@
 // Mantid Repository : https://github.com/mantidproject/mantid
 //
-// Copyright &copy; 2024 ISIS Rutherford Appleton Laboratory UKRI,
+// Copyright &copy; 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 +
@@ -29,11 +29,9 @@ class MANTID_ALGORITHMS_DLL CreateMonteCarloWorkspace : public API::Algorithm {
   const string summary() const override;
 
   Mantid::HistogramData::HistogramY fillHistogramWithRandomData(const std::vector<double> &cdf, int numIterations,
-                                                                int seed_input, API::Progress &progress);
+                                                                int seedInput, API::Progress &progress);
   std::vector<double> computeNormalizedCDF(const Mantid::HistogramData::HistogramY &yData);
-  int computeNumberOfIterations(const Mantid::HistogramData::HistogramY &yData, int userMCEvents);
-  Mantid::HistogramData::HistogramY scaleInputToMatchMCEvents(const Mantid::HistogramData::HistogramY &yData,
-                                                              int targetMCEvents);
+  int integrateYData(const Mantid::HistogramData::HistogramY &yData);
 
 private:
   void init() override;

diff --git a/Framework/Algorithms/src/CreateMonteCarloWorkspace.cpp b/Framework/Algorithms/src/CreateMonteCarloWorkspace.cpp
@@ -1,6 +1,6 @@
 // Mantid Repository : https://github.com/mantidproject/mantid
 //
-// Copyright &copy; 2024 ISIS Rutherford Appleton Laboratory UKRI,
+// Copyright &copy; 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 +
@@ -57,9 +57,9 @@ void CreateMonteCarloWorkspace::init() {
   mustBePositive->setLower(0);
 
   declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
-                  "Input Workspace containing data to be fitted");
+                  "Input Workspace containing data to be simulated");
   declareProperty("Seed", 32, mustBePositive,
-                  "Integer that initializes a random-number generator, good for reproducibility");
+                  "Integer that seeds 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<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
@@ -70,11 +70,11 @@ void CreateMonteCarloWorkspace::init() {
 
 Mantid::HistogramData::HistogramY CreateMonteCarloWorkspace::fillHistogramWithRandomData(const std::vector<double> &cdf,
                                                                                          int numIterations,
-                                                                                         int seed_input,
+                                                                                         int seedInput,
                                                                                          API::Progress &progress) {
 
   Mantid::HistogramData::HistogramY outputY(cdf.size(), 0.0);
-  std::mt19937 gen(seed_input);
+  std::mt19937 gen(seedInput);
   std::uniform_real_distribution<> dis(0.0, 1.0);
 
   int progressInterval = std::max(1, numIterations / 100); // Update progress every 1%
@@ -95,83 +95,58 @@ Mantid::HistogramData::HistogramY CreateMonteCarloWorkspace::fillHistogramWithRa
   return outputY;
 }
 
+/**
+ *  Compute a normalized CDF [0..1] from the given histogram data.
+ */
 std::vector<double> CreateMonteCarloWorkspace::computeNormalizedCDF(const Mantid::HistogramData::HistogramY &yData) {
   std::vector<double> cdf(yData.size());
   std::partial_sum(yData.begin(), yData.end(), cdf.begin());
-  double total_counts = cdf.back();
-  // Normalize the CDF
-  std::transform(cdf.begin(), cdf.end(), cdf.begin(), [total_counts](double val) { return val / total_counts; });
-  return cdf;
-}
+  double totalCounts = cdf.back();
 
-int CreateMonteCarloWorkspace::computeNumberOfIterations(const Mantid::HistogramData::HistogramY &yData,
-                                                         int userMCEvents) {
-  if (userMCEvents > 0) {
-    return userMCEvents;
+  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;
+}
 
-  // Default: Integral of Input
+/**
+ *  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 total_counts = std::accumulate(yData.begin(), yData.end(), 0.0);
   int iterations = static_cast<int>(std::round(total_counts));
 
   if (iterations == 0) {
     g_log.warning("Total counts in the input workspace round to 0. No Monte Carlo events will be generated.");
   }
-
   return iterations;
 }
 
-Mantid::HistogramData::HistogramY
-CreateMonteCarloWorkspace::scaleInputToMatchMCEvents(const Mantid::HistogramData::HistogramY &yData,
-                                                     int targetMCEvents) {
-  double total_counts = std::accumulate(yData.begin(), yData.end(), 0.0);
-  if (total_counts == 0) {
-    g_log.warning("Total counts in the input workspace are 0. Scaling cannot be performed.");
-    return yData;
-  }
-
-  double scaleFactor = static_cast<double>(targetMCEvents) / total_counts;
-  Mantid::HistogramData::HistogramY scaledY(yData.size());
-
-  std::transform(yData.begin(), yData.end(), scaledY.begin(),
-                 [scaleFactor](double count) { return count * scaleFactor; });
-
-  return scaledY;
-}
-
 //----------------------------------------------------------------------------------------------
 /** Execute the algorithm.
  */
-
-// Using Cumulative Distribution Function
 void CreateMonteCarloWorkspace::exec() {
   MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace");
-  int seed_input = getProperty("Seed");
+  int seedInput = getProperty("Seed");
   int userMCEvents = getProperty("MonteCarloEvents");
 
-  MatrixWorkspace_sptr scaledWs = inputWs;
+  const auto &originalYData = inputWs->y(0); // Counts in each bin
 
-  const Mantid::HistogramData::HistogramY &originalYData = inputWs->y(0); // Counts in each bin
-  Mantid::HistogramData::HistogramY yData = originalYData;
+  int numIterations = (userMCEvents > 0) ? userMCEvents : integrateYData(originalYData);
 
-  // Scale input workspace if user has specified MC events
-  if (userMCEvents > 0) {
-    g_log.warning() << "Custom Monte Carlo events: " << userMCEvents << ". Input workspace scaled accordingly.\n";
-    yData = scaleInputToMatchMCEvents(originalYData, userMCEvents);
-  }
-
-  int numIterations = computeNumberOfIterations(yData, userMCEvents);
-
-  std::vector<double> cdf = computeNormalizedCDF(yData);
-  // Set up progress bar
   API::Progress progress(this, 0.0, 1.0, 2);
   progress.report("Computing normalized CDF...");
+  std::vector<double> cdf = computeNormalizedCDF(originalYData);
 
   MatrixWorkspace_sptr outputWs = WorkspaceFactory::Instance().create(inputWs, 1);
-
-  // Copy the bin boundaries (X-values) from the input to the output
   outputWs->setSharedX(0, inputWs->sharedX(0));
-  Mantid::HistogramData::HistogramY outputY = fillHistogramWithRandomData(cdf, numIterations, seed_input, progress);
+
+  // 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

diff --git a/Framework/Algorithms/test/CreateMonteCarloWorkspaceTest.h b/Framework/Algorithms/test/CreateMonteCarloWorkspaceTest.h
@@ -5,6 +5,7 @@
 #include "MantidHistogramData/Histogram.h"
 #include <cxxtest/TestSuite.h>
 
+#include <numeric>
 #include <random>
 
 using namespace Mantid::Algorithms;
@@ -24,23 +25,19 @@ class CreateMonteCarloWorkspaceTest : public CxxTest::TestSuite {
     TS_ASSERT(alg.isInitialized());
   }
 
-  void test_computeNumberOfIterations() {
+  void test_integrateYData() {
     CreateMonteCarloWorkspace alg;
     Mantid::HistogramData::HistogramY yData = {1.0, 2.0, 3.0, 4.0};
-    int iterations = alg.computeNumberOfIterations(yData, 0);
-    TS_ASSERT_EQUALS(iterations, 10); // Verify yData rounds correctly
-
-    // Test custom Monte Carlo events
-    iterations = alg.computeNumberOfIterations(yData, 20);
-    TS_ASSERT_EQUALS(iterations, 20); // Should override with the custom number
+    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<double> cdf = alg.computeNormalizedCDF(yData);
     TS_ASSERT_EQUALS(cdf.size(), yData.size());
-    TS_ASSERT_DELTA(cdf.back(), 1.0, 1e-6); // Check the last element is normalized to 1.0
+    TS_ASSERT_DELTA(cdf.back(), 1.0, 1e-6); // Check last element is normalized to 1.0
   }
 
   void test_fillHistogramWithRandomData() {
@@ -50,81 +47,37 @@ class CreateMonteCarloWorkspaceTest : public CxxTest::TestSuite {
     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 the total number of counts is correct
+    TS_ASSERT_EQUALS(sumCounts, 100); // Ensure total number of counts is correct
   }
 
-  void test_scaleInputToMatchMCEvents() {
-    CreateMonteCarloWorkspace alg;
-    Mantid::HistogramData::HistogramY yData = {1.0, 2.0, 3.0, 4.0};
-
-    // Test scaling
-    int targetMCEvents = 20;
-    Mantid::HistogramData::HistogramY scaledY = alg.scaleInputToMatchMCEvents(yData, targetMCEvents);
-
-    double totalScaledCounts = std::accumulate(scaledY.begin(), scaledY.end(), 0.0);
-    TS_ASSERT_DELTA(totalScaledCounts, 20.0, 1e-6); // Verify the scaled sum matches targetMCEvents
-  }
-
-  MatrixWorkspace_sptr 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 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<MatrixWorkspace>(outputName);
-  }
-
-  void removeWorkspace(const std::string &workspaceName) {
-    using Mantid::API::AnalysisDataService;
-    AnalysisDataService::Instance().remove(workspaceName);
-  }
-
-  void test_exec_with_scaling() {
-    auto inputWS = createInputWorkspace(10, 5.0); // 10 bins, initial value 5.0
-    auto outputWS = runMonteCarloWorkspace(inputWS, 32, 100, "MonteCarloTest_Scaled");
+  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);
-
-    // Verify the output data
     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
 
-    // Clean up
-    removeWorkspace("MonteCarloTest_Scaled");
+    removeWorkspace("MonteCarloTest_CustomMC");
   }
 
-  void test_exec_without_scaling() {
-    auto inputWS = createInputWorkspace(10, 5.0); // 10 bins, initial value 5.0
+  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);
-
-    // Verify the output data
     const auto &outputY = outputWS->y(0);
     auto sumOutput = std::accumulate(outputY.begin(), outputY.end(), 0.0);
-    TS_ASSERT_DELTA(sumOutput, 50.0, 1e-6); // Verify sum matches input data's total counts
+    TS_ASSERT_DELTA(sumOutput, 50.0, 1e-6); // Sum matches input data's total counts
 
-    // Clean up
     removeWorkspace("MonteCarloTest_Default");
   }
 
   void test_reproducibility_with_seed() {
-    auto inputWS = createInputWorkspace(10, 5.0); // 10 bins, initial value 5.0
+    // 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");
@@ -140,16 +93,14 @@ class CreateMonteCarloWorkspaceTest : public CxxTest::TestSuite {
       TS_ASSERT_EQUALS(outputY1[i], outputY2[i]);
     }
 
-    // Clean up
     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);
-
-    // Using perfect squares (for sqrt testing)
     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");
@@ -163,7 +114,43 @@ class CreateMonteCarloWorkspaceTest : public CxxTest::TestSuite {
       TS_ASSERT_DELTA(outputE[i], std::sqrt(outputY[i]), 1e-6);
     }
 
-    // Clean up
     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<MatrixWorkspace>(outputName);
+}
+
+void CreateMonteCarloWorkspaceTest::removeWorkspace(const std::string &workspaceName) {
+  using Mantid::API::AnalysisDataService;
+  AnalysisDataService::Instance().remove(workspaceName);
+}