Correct error propagation on PolarizationEfficienciesWildes algorithm (#39300)

### Description of work

#### Summary of work
This PR corrects the error propagation for the
`PolarizationEfficienciesWildes` algorithm.

As part of this work a class `ErrorPropagation` has been added to
`PolarizationCorrectionsHelpers.h`. This generic class utilizes the
eigen `AutoDiff` module to automatically calculate the partials for
different input variables. If there is wider interest in this class I
will move it into the `kernel` in a separate PR.

The `PolarizationEfficienciesWildes` algorithm has been refactored
slightly to enable use of the `ErrorPropagation` class and to remove the
use of derived quantities in workspace arithmetic.

The theory behind the approach to error propagation can be read: #37768

Fixes #37768

### To test:

1) Check all unit tests pass. As unit tests regarding the values of the
derived parameters were implemented, this proves that the calculation of
the values of these parameters is unchanged.
2) Check unit tests for the `ErrorPropagation` class in
`PolarizationCorrectionsHelpersTest`. These include simple examples
`testErrorPropagation_linear`, `testErrorPropagation_quad` and
`testErrorPropagation_mult_vars` that prove the function of the class.
You can check these on paper/using wolfram alpha. The latter test I
copied from this blog, obtaining the same results:
https://michaelspieler.com/post/eigen-autodiff/
3) Script to run `PolarizationEfficienciesWildes` on some real (system
test) data:
```
from mantid.simpleapi import *

Load(Filename=r'POLREF00032130.nxs', OutputWorkspace='POLREF00032130')
ConvertUnits(InputWorkspace='POLREF00032130', OutputWorkspace='POLREF00032130_y', Target='Wavelength', ConvertFromPointData=False)
Rebin(InputWorkspace='POLREF00032130_y', OutputWorkspace='POLREF00032130_y_rb', Params='0.01', PreserveEvents=False)
SumSpectra(InputWorkspace='POLREF00032130_y_rb', OutputWorkspace='POLREF00032130_y_rb_sum')
PolarizationEfficienciesWildes(InputNonMagWorkspace='POLREF00032130_y_rb_sum', OutputFpEfficiency='fp', OutputFaEfficiency='fa', IncludeDiagnosticOutputs=True)

```
4) Check error values output are sensible (however you determine that!)

Author: MialLewis

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