Use Minus() to subtract masses from total ncp (#150)

* Use Minus to subtract masses from total ncp

Using Minus algorithm makes sure that workspaces are compatible
with each other, which is important for scientists when sometimes
they have to subtract one workspace from another.

It also improves readability of the code.

* Update system tests for y-space fit

Before overwriting the benchmarks I checked that my changes
did not fail the previous system tests that were in place.
Updated the inputs and the benchmarks of the system tests,
because the previous inputs had workspaces that had unmatching
units of Y and X.

.gitignore

@@ -2,6 +2,5 @@
 __pycache__/
 *.py[cod]
 
-src/vesuvio_analysis/core_functions/bootstrap_ws/
 figures/
-tests/analysis/data/inputs/sample_test/output_files/
+tests/**/output_files/

src/mvesuvio/analysis_fitting.py

@@ -23,18 +23,15 @@ def fitInYSpaceProcedure(yFitIC, IC, wsTOF):
         print(f"Workspace to fit {wsTOF} not found.")
         return
 
-    ncpForEachMass = extractNCPFromWorkspaces(wsTOF, IC)
     wsResSum, wsRes = calculateMantidResolutionFirstMass(IC, yFitIC, wsTOF)
 
-    wsTOFMass0 = subtractAllMassesExceptFirst(IC, wsTOF, ncpForEachMass)
+    wsTOFMass0 = subtract_profiles_except_lightest(IC, wsTOF)
 
     if yFitIC.subtractFSE:
         wsFSEMass0 = find_ws_name_fse_first_mass(IC)
         wsTOFMass0 = Minus(wsTOFMass0, wsFSEMass0, OutputWorkspace=wsTOFMass0.name() + "_fse")
 
-    wsJoY, wsJoYAvg = ySpaceReduction(
-        wsTOFMass0, IC.masses[0], yFitIC, ncpForEachMass[:, 0, :]
-    )
+    wsJoY, wsJoYAvg = ySpaceReduction(wsTOFMass0, IC.masses[0], yFitIC, IC)
 
     if yFitIC.symmetrisationFlag:
         wsJoYAvg = symmetrizeWs(wsJoYAvg)
@@ -69,34 +66,6 @@ def find_ws_name_fse_first_mass(ic):
     return ws_names_fse[np.argmin(ws_masses)]
 
 
-def extractNCPFromWorkspaces(wsFinal, ic):
-    """Extra function to extract ncps from loaded ws in mantid."""
-
-    for ws_name in mtd.getObjectNames():
-        if ws_name.startswith(ic.name+'_'+str(ic.noOfMSIterations)) and ws_name.endswith('ncp'):
-
-            if 'total' in ws_name:
-                continue
-
-            ws = mtd[ws_name]
-            dataY = ws.extractY()[np.newaxis, :, :]
-            try:
-                ncpForEachMass = np.append(ncpForEachMass, dataY, axis=0)
-            except UnboundLocalError:
-                ncpForEachMass = dataY 
-
-    # Ensure shape of ncp matches data
-    shape = ncpForEachMass.shape
-    assert shape[0] == len(ic.masses)
-    assert shape[1] == wsFinal.getNumberHistograms()
-    # Final dimension can be missing last col or not
-    assert (shape[2] == wsFinal.blocksize()) | (shape[2] == wsFinal.blocksize() - 1)
-
-    ncpForEachMass = switchFirstTwoAxis(ncpForEachMass)  # Organizes ncp by spectra
-    print("\nExtracted NCP profiles from workspaces.\n")
-    return ncpForEachMass
-
-
 def calculateMantidResolutionFirstMass(IC, yFitIC, ws):
     mass = IC.masses[0]
 
@@ -125,36 +94,26 @@ def calculateMantidResolutionFirstMass(IC, yFitIC, ws):
     return wsResSum, mtd[resName]
 
 
-def subtractAllMassesExceptFirst(ic, ws, ncpForEachMass):
+def subtract_profiles_except_lightest(ic, ws):
     """
     Isolates TOF data from first mass only.
     Input: ws containing TOF values, NCP for all mass profiles.
     Output: ws with all profiles except first subtracted.
     """
+    if len(ic.masses) == 1:
+        return
 
-    ncpForEachMass = switchFirstTwoAxis(ncpForEachMass)
-    # Select all masses other than the first one
-    ncpForEachMassExceptFirst = ncpForEachMass[1:, :, :]
-    # Sum the ncpTotal for remaining masses
-    ncpTotalExceptFirst = np.sum(ncpForEachMassExceptFirst, axis=0)
+    ws_name_lightest_mass = ic.name + '_' + str(ic.noOfMSIterations) + '_' + str(min(ic.masses)) + '_ncp'
+    ws_name_profiles = ic.name + '_' + str(ic.noOfMSIterations) + '_total_ncp'
 
-    dataX, dataY, dataE = extractWS(ws)
+    wsNcpExceptFirst = Minus(mtd[ws_name_profiles], mtd[ws_name_lightest_mass], 
+                             OutputWorkspace=ws_name_profiles + '_except_lightest')
 
-    # Adjust for last column missing or not
-    dataY[:, : ncpTotalExceptFirst.shape[1]] -= ncpTotalExceptFirst
+    SumSpectra(wsNcpExceptFirst.name(), OutputWorkspace=wsNcpExceptFirst.name() + "_sum")
 
-    # Ignore any masked bins (columns) from initial ws
-    mask = np.all(ws.extractY() == 0, axis=0)
-    dataY[:, mask] = 0
-
-    wsSubMass = CloneWorkspace(InputWorkspace=ws, OutputWorkspace=ws.name() + "_m0")
-    passDataIntoWS(dataX, dataY, dataE, wsSubMass)
-    wsMask, maskList = ExtractMask(ws)
-    MaskDetectors(Workspace=wsSubMass, MaskedWorkspace=wsMask)
-    SumSpectra(
-        InputWorkspace=wsSubMass.name(), OutputWorkspace=wsSubMass.name() + "_Sum"
-    )
-    return wsSubMass
+    wsFirstMass = Minus(ws, wsNcpExceptFirst, OutputWorkspace=ws.name()+"_m0")
+    SumSpectra(wsFirstMass.name(), OutputWorkspace=wsFirstMass.name() + "_sum")
+    return wsFirstMass
 
 
 def switchFirstTwoAxis(A):
@@ -164,9 +123,12 @@ def switchFirstTwoAxis(A):
     return np.stack(np.split(A, len(A), axis=0), axis=2)[0]
 
 
-def ySpaceReduction(wsTOF, mass0, yFitIC, ncp):
+def ySpaceReduction(wsTOF, mass0, yFitIC, ic):
     """Seperate procedures depending on masking specified."""
 
+    ws_name_lightest_mass = ic.name + '_' + str(ic.noOfMSIterations) + '_' + str(min(ic.masses)) + '_ncp'
+    ncp = mtd[ws_name_lightest_mass].extractY()
+
     rebinPars = yFitIC.rebinParametersForYSpaceFit
 
     if np.any(np.all(wsTOF.extractY() == 0, axis=0)):  # Masked columns present

tests/data/analysis/inputs/analysis_test.py

@@ -94,3 +94,10 @@ class YSpaceFitInputs:
     globalFit = None
     nGlobalFitGroups = 4
     maskTypeProcedure = None
+
+
+if (__name__ == "__main__") or (__name__ == "mantidqt.widgets.codeeditor.execution"):
+    import mvesuvio
+    from pathlib import Path
+    mvesuvio.set_config(inputs_file=Path(__file__))
+    mvesuvio.run()

tests/data/analysis/inputs/yspace_gauss_test.py

@@ -6,8 +6,15 @@
 ForwardAnalysisInputs.noOfMSIterations = 1
 ForwardAnalysisInputs.firstSpec = 164
 ForwardAnalysisInputs.lastSpec = 175
+ForwardAnalysisInputs.maskedSpecAllNo = [173, 174]
 ForwardAnalysisInputs.fit_in_y_space = True 
 BackwardAnalysisInputs.fit_in_y_space = False
 ForwardAnalysisInputs.run_this_scattering_type = True 
 BackwardAnalysisInputs.run_this_scattering_type = False
 YSpaceFitInputs.fitModel = "SINGLE_GAUSSIAN"
+
+if (__name__ == "__main__") or (__name__ == "mantidqt.widgets.codeeditor.execution"):
+    import mvesuvio
+    from pathlib import Path
+    mvesuvio.set_config(inputs_file=Path(__file__))
+    mvesuvio.run()

tests/data/analysis/inputs/yspace_gc_test.py

@@ -7,9 +7,17 @@
 ForwardAnalysisInputs.noOfMSIterations = 1
 ForwardAnalysisInputs.firstSpec = 164
 ForwardAnalysisInputs.lastSpec = 175
+ForwardAnalysisInputs.maskedSpecAllNo = [173, 174]
 ForwardAnalysisInputs.fit_in_y_space = True 
 BackwardAnalysisInputs.fit_in_y_space = False
 ForwardAnalysisInputs.run_this_scattering_type = True 
 BackwardAnalysisInputs.run_this_scattering_type = False
 YSpaceFitInputs.fitModel = "GC_C4_C6"
 YSpaceFitInputs.symmetrisationFlag = False 
+
+
+if (__name__ == "__main__") or (__name__ == "mantidqt.widgets.codeeditor.execution"):
+    import mvesuvio
+    from pathlib import Path
+    mvesuvio.set_config(inputs_file=Path(__file__))
+    mvesuvio.run()

tests/system/analysis/test_yspace_fit.py

@@ -1,5 +1,5 @@
 from mvesuvio.main.run_routine import Runner 
-from mantid.simpleapi import Load
+from mantid.simpleapi import Load, Plus, mtd, CreateWorkspace, CloneWorkspace
 from mantid.api import AnalysisDataService
 from pathlib import Path
 import numpy as np
@@ -43,17 +43,19 @@ def get_current_result(cls):
     @classmethod
     def _load_workspaces(cls):
         AnalysisDataService.clear()
-        scriptName = handle_config.get_script_name()
-        wsFinal = Load(
-            str(cls._input_data_path / "wsFinal.nxs"),
-            OutputWorkspace=scriptName + "_fwd_1",
+        Load(
+            str(cls._input_data_path / "yspace_tests_fwd_1.nxs"),
+            OutputWorkspace="yspace_gauss_test_fwd_1"
         )
-        for i in range(4):
-            fileName = "wsFinal_ncp_" + str(i) + ".nxs"
-            Load(
-                str(cls._input_data_path / fileName),
-                OutputWorkspace=wsFinal.name() + "_label" + str(i) +"_ncp",
-            )
+        Load(
+            str(cls._input_data_path / "yspace_tests_fwd_1_total_ncp.nxs"),
+            OutputWorkspace="yspace_gauss_test_fwd_1_total_ncp"
+        )
+        Load(
+            str(cls._input_data_path / "yspace_tests_fwd_1_1.0079_ncp.nxs"),
+            OutputWorkspace="yspace_gauss_test_fwd_1_1.0079_ncp"
+        )
+        return
 
     @classmethod
     def _run(cls):
@@ -63,7 +65,7 @@ def _run(cls):
 
     @classmethod
     def _load_benchmark_results(cls):
-        cls._benchmarkResults = np.load(str(cls._benchmark_path / "stored_yspace_fit.npz"))
+        cls._benchmarkResults = np.load(str(cls._benchmark_path / "yspace_gauss_test.npz"))
 
 
 class TestSymSumYSpace(unittest.TestCase):

tests/system/analysis/test_yspace_fit_GC.py

@@ -43,17 +43,19 @@ def get_current_result(cls):
     @classmethod
     def _load_workspaces(cls):
         AnalysisDataService.clear()
-        scriptName = handle_config.get_script_name()
-        wsFinal = Load(
-            str(cls._input_data_path / "wsFinal.nxs"),
-            OutputWorkspace=scriptName + "_fwd_1",
+        Load(
+            str(cls._input_data_path / "yspace_tests_fwd_1.nxs"),
+            OutputWorkspace="yspace_gc_test_fwd_1"
         )
-        for i in range(4):
-            fileName = "wsFinal_ncp_" + str(i) + ".nxs"
-            Load(
-                str(cls._input_data_path / fileName),
-                OutputWorkspace=wsFinal.name() + "_label" + str(i) +"_ncp",
-            )
+        Load(
+            str(cls._input_data_path / "yspace_tests_fwd_1_total_ncp.nxs"),
+            OutputWorkspace="yspace_gc_test_fwd_1_total_ncp"
+        )
+        Load(
+            str(cls._input_data_path / "yspace_tests_fwd_1_1.0079_ncp.nxs"),
+            OutputWorkspace="yspace_gc_test_fwd_1_1.0079_ncp"
+        )
+        return
 
     @classmethod
     def _run(cls):
@@ -64,7 +66,7 @@ def _run(cls):
     @classmethod
     def _load_benchmark_results(cls):
         cls._benchmarkResults = np.load(
-            str(cls._benchmark_path / "stored_yspace_fit_GC.npz")
+            str(cls._benchmark_path / "yspace_gc_test.npz")
         )