Correct error propagation on PolarizationEfficienciesWildes algorithm

Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrections/PolarizationCorrectionsHelpers.h

@@ -7,10 +7,14 @@
 
 #pragma once
 
+#include "MantidAPI/AlgorithmManager.h"
 #include "MantidAPI/MatrixWorkspace.h"
 #include "MantidAPI/WorkspaceGroup.h"
 #include "MantidAlgorithms/DllConfig.h"
+#include "MantidKernel/MultiThreaded.h"
+#include <Eigen/Dense>
 #include <optional>
+#include <unsupported/Eigen/AutoDiff>
 
 namespace Mantid::Algorithms {
 namespace PolarizationCorrectionsHelpers {
@@ -66,4 +70,109 @@ MANTID_ALGORITHMS_DLL const std::string &getORSONotationForSpinState(const std::
 MANTID_ALGORITHMS_DLL void addORSOLogForSpinState(const Mantid::API::MatrixWorkspace_sptr &ws,
                                                   const std::string &spinState);
 } // namespace SpinStatesORSO
+
+namespace Arithmetic {
+
+template <size_t N> class ErrorTypeHelper {
+public:
+  using DerType = Eigen::Matrix<double, N, 1>;
+  using InputArray = DerType;
+  using ADScalar = Eigen::AutoDiffScalar<DerType>;
+};
+
+template <size_t N, typename Func> class ErrorPropagation {
+public:
+  using Types = ErrorTypeHelper<N>;
+  using DerType = Types::DerType;
+  using ADScalar = Types::ADScalar;
+  using InputArray = Types::InputArray;
+  ErrorPropagation(Func func) : computeFunc(std::move(func)) {}
+
+  struct AutoDevResult {
+    double value;
+    double error;
+    Eigen::Array<double, N, 1> derivatives;
+  };
+
+  AutoDevResult evaluate(const InputArray &values, const InputArray &errors) const {
+    std::array<ADScalar, N> x;
+    for (size_t i = 0; i < N; ++i) {
+      x[i] = ADScalar(values[i], DerType::Unit(N, i));
+    }
+    const ADScalar y = computeFunc(x);
+    const auto &derivatives = y.derivatives();
+    return {y.value(), std::sqrt((derivatives.array().square() * errors.array().square()).sum()), derivatives};
+  }
+
+  template <std::same_as<API::MatrixWorkspace_sptr>... Ts>
+  API::MatrixWorkspace_sptr evaluateWorkspaces(const bool outputWorkspaceDistribution, Ts... args) const {
+    return evaluateWorkspacesImpl(outputWorkspaceDistribution, std::forward<Ts>(args)...);
+  }
+
+  template <std::same_as<API::MatrixWorkspace_sptr>... Ts>
+  API::MatrixWorkspace_sptr evaluateWorkspaces(Ts... args) const {
+    return evaluateWorkspacesImpl(std::nullopt, std::forward<Ts>(args)...);
+  }
+
+private:
+  Func computeFunc;
+
+  template <std::same_as<API::MatrixWorkspace_sptr>... Ts>
+  API::MatrixWorkspace_sptr evaluateWorkspacesImpl(std::optional<bool> outputWorkspaceDistribution, Ts... args) const {
+    const auto firstWs = std::get<0>(std::forward_as_tuple(args...));
+    API::MatrixWorkspace_sptr outWs = firstWs->clone();
+
+    if (outWs->id() == "EventWorkspace") {
+      outWs = convertToWorkspace2D(outWs);
+    }
+
+    const size_t numSpec = outWs->getNumberHistograms();
+    const size_t specSize = outWs->blocksize();
+
+    // cppcheck-suppress unreadVariable
+    const bool isThreadSafe = Kernel::threadSafe((*args)..., *outWs);
+    // cppcheck-suppress unreadVariable
+    const bool specOverBins = numSpec > specSize;
+
+    PARALLEL_FOR_IF(isThreadSafe && specOverBins)
+    for (int64_t i = 0; i < static_cast<int64_t>(numSpec); i++) {
+      auto &yOut = outWs->mutableY(i);
+      auto &eOut = outWs->mutableE(i);
+
+      PARALLEL_FOR_IF(isThreadSafe && !specOverBins)
+      for (int64_t j = 0; j < static_cast<int64_t>(specSize); ++j) {
+        const auto result = evaluate(InputArray{args->y(i)[j]...}, InputArray(args->e(i)[j]...));
+        yOut[j] = result.value;
+        eOut[j] = result.error;
+      }
+    }
+
+    if (outputWorkspaceDistribution.has_value()) {
+      outWs->setDistribution(outputWorkspaceDistribution.value());
+    }
+    return outWs;
+  }
+
+  API::MatrixWorkspace_sptr runWorkspaceConversionAlg(const API::MatrixWorkspace_sptr &workspace,
+                                                      const std::string &algName) const {
+    auto conversionAlg = API::AlgorithmManager::Instance().create(algName);
+    conversionAlg->initialize();
+    conversionAlg->setChild(true);
+    conversionAlg->setProperty("InputWorkspace", workspace);
+    conversionAlg->setProperty("OutputWorkspace", workspace->getName());
+    conversionAlg->execute();
+    return conversionAlg->getProperty("OutputWorkspace");
+  }
+
+  API::MatrixWorkspace_sptr convertToWorkspace2D(const API::MatrixWorkspace_sptr &workspace) const {
+    runWorkspaceConversionAlg(workspace, "ConvertToHistogram");
+    return runWorkspaceConversionAlg(workspace, "ConvertToMatrixWorkspace");
+  }
+};
+
+template <size_t N, typename Func> auto makeErrorPropagation(Func &&func) {
+  return ErrorPropagation<N, std::decay_t<Func>>(std::forward<Func>(func));
+}
+
+} // namespace Arithmetic
 } // namespace Mantid::Algorithms

Framework/Algorithms/inc/MantidAlgorithms/PolarizationCorrections/PolarizationEfficienciesWildes.h

@@ -10,6 +10,7 @@
 #include "MantidAPI/MatrixWorkspace.h"
 #include "MantidAPI/WorkspaceGroup.h"
 #include "MantidAlgorithms/DllConfig.h"
+#include <unordered_map>
 
 namespace Mantid::Algorithms {
 
@@ -45,19 +46,12 @@ class MANTID_ALGORITHMS_DLL PolarizationEfficienciesWildes : public API::Algorit
   /// Calculate Fp, Fa and Phi
   void calculateFlipperEfficienciesAndPhi();
 
-  /// Calculate (2p-1) from Phi, Fp, Fa and the magnetic workspace intensities
-  MatrixWorkspace_sptr calculateTPMOFromPhi(const WorkspaceGroup_sptr &magWsGrp);
+  /// Calculate (2p-1) from intensities
+  MatrixWorkspace_sptr calculateTPMO();
 
   /// Calculate the polarizer and/or analyser efficiencies, as requested
   void calculatePolarizerAndAnalyserEfficiencies(const bool solveForP, const bool solveForA);
 
-  /// If either the polarizer or the analyser efficiency is known, use the relationship Phi = (2p-1)(2a-1) to solve for
-  /// the other efficiency
-  MatrixWorkspace_sptr solveForUnknownEfficiency(const MatrixWorkspace_sptr &knownEfficiency);
-
-  /// Solve for the unknown efficiency from either (2p-1) or (2a-1) using the relationship Phi = (2p-1)(2a-1)
-  MatrixWorkspace_sptr solveUnknownEfficiencyFromTXMO(const MatrixWorkspace_sptr &wsTXMO);
-
   /// Set the algorithm outputs
   void setOutputs();
 
@@ -68,10 +62,29 @@ class MANTID_ALGORITHMS_DLL PolarizationEfficienciesWildes : public API::Algorit
   /// held by the property
   void resetPropertyValue(const std::string &propertyName);
 
+  /// Populates the spin state workspaces map from ws group given key prefix.
+  void populateSpinStateWorkspaces(const WorkspaceGroup_sptr &wsGrp, const std::string &keyPrefix = "");
+
+  /// Populates the spin state workspaces map
+  void mapSpinStateWorkspaces();
+
+  // Convenience struct for handling of flipper workspaces
+  struct FlipperWorkspaces {
+    const MatrixWorkspace_sptr &ws00;
+    const MatrixWorkspace_sptr &ws01;
+    const MatrixWorkspace_sptr &ws10;
+    const MatrixWorkspace_sptr &ws11;
+  };
+
+  /// Access flipper workspaces in the spin state workspaces map
+  FlipperWorkspaces getFlipperWorkspaces(const bool mag = false);
+
   MatrixWorkspace_sptr m_wsFp = nullptr;
   MatrixWorkspace_sptr m_wsFa = nullptr;
   MatrixWorkspace_sptr m_wsPhi = nullptr;
   MatrixWorkspace_sptr m_wsP = nullptr;
   MatrixWorkspace_sptr m_wsA = nullptr;
+  std::unordered_map<std::string, MatrixWorkspace_sptr> m_spinStateWorkspaces;
+  bool m_magWsProvided = false;
 };
 } // namespace Mantid::Algorithms