Use table of profiles passed as an argument

Previously the profiles were accessed through a list of bespoke dataclass
objects called NeutronComptonProfile.
Made the change to use a table workspace instead, to conform with the new
input argument.
Still need to change the way the Joint routine links profiles together.

Author: GuiMacielPereira

src/mvesuvio/analysis_reduction.py

@@ -52,6 +52,11 @@ def PyInit(self):
                                                direction=Direction.Input),
                                                doc="Workspace to fit Neutron Compton Profiles.")
 
+        self.declareProperty(TableWorkspaceProperty(name="InputProfiles",
+                                        defaultValue="",
+                                        direction=Direction.Input),
+                                        doc="Table workspace containing starting parameters for profiles")
+
         self.declareProperty(FileProperty(name='InstrumentParametersFile', 
                                       defaultValue='', 
                                       action=FileAction.Load, 
@@ -90,17 +95,14 @@ def PyInit(self):
         self.declareProperty("ResultsPath", "", doc="Directory to store results, to be deleted later")
         self.declareProperty("FiguresPath", "", doc="Directory to store figures, to be deleted later")
 
-        self.declareProperty(TableWorkspaceProperty(name="InputProfiles",
-                                                    defaultValue="",
-                                                    direction=Direction.Input),
-                                                    doc="Table workspace containing starting parameters for profiles")
 
     def PyExec(self):
         self._setup()
+        self.run()
 
     def _setup(self):
         self._name = self.getProperty("InputWorkspace").value.name()
-        self._ip_file = self.getPropertY("InstrumentParametersFile").value
+        self._ip_file = self.getProperty("InstrumentParametersFile").value
         self._number_of_iterations = self.getProperty("NumberOfIterations").value
         self._mask_spectra = self.getProperty("InvalidDetectors").value 
         self._transmission_guess = self.getProperty("TransmissionGuess").value 
@@ -133,76 +135,32 @@ def _setup(self):
         self._run_norm_voigt = False
 
 
-    def add_profiles(self, *args: NeutronComptonProfile):
-        for profile in args:
-            self._profiles[profile.label] = profile
-
-
-    @property
-    def profiles(self):
-        return self._profiles
-
-
-    def set_initial_profiles_from(self, source: 'AnalysisRoutine'):
-        
-        # Set intensities
-        for p in self._profiles.values():
-            if np.isclose(p.mass, 1, atol=0.1):    # Hydrogen present
-                p.intensity = source._h_ratio * source._get_lightest_profile().mean_intensity
-                continue
-            p.intensity = source.profiles[p.label].mean_intensity
-
-        # Normalise intensities
-        sum_intensities = sum([p.intensity for p in self._profiles.values()])
-        for p in self._profiles.values():
-            p.intensity /= sum_intensities
-            
-        # Set widths
-        for p in self._profiles.values():
-            try:
-                p.width = source.profiles[p.label].mean_width
-            except KeyError:
-                continue
-
-        # Fix all widths except lightest mass
-        for p in self._profiles.values():
-            if p == self._get_lightest_profile():
-                continue
-            p.width_bounds = [p.width, p.width]
-
-        return
-
-
-    def _get_lightest_profile(self):
-        profiles = [p for p in self._profiles.values()]
-        masses = [p.mass for p in self._profiles.values()]
-        return profiles[np.argmin(masses)]
-
-
     def calculate_h_ratio(self):
 
-        if not np.isclose(self._get_lightest_profile().mass, 1, atol=0.1):    # Hydrogen present
+        if not np.isclose(min(self._profiles_table.column("mass")), 1, atol=0.1):    # Hydrogen not present
             return None
 
-        # Hydrogen is present 
-        intensities = np.array([p.mean_intensity for p in self._profiles.values()])
-        masses = np.array([p.mass for p in self._profiles.values()])
+        # Hydrogen present 
+        # intensities = np.array([p.mean_intensity for p in self._profiles.values()])
+        intensities = self.mean_intensity_ratios
+        # masses = np.array([p.mass for p in self._profiles.values()])
+        masses = self._profiles_table.column("mass")
 
         sorted_intensities = intensities[np.argsort(masses)]
         return sorted_intensities[0] / sorted_intensities[1] 
 
 
-    @property
-    def profiles(self):
-        return self._profiles
+    # @property
+    # def profiles(self):
+    #     return self._profiles
 
 
-    @profiles.setter
-    def profiles(self, incoming_profiles):
-        assert(isinstance(incoming_profiles, dict))
-        common_keys = self._profiles.keys() & incoming_profiles.keys() 
-        common_keys_profiles = {k: incoming_profiles[k] for k in common_keys}
-        self._profiles = {**self._profiles, **common_keys_profiles}
+    # @profiles.setter
+    # def profiles(self, incoming_profiles):
+    #     assert(isinstance(incoming_profiles, dict))
+    #     common_keys = self._profiles.keys() & incoming_profiles.keys() 
+    #     common_keys_profiles = {k: incoming_profiles[k] for k in common_keys}
+    #     self._profiles = {**self._profiles, **common_keys_profiles}
 
 
     def _update_workspace_data(self):
@@ -216,13 +174,13 @@ def _update_workspace_data(self):
 
         self._set_up_kinematic_arrays()
 
-        self._fit_parameters = np.zeros((len(self._dataY), 3 * len(self._profiles) + 3))
+        self._fit_parameters = np.zeros((len(self._dataY), 3 * self._profiles_table.rowCount() + 3))
         self._row_being_fit = 0 
         self._table_fit_results = self._initialize_table_fit_parameters()
 
         # Initialise workspaces for fitted ncp 
         self._fit_profiles_workspaces = {}
-        for element, p in self._profiles.items():
+        for element in self._profiles_table.column("label"):
             self._fit_profiles_workspaces[element] = self._create_emtpy_ncp_workspace(f'_{element}_ncp')
         self._fit_profiles_workspaces['total'] = self._create_emtpy_ncp_workspace(f'_total_ncp')
 
@@ -239,10 +197,10 @@ def _initialize_table_fit_parameters(self):
         )
         table.setTitle("SciPy Fit Parameters")
         table.addColumn(type="float", name="Spectrum")
-        for key in self._profiles.keys():
-            table.addColumn(type="float", name=f"{key} Intensity")
-            table.addColumn(type="float", name=f"{key} Width")
-            table.addColumn(type="float", name=f"{key} Center ")
+        for label in self._profiles_table.column("label"):
+            table.addColumn(type="float", name=f"{label} Intensity")
+            table.addColumn(type="float", name=f"{label} Width")
+            table.addColumn(type="float", name=f"{label} Center ")
         table.addColumn(type="float", name="Normalised Chi2")
         table.addColumn(type="float", name="Number of Iterations")
         return table
@@ -256,8 +214,8 @@ def mean_widths(self):
     @mean_widths.setter
     def mean_widths(self, value):
         self._mean_widths = value
-        for i, p in enumerate(self._profiles.values()):
-            p.mean_width = self._mean_widths[i]
+        # for i, p in enumerate(self._profiles.values()):
+        #     p.mean_width = self._mean_widths[i]
         return
 
 
@@ -268,8 +226,8 @@ def mean_intensity_ratios(self):
     @mean_intensity_ratios.setter
     def mean_intensity_ratios(self, value):
         self._mean_intensity_ratios = value
-        for i, p in enumerate(self.profiles.values()):
-            p.mean_intensity = self._mean_intensity_ratios[i]
+        # for i, p in enumerate(self.profiles.values()):
+        #     p.mean_intensity = self._mean_intensity_ratios[i]
         return
 
 
@@ -294,7 +252,7 @@ def _set_up_kinematic_arrays(self):
 
     def run(self):
 
-        assert len(self.profiles) > 0, "Add profiles before attempting to run the routine!"
+        assert self._profiles_table.rowCount() > 0, "Need at least one profile to run the routine!"
 
         # Legacy code from Bootstrap
         # if self.runningSampleWS:
@@ -442,7 +400,7 @@ def convertDataXToYSpacesForEachMass(self, dataX, delta_Q, delta_E):
         delta_E = delta_E[np.newaxis, :, :]
 
         mN, Ef, en_to_vel, vf, hbar = loadConstants()
-        masses = np.array([p.mass for p in self._profiles.values()]).reshape(-1, 1, 1)
+        masses = np.array(self._profiles_table.column("mass")).reshape(-1, 1, 1)
 
         energyRecoil = np.square(hbar * delta_Q) / 2.0 / masses
         ySpacesForEachMass = (
@@ -475,7 +433,7 @@ def _save_plots(self):
         ax.legend()
 
         fileName = self._workspace_being_fit.name() + "_profiles_sum.pdf"
-        savePath = self._save_figures_path / fileName
+        savePath = self._save_figures_path + fileName
         plt.savefig(savePath, bbox_inches="tight")
         plt.close(fig)
         return
@@ -497,11 +455,11 @@ def _set_means_and_std(self):
         """Calculate mean widths and intensities from tableWorkspace"""
 
         fitParsTable = self._table_fit_results
-        widths = np.zeros((len(self._profiles), fitParsTable.rowCount()))
+        widths = np.zeros((self._profiles_table.rowCount(), fitParsTable.rowCount()))
         intensities = np.zeros(widths.shape)
-        for i, p in enumerate(self._profiles.values()):
-            widths[i] = fitParsTable.column(f"{p.label} Width")
-            intensities[i] = fitParsTable.column(f"{p.label} Intensity")
+        for i, label in enumerate(self._profiles_table.column("label")):
+            widths[i] = fitParsTable.column(f"{label} Width")
+            intensities[i] = fitParsTable.column(f"{label} Intensity")
         (
             meanWidths,
             stdWidths,
@@ -510,7 +468,7 @@ def _set_means_and_std(self):
         ) = self.calculateMeansAndStds(widths, intensities)
 
         assert (
-            len(meanWidths) == len(self._profiles) 
+            len(meanWidths) == self._profiles_table.rowCount()
         ), "Number of mean widths must match number of profiles!"
 
         self.mean_widths = meanWidths     # Use setter
@@ -534,15 +492,15 @@ def createMeansAndStdTableWS(self):
 
         print("\nCreated Table with means and std:")
         print("\nMass    Mean \u00B1 Std Widths    Mean \u00B1 Std Intensities\n")
-        for p, mean_width, std_width, mean_intensity, std_intensity in zip(
-            self._profiles.values(),
+        for label, mean_width, std_width, mean_intensity, std_intensity in zip(
+            self._profiles_table.column("label"),
             self._mean_widths,
             self._std_widths,
             self._mean_intensity_ratios,
             self._std_intensity_ratios,
         ):
-            meansTableWS.addRow([p.label, mean_width, std_width, mean_intensity, std_intensity])
-            print(f"{p.label:5s}  {mean_width:10.5f} \u00B1 {std_width:7.5f}  \
+            meansTableWS.addRow([label, mean_width, std_width, mean_intensity, std_intensity])
+            print(f"{label:5s}  {mean_width:10.5f} \u00B1 {std_width:7.5f}  \
                     {mean_intensity:10.5f} \u00B1 {std_intensity:7.5f}\n")
         return
 
@@ -613,17 +571,31 @@ def filterWidthsAndIntensities(self, widthsIn, intensitiesIn):
     def _fit_neutron_compton_profiles_to_row(self):
 
         if np.all(self._dataY[self._row_being_fit] == 0):
-            self._table_fit_results.addRow(np.zeros(3*len(self._profiles)+3))
+            self._table_fit_results.addRow(np.zeros(3*self._profiles_table.rowCount()+3))
             return
 
         # Need to transform profiles into parameter array for minimize
         initial_parameters = []
+        for intensity, width, center in zip(
+            self._profiles_table.column("intensity"),
+            self._profiles_table.column("width"),
+            self._profiles_table.column("center")
+        ):
+            initial_parameters.extend([intensity, width, center])
+
         bounds = []
-        for p in self._profiles.values():
-            for attr in ['intensity', 'width', 'center']:
-                initial_parameters.append(getattr(p, attr))
-            for attr in ['intensity_bounds', 'width_bounds', 'center_bounds']:
-                bounds.append(getattr(p, attr))
+        for intensity_bounds, width_bounds, center_bounds in zip(
+            self._profiles_table.column("intensity_bounds"),
+            self._profiles_table.column("width_bounds"),
+            self._profiles_table.column("center_bounds")
+        ):
+            bounds.extend([eval(intensity_bounds), eval(width_bounds), eval(center_bounds)])
+
+        # for p in self._profiles.values():
+        #     for attr in ['intensity', 'width', 'center']:
+        #         initial_parameters.append(getattr(p, attr))
+        #     for attr in ['intensity_bounds', 'width_bounds', 'center_bounds']:
+        #         bounds.append(getattr(p, attr))
 
         result = optimize.minimize(
             self.errorFunction,
@@ -649,7 +621,7 @@ def _fit_neutron_compton_profiles_to_row(self):
 
         # Pass fit profiles into workspaces
         individual_ncps = self._neutron_compton_profiles(fitPars)
-        for ncp, element in zip(individual_ncps, self._profiles.keys()):
+        for ncp, element in zip(individual_ncps, self._profiles_table.column("label")):
             self._fit_profiles_workspaces[element].dataY(self._row_being_fit)[:] = ncp
 
         self._fit_profiles_workspaces['total'].dataY(self._row_being_fit)[:] = np.sum(individual_ncps, axis=0)
@@ -683,7 +655,7 @@ def _neutron_compton_profiles(self, pars):
         intensities = pars[::3].reshape(-1, 1)
         widths = pars[1::3].reshape(-1, 1)
         centers = pars[2::3].reshape(-1, 1)
-        masses = np.array([p.mass for p in self._profiles.values()]).reshape(-1, 1)
+        masses = np.array(self._profiles_table.column("mass")).reshape(-1, 1)
 
         v0, E0, deltaE, deltaQ = self._kinematic_arrays[self._row_being_fit]
 
@@ -734,7 +706,7 @@ def kinematicsAtYCenters(self, centers):
 
 
     def calcGaussianResolution(self, centers):
-        masses = np.array([p.mass for p in self._profiles.values()]).reshape(-1, 1)
+        masses = np.array(self._profiles_table.column("mass")).reshape(-1, 1)
         v0, E0, delta_E, delta_Q = self.kinematicsAtYCenters(centers)
         det, plick, angle, T0, L0, L1 = self._instrument_params[self._row_being_fit]
         dE1, dTOF, dTheta, dL0, dL1, dE1_lorz = self._resolution_params[self._row_being_fit]
@@ -780,7 +752,7 @@ def calcGaussianResolution(self, centers):
 
 
     def calcLorentzianResolution(self, centers):
-        masses = np.array([p.mass for p in self._profiles.values()]).reshape(-1, 1)
+        masses = np.array(self._profiles_table.column("mass")).reshape(-1, 1)
         v0, E0, delta_E, delta_Q = self.kinematicsAtYCenters(centers)
         det, plick, angle, T0, L0, L1 = self._instrument_params[self._row_being_fit]
         dE1, dTOF, dTheta, dL0, dL1, dE1_lorz = self._resolution_params[self._row_being_fit]
@@ -949,7 +921,7 @@ def createSlabGeometry(self, wsNCPM):
 
 
     def calcMSCorrectionSampleProperties(self, meanWidths, meanIntensityRatios):
-        masses = [p.mass for p in self._profiles.values()]
+        masses = np.array(self._profiles_table.column("mass"))
 
         # If Backscattering mode and H is present in the sample, add H to MS properties
         if self._mode_running == "BACKWARD":
@@ -1040,7 +1012,7 @@ def create_gamma_workspaces(self):
 
 
     def calcGammaCorrectionProfiles(self, meanWidths, meanIntensityRatios):
-        masses = [p.mass for p in self._profiles.values()]
+        masses = np.array(self._profiles_table.column("mass"))
         profiles = ""
         for mass, width, intensity in zip(masses, meanWidths, meanIntensityRatios):
             profiles += (
@@ -1091,9 +1063,9 @@ def _set_results(self):
 
                 # Extract individual ncp
                 allNCP = []
-                for p in self._profiles.values():
+                for label in self._profiles_table.column("label"):
                     ncpWsToAppend = mtd[
-                        wsIterName + f"_{p.label}_ncp"
+                        wsIterName + f"_{label}_ncp"
                     ]
                     allNCP.append(ncpWsToAppend.extractY())
                 allNCP = np.swapaxes(np.array(allNCP), 0, 1)