Merge pull request #133 from mantidproject/update-system-benchmark

Update benchmark for analysis system test

Author: GuiMacielPereira

src/mvesuvio/analysis_reduction.py

@@ -1,14 +1,13 @@
 import numpy as np 
 import matplotlib.pyplot as plt
-from scipy import optimize
+import scipy
 from mantid.simpleapi import mtd, CreateEmptyTableWorkspace, SumSpectra, \
                             CloneWorkspace, DeleteWorkspace, VesuvioCalculateGammaBackground, \
                             VesuvioCalculateMS, Scale, RenameWorkspace, Minus, CreateSampleShape, \
                             VesuvioThickness, Integration, Divide, Multiply, DeleteWorkspaces, \
                             CreateWorkspace
 
-from mvesuvio.util.analysis_helpers import histToPointData, loadConstants, \
-                                        gaussian, lorentzian, numericalThirdDerivative
+from mvesuvio.util.analysis_helpers import loadConstants, numericalThirdDerivative
 
 from dataclasses import dataclass
 
@@ -69,10 +68,6 @@ def __init__(self, workspace, ip_file, h_ratio_to_lowest_mass, number_of_iterati
         self._table_fit_results = None
         self._fit_profiles_workspaces = {}
 
-        # Only used for system tests, remove once tests are updated
-        self._run_hist_data = True 
-        self._run_norm_voigt = False
-
 
     def add_profiles(self, *args: NeutronComptonProfile):
         for profile in args:
@@ -152,9 +147,6 @@ def _update_workspace_data(self):
         self._dataY = self._workspace_being_fit.extractY()
         self._dataE = self._workspace_being_fit.extractE()
 
-        if self._run_hist_data:  # Converts point data from workspaces to histogram data
-            self._dataY, self._dataX, self._dataE = histToPointData(self._dataY, self._dataX, self._dataE)
-
         self._set_up_kinematic_arrays()
 
         self._fit_parameters = np.zeros((len(self._dataY), 3 * len(self._profiles) + 3))
@@ -593,7 +585,7 @@ def _fit_neutron_compton_profiles_to_row(self):
             for attr in ['intensity_bounds', 'width_bounds', 'center_bounds']:
                 bounds.append(getattr(p, attr))
 
-        result = optimize.minimize(
+        result = scipy.optimize.minimize(
             self.errorFunction,
             initial_parameters,
             method="SLSQP",
@@ -658,7 +650,7 @@ def _neutron_compton_profiles(self, pars):
         gaussRes, lorzRes = self.caculateResolutionForEachMass(centers)
         totalGaussWidth = np.sqrt(widths**2 + gaussRes**2)
 
-        JOfY = self.pseudoVoigt(self._y_space_arrays[self._row_being_fit] - centers, totalGaussWidth, lorzRes)
+        JOfY = scipy.special.voigt_profile(self._y_space_arrays[self._row_being_fit] - centers, totalGaussWidth, lorzRes)
 
         FSE = (
             -numericalThirdDerivative(self._y_space_arrays[self._row_being_fit], JOfY)
@@ -771,20 +763,6 @@ def calcLorentzianResolution(self, centers):
         return lorentzianResWidth
 
 
-    def pseudoVoigt(self, x, sigma, gamma):
-        """Convolution between Gaussian with std sigma and Lorentzian with HWHM gamma"""
-        fg, fl = 2.0 * sigma * np.sqrt(2.0 * np.log(2.0)), 2.0 * gamma
-        f = 0.5346 * fl + np.sqrt(0.2166 * fl**2 + fg**2)
-        eta = 1.36603 * fl / f - 0.47719 * (fl / f) ** 2 + 0.11116 * (fl / f) ** 3
-        sigma_v, gamma_v = f / (2.0 * np.sqrt(2.0 * np.log(2.0))), f / 2.0
-        pseudo_voigt = eta * lorentzian(x, gamma_v) + (1.0 - eta) * gaussian(x, sigma_v)
-
-        norm = (
-            np.abs(np.trapz(pseudo_voigt, x, axis=1))[:, np.newaxis] if self._run_norm_voigt else 1
-        )
-        return pseudo_voigt / norm
-
-
     def _get_parsed_constraints(self):
 
         parsed_constraints = []
@@ -836,8 +814,7 @@ def _replace_zero_columns_with_ncp_fit(self):
             OutputWorkspace=self._workspace_for_corrections.name() + "_CorrectionsInput"
         )
         for row in range(self._workspace_for_corrections.getNumberHistograms()):
-            # TODO: Once the option to change point to hist is removed, remove [:len(ncp)]
-            self._workspace_for_corrections.dataY(row)[self._zero_columns_boolean_mask] = ncp[row, self._zero_columns_boolean_mask[:len(ncp[row])]]
+            self._workspace_for_corrections.dataY(row)[self._zero_columns_boolean_mask] = ncp[row, self._zero_columns_boolean_mask]
 
         SumSpectra(
             InputWorkspace=self._workspace_for_corrections.name(), 
@@ -1082,7 +1059,6 @@ def _set_results(self):
         self.all_spec_best_par_chi_nit = np.array(allBestPar)
         self.all_tot_ncp = np.array(allTotNcp)
         self.all_ncp_for_each_mass = np.array(allIterNcp)
-
         self.all_mean_widths = np.array(allMeanWidhts)
         self.all_mean_intensities = np.array(allMeanIntensities)
         self.all_std_widths = np.array(allStdWidths)
@@ -1091,14 +1067,6 @@ def _set_results(self):
     def _save_results(self):
         """Saves all of the arrays stored in this object"""
 
-        maskedDetectorIdx = np.array(self._mask_spectra) - min(self._workspace_being_fit.getSpectrumNumbers())
-
-        # TODO: Take out nans next time when running original results
-        # Because original results were recently saved with nans, mask spectra with nans
-        self.all_spec_best_par_chi_nit[:, maskedDetectorIdx, :] = np.nan
-        self.all_ncp_for_each_mass[:, maskedDetectorIdx, :, :] = np.nan
-        self.all_tot_ncp[:, maskedDetectorIdx, :] = np.nan
-
         if not self._save_results_path:
             return 
 

src/mvesuvio/util/analysis_helpers.py

@@ -112,25 +112,6 @@ def extractWS(ws):
     return ws.extractX(), ws.extractY(), ws.extractE()
 
 
-def histToPointData(dataY, dataX, dataE):
-    """
-    Used only when comparing with original results.
-    Sets each dataY point to the center of bins.
-    Last column of data is removed.
-    Removed original scaling by bin widths
-    """
-
-    histWidths = dataX[:, 1:] - dataX[:, :-1]
-    assert np.min(histWidths) == np.max(
-        histWidths
-    ), "Histogram widths need to be the same length"
-
-    dataYp = dataY[:, :-1]
-    dataEp = dataE[:, :-1]
-    dataXp = dataX[:, :-1] + histWidths[0, 0] / 2
-    return dataYp, dataXp, dataEp
-
-
 def loadConstants():
     """Output: the mass of the neutron, final energy of neutrons (selected by gold foil),
     factor to change energies into velocities, final velocity of neutron and hbar"""
@@ -143,18 +124,6 @@ def loadConstants():
     return constants
 
 
-def gaussian(x, sigma):
-    """Gaussian centered at zero"""
-    gauss = np.exp(-(x**2) / 2 / sigma**2) 
-    gauss /= np.sqrt(2.0 * np.pi) * sigma
-    return gauss
-
-
-def lorentzian(x, gamma):
-    """Lorentzian centered at zero"""
-    return gamma / np.pi / (x**2 + gamma**2)
-
-
 def numericalThirdDerivative(x, y):
     k6 = (- y[:, 12:] + y[:, :-12]) * 1
     k5 = (+ y[:, 11:-1] - y[:, 1:-11]) * 24

tests/data/analysis/benchmark/stored_spec_144-182_iter_3_GC_MS.npz

@@ -100,13 +100,13 @@ def test_nit(self):
         nptest.assert_almost_equal(self.orinit, self.optnit, decimal=-2)
 
     def test_intensities(self):
-        nptest.assert_almost_equal(self.oriintensities, self.optintensities, decimal=2)
+        nptest.assert_almost_equal(self.oriintensities, self.optintensities, decimal=4)
 
     def test_widths(self):
-        nptest.assert_almost_equal(self.oriwidths, self.optwidths, decimal=2)
+        nptest.assert_almost_equal(self.oriwidths, self.optwidths, decimal=4)
 
     def test_centers(self):
-        nptest.assert_almost_equal(self.oricenters, self.optcenters, decimal=1)
+        nptest.assert_almost_equal(self.oricenters, self.optcenters, decimal=2)
 
 
 class TestNcp(unittest.TestCase):
@@ -116,15 +116,11 @@ def setUpClass(cls):
         cls.currentResults = AnalysisRunner.get_current_result()
 
     def setUp(self):
-        self.orincp = self.benchmarkResults["all_tot_ncp"][:, :, :-1]
-
+        self.orincp = self.benchmarkResults["all_tot_ncp"]
         self.optncp = self.currentResults.all_tot_ncp
 
     def test_ncp(self):
-        correctNansOri = np.where(
-            (self.orincp == 0) & np.isnan(self.optncp), np.nan, self.orincp
-        )
-        nptest.assert_almost_equal(correctNansOri, self.optncp, decimal=4)
+        nptest.assert_almost_equal(self.orincp, self.optncp, decimal=5)
 
 
 class TestMeanWidths(unittest.TestCase):
@@ -138,7 +134,7 @@ def setUp(self):
         self.optmeanwidths = self.currentResults.all_mean_widths
 
     def test_widths(self):
-        nptest.assert_almost_equal(self.orimeanwidths, self.optmeanwidths, decimal=5)
+        nptest.assert_almost_equal(self.orimeanwidths, self.optmeanwidths, decimal=6)
 
 
 class TestMeanIntensities(unittest.TestCase):