Merge pull request #23777 from mantidproject/23776_Indirect_ReplaceOldMsdYiFile

Indirect - Replace old MsdYi.py file with updated version

Framework/PythonInterface/mantid/fitfunctions.py

@@ -23,7 +23,7 @@ class FunctionWrapper(object):
                 return wrapper(fun, *args, **kwargs)
         return FunctionWrapper(fun, **kwargs)
 
-    def __init__ (self, name, **kwargs):
+    def __init__(self, name, **kwargs):
         """
         Called when creating an instance
 

Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/MSDFit.py

@@ -31,7 +31,8 @@ class MSDFit(DataProcessorAlgorithm):
         self.declareProperty(MatrixWorkspaceProperty('InputWorkspace', '', direction=Direction.Input),
                              doc='Sample input workspace')
         self.declareProperty(name='Model', defaultValue='Gauss',
-                             validator=StringListValidator(['Gauss', 'Peters', 'Yi']),
+                             validator=StringListValidator(
+                                 ['Gauss', 'Peters', 'Yi']),
                              doc='Model options : Gauss, Peters, Yi')
 
         self.declareProperty(name='XStart', defaultValue=0.0,
@@ -100,9 +101,12 @@ class MSDFit(DataProcessorAlgorithm):
         progress = Progress(self, 0.0, 0.05, 3)
         self._original_ws = self._input_ws
 
-        rename_alg = self.createChildAlgorithm("RenameWorkspace", enableLogging=False)
+        rename_alg = self.createChildAlgorithm(
+            "RenameWorkspace", enableLogging=False)
         rename_alg.setProperty("InputWorkspace", self._input_ws)
-        rename_alg.setProperty("OutputWorkspace", self._input_ws + "_" + self._model)
+        rename_alg.setProperty(
+            "OutputWorkspace",
+            self._input_ws + "_" + self._model)
         rename_alg.execute()
         self._input_ws = self._input_ws + "_" + self._model
         input_params = [self._input_ws + ',i%d' % i for i in range(self._spec_range[0],
@@ -111,19 +115,19 @@ class MSDFit(DataProcessorAlgorithm):
         # Fit line to each of the spectra
         if self._model == 'Gauss':
             logger.information('Model : Gaussian approximation')
-            function = 'name=MsdGauss, Height=1.0, MSD=0.1'
-            function += ',constraint=(Height>0.0, MSD>0.0)'
-            params_list = ['Height', 'MSD']
+            function = 'name=MsdGauss, Height=1.0, Msd=0.1'
+            function += ',constraint=(Height>0.0, Msd>0.0)'
+            params_list = ['Height', 'Msd']
         elif self._model == 'Peters':
             logger.information('Model : Peters & Kneller')
-            function = 'name=MsdPeters, Height=1.0, MSD=1.0, Beta=1.0'
-            function += ',constraint=(Height>0.0, MSD>0.0, 100.0>Beta>0.3)'
-            params_list = ['Height', 'MSD', 'Beta']
+            function = 'name=MsdPeters, Height=1.0, Msd=1.0, Beta=1.0'
+            function += ',constraint=(Height>0.0, Msd>0.0, 100.0>Beta>0.3)'
+            params_list = ['Height', 'Msd', 'Beta']
         elif self._model == 'Yi':
             logger.information('Model : Yi et al')
-            function = 'name=MsdYi, Height=1.0, MSD=1.0, Sigma=0.1'
-            function += ',constraint=(Height>0.0, MSD>0.0, Sigma>0.0)'
-            params_list = ['Height', 'MSD', 'Sigma']
+            function = 'name=MsdYi, Height=1.0, Msd=1.0, Sigma=0.1'
+            function += ',constraint=(Height>0.0, Msd>0.0, Sigma>0.0)'
+            params_list = ['Height', 'Msd', 'Sigma']
         else:
             raise ValueError('No Model defined')
 
@@ -137,12 +141,20 @@ class MSDFit(DataProcessorAlgorithm):
                            FitType='Sequential',
                            CreateOutput=True)
 
-        delete_alg = self.createChildAlgorithm("DeleteWorkspace", enableLogging=False)
-        delete_alg.setProperty("Workspace", self._output_msd_ws + '_NormalisedCovarianceMatrices')
+        delete_alg = self.createChildAlgorithm(
+            "DeleteWorkspace", enableLogging=False)
+        delete_alg.setProperty(
+            "Workspace",
+            self._output_msd_ws +
+            '_NormalisedCovarianceMatrices')
         delete_alg.execute()
-        delete_alg.setProperty("Workspace", self._output_msd_ws + '_Parameters')
+        delete_alg.setProperty(
+            "Workspace",
+            self._output_msd_ws +
+            '_Parameters')
         delete_alg.execute()
-        rename_alg = self.createChildAlgorithm("RenameWorkspace", enableLogging=False)
+        rename_alg = self.createChildAlgorithm(
+            "RenameWorkspace", enableLogging=False)
         rename_alg.setProperty("InputWorkspace", self._output_msd_ws)
         rename_alg.setProperty("OutputWorkspace", self._output_param_ws)
         rename_alg.execute()
@@ -154,31 +166,50 @@ class MSDFit(DataProcessorAlgorithm):
         for par in params_list:
             ws_name = self._output_msd_ws + '_' + par
             parameter_ws_group.append(ws_name)
-            convert_alg = self.createChildAlgorithm("ConvertTableToMatrixWorkspace", enableLogging=False)
+            convert_alg = self.createChildAlgorithm(
+                "ConvertTableToMatrixWorkspace", enableLogging=False)
             convert_alg.setProperty("InputWorkspace", self._output_param_ws)
             convert_alg.setProperty("OutputWorkspace", ws_name)
             convert_alg.setProperty("ColumnX", 'axis-1')
             convert_alg.setProperty("ColumnY", par)
             convert_alg.setProperty("ColumnE", par + '_Err')
             convert_alg.execute()
-            mtd.addOrReplace(ws_name, convert_alg.getProperty("OutputWorkspace").value)
-
-        append_alg = self.createChildAlgorithm("AppendSpectra", enableLogging=False)
-        append_alg.setProperty("InputWorkspace1", self._output_msd_ws + '_' + params_list[0])
-        append_alg.setProperty("InputWorkspace2", self._output_msd_ws + '_' + params_list[1])
+            mtd.addOrReplace(
+                ws_name, convert_alg.getProperty("OutputWorkspace").value)
+
+        append_alg = self.createChildAlgorithm(
+            "AppendSpectra", enableLogging=False)
+        append_alg.setProperty(
+            "InputWorkspace1",
+            self._output_msd_ws +
+            '_' +
+            params_list[0])
+        append_alg.setProperty(
+            "InputWorkspace2",
+            self._output_msd_ws +
+            '_' +
+            params_list[1])
         append_alg.setProperty("ValidateInputs", False)
         append_alg.setProperty("OutputWorkspace", self._output_msd_ws)
         append_alg.execute()
-        mtd.addOrReplace(self._output_msd_ws, append_alg.getProperty("OutputWorkspace").value)
+        mtd.addOrReplace(
+            self._output_msd_ws,
+            append_alg.getProperty("OutputWorkspace").value)
         if len(params_list) > 2:
             append_alg.setProperty("InputWorkspace1", self._output_msd_ws)
-            append_alg.setProperty("InputWorkspace2", self._output_msd_ws + '_' + params_list[2])
+            append_alg.setProperty(
+                "InputWorkspace2",
+                self._output_msd_ws +
+                '_' +
+                params_list[2])
             append_alg.setProperty("ValidateInputs", False)
             append_alg.setProperty("OutputWorkspace", self._output_msd_ws)
             append_alg.execute()
-            mtd.addOrReplace(self._output_msd_ws, append_alg.getProperty("OutputWorkspace").value)
+            mtd.addOrReplace(self._output_msd_ws,
+                             append_alg.getProperty("OutputWorkspace").value)
         for par in params_list:
-            delete_alg.setProperty("Workspace", self._output_msd_ws + '_' + par)
+            delete_alg.setProperty(
+                "Workspace", self._output_msd_ws + '_' + par)
             delete_alg.execute()
 
         progress.report('Change axes')
@@ -187,7 +218,9 @@ class MSDFit(DataProcessorAlgorithm):
         sort_alg.setProperty("InputWorkspace", self._output_msd_ws)
         sort_alg.setProperty("OutputWorkspace", self._output_msd_ws)
         sort_alg.execute()
-        mtd.addOrReplace(self._output_msd_ws, sort_alg.getProperty("OutputWorkspace").value)
+        mtd.addOrReplace(
+            self._output_msd_ws,
+            sort_alg.getProperty("OutputWorkspace").value)
         # Create a new x axis for the Q and Q**2 workspaces
         xunit = mtd[self._output_msd_ws].getAxis(0).setUnit('Label')
         xunit.setLabel('Temperature', 'K')
@@ -200,7 +233,10 @@ class MSDFit(DataProcessorAlgorithm):
         # Rename fit workspace group
         original_fit_ws_name = self._output_msd_ws + '_Workspaces'
         if original_fit_ws_name != self._output_fit_ws:
-            rename_alg.setProperty("InputWorkspace", self._output_msd_ws + '_Workspaces')
+            rename_alg.setProperty(
+                "InputWorkspace",
+                self._output_msd_ws +
+                '_Workspaces')
             rename_alg.setProperty("OutputWorkspace", self._output_fit_ws)
             rename_alg.execute()
 
@@ -213,7 +249,8 @@ class MSDFit(DataProcessorAlgorithm):
         copy_alg.setProperty("OutputWorkspace", self._output_fit_ws)
         copy_alg.execute()
 
-        rename_alg = self.createChildAlgorithm("RenameWorkspace", enableLogging=False)
+        rename_alg = self.createChildAlgorithm(
+            "RenameWorkspace", enableLogging=False)
         rename_alg.setProperty("InputWorkspace", self._input_ws)
         rename_alg.setProperty("OutputWorkspace", self._original_ws)
         rename_alg.execute()

Framework/PythonInterface/plugins/functions/MsdGauss.py

@@ -27,11 +27,11 @@ class MsdGauss(IFunction1D):
     def init(self):
         # Active fitting parameters
         self.declareParameter("Height", 1.0, 'Height')
-        self.declareParameter("MSD", 0.05, 'Mean square displacement')
+        self.declareParameter("Msd", 0.05, 'Mean square displacement')
 
     def function1D(self, xvals):
         height = self.getParameterValue("Height")
-        msd = self.getParameterValue("MSD")
+        msd = self.getParameterValue("Msd")
 
         xvals = np.array(xvals)
         intensity = height * np.exp(-msd * xvals**2)
@@ -40,7 +40,7 @@ class MsdGauss(IFunction1D):
 
     def functionDeriv1D(self, xvals, jacobian):
         height = self.getParameterValue("Height")
-        msd = self.getParameterValue("MSD")
+        msd = self.getParameterValue("Msd")
 
         for i, x in enumerate(xvals):
             e = math.exp(-msd * x**2)

Framework/PythonInterface/plugins/functions/MsdPeters.py

@@ -28,12 +28,12 @@ class MsdPeters(IFunction1D):
     def init(self):
         # Active fitting parameters
         self.declareParameter("Height", 1.0, 'Height')
-        self.declareParameter("MSD", 0.05, 'Mean square displacement')
+        self.declareParameter("Msd", 0.05, 'Mean square displacement')
         self.declareParameter("Beta", 1.0, 'beta')
 
     def function1D(self, xvals):
         height = self.getParameterValue("Height")
-        msd = self.getParameterValue("MSD")
+        msd = self.getParameterValue("Msd")
         beta = self.getParameterValue("Beta")
 
         xvals = np.array(xvals)
@@ -44,7 +44,7 @@ class MsdPeters(IFunction1D):
 
     def functionDeriv1D(self, xvals, jacobian):
         height = self.getParameterValue("Height")
-        msd = self.getParameterValue("MSD")
+        msd = self.getParameterValue("Msd")
         beta = self.getParameterValue("Beta")
 
         for i, x in enumerate(xvals):

Framework/PythonInterface/plugins/functions/MsdYi.py

-# Mantid Repository : https://github.com/mantidproject/mantid
-#
-# Copyright © 2007 ISIS Rutherford Appleton Laboratory UKRI,
-#     NScD Oak Ridge National Laboratory, European Spallation Source
-#     & Institut Laue - Langevin
-# SPDX - License - Identifier: GPL - 3.0 +
-# pylint: disable=no-init,invalid-name
 '''
 @author Spencer Howells, ISIS
 @date December 05, 2013
+
+Copyright © 2007-8 ISIS Rutherford Appleton Laboratory,
+NScD Oak Ridge National Laboratory & European Spallation Source
+
+pylint: disable=no-init,invalid-name
+
+This file is part of Mantid.
+Mantid Repository : https://github.com/mantidproject/mantid
+File change history is stored at: <https://github.com/mantidproject/mantid>
+Code Documentation is available at: <http://doxygen.mantidproject.org>
 '''
 from __future__ import (absolute_import, division, print_function)
-
 import math
 import numpy as np
-
 from mantid.api import IFunction1D, FunctionFactory
 
-
 # The model of Yi et al(J Phys Chem B 1316 5029 2012) takes into account motional heterogeneity.
-# The elastic intensity is propotional to exp(-1/6*Q^2*msd)*(1+Q^4*sigma/72)
+# The elastic intensity is proportional to exp(-(1/6)*Q^2*msd)*(1+Q^4*sigma/72)
 # where the mean square displacement msd = <r^2> and sigma^2 is the variance of the msd.
+# In the limit sigma = 0 the model becomes Gaussian.
+
 
 class MsdYi(IFunction1D):
+
     def category(self):
         return "QuasiElastic"
 
     def init(self):
         # Active fitting parameters
         self.declareParameter("Height", 1.0, 'Height')
-        self.declareParameter("MSD", 0.05, 'Mean square displacement')
+        self.declareParameter("Msd", 0.05, 'Mean square displacement')
         self.declareParameter("Sigma", 1.0, 'Sigma')
 
     def function1D(self, xvals):
         height = self.getParameterValue("Height")
-        msd = self.getParameterValue("MSD")
+        msd = self.getParameterValue("Msd")
         sigma = self.getParameterValue("Sigma")
 
         xvals = np.array(xvals)
-        i1 = np.exp(-1.0 / (6.0 * xvals**2 * msd))
+        i1 = np.exp((-1.0 / 6.0) * xvals * xvals * msd)
         i2 = 1.0 + (np.power(xvals, 4) * sigma / 72.0)
         intensity = height * i1 * i2
 
@@ -45,13 +48,12 @@ class MsdYi(IFunction1D):
 
     def functionDeriv1D(self, xvals, jacobian):
         height = self.getParameterValue("Height")
-        msd = self.getParameterValue("MSD")
+        msd = self.getParameterValue("Msd")
         sigma = self.getParameterValue("Sigma")
 
         for i, x in enumerate(xvals):
-            q = msd * x**2
-            f1 = math.exp(-1.0 / (6.0 * q))
-            df1 = f1 / (6.0 * x * q)
+            f1 = math.exp((-1.0 / 6.0) * x * x * msd)
+            df1 = -f1 * ((2.0 / 6.0) * x * msd)
             x4 = math.pow(x, 4)
             f2 = 1.0 + (x4 * sigma / 72.0)
             df2 = x4 / 72.0

Framework/PythonInterface/test/python/plugins/functions/MsdGaussTest.py

@@ -23,15 +23,15 @@ class MsdGaussTest(unittest.TestCase):
         input = np.array([[0, 1], [2, 3]])
         expected = np.array([[1., 0.95122942], [0.81873075, 0.63762815]])
         tolerance = 0.000001
-        status, output = check_output("MsdGauss", input, expected, tolerance, Height=1.0, MSD=0.05)
+        status, output = check_output("MsdGauss", input, expected, tolerance, Height=1.0, Msd=0.05)
 
         if not status:
             self.fail("Computed output " + str(output) + " from input " + str(input) +
                       " is not equal to the expected output: " + str(expected))
 
     def test_use_in_fit(self):
-        workspace = create_test_workspace(create_model("MsdGauss", Height=1.0, MSD=0.05), 1000)
-        function_string = create_function_string("MsdGauss", Height=1.0, MSD=0.05)
+        workspace = create_test_workspace(create_model("MsdGauss", Height=1.0, Msd=0.05), 1000)
+        function_string = create_function_string("MsdGauss", Height=1.0, Msd=0.05)
         Fit(Function=function_string, InputWorkspace=workspace, StartX=1.2, EndX=1200)
 
 

Framework/PythonInterface/test/python/plugins/functions/MsdPetersTest.py

@@ -23,15 +23,15 @@ class MsdPetersTest(unittest.TestCase):
         input = np.array([[0, 1], [2, 3]])
         expected = np.array([[1., 0.99173554], [0.96774194, 0.93023256]])
         tolerance = 0.000001
-        status, output = check_output("MsdPeters", input, expected, tolerance, Height=1.0, MSD=0.05, Beta=1.0)
+        status, output = check_output("MsdPeters", input, expected, tolerance, Height=1.0, Msd=0.05, Beta=1.0)
 
         if not status:
             self.fail("Computed output " + str(output) + " from input " + str(input) +
                       " is not equal to the expected output: " + str(expected))
 
     def test_use_in_fit(self):
-        workspace = create_test_workspace(create_model("MsdPeters", Height=1.0, MSD=0.05, Beta=1.0), 1000)
-        function_string = create_function_string("MsdPeters", Height=1.0, MSD=0.05, Beta=1.0)
+        workspace = create_test_workspace(create_model("MsdPeters", Height=1.0, Msd=0.05, Beta=1.0), 1000)
+        function_string = create_function_string("MsdPeters", Height=1.0, Msd=0.05, Beta=1.0)
         Fit(Function=function_string, InputWorkspace=workspace, StartX=1.2, EndX=1200)
 
 

Framework/PythonInterface/test/python/plugins/functions/MsdYiTest.py

@@ -21,18 +21,30 @@ class MsdYiTest(unittest.TestCase):
 
     def test_function_output(self):
         input = np.array([[1, 2], [3, 4]])
-        expected = np.array([[0.03616947, 0.53117559], [1.46726692, 3.69882113]])
+        expected = np.array(
+            [[1.00547492, 1.18215301], [1.97145491, 3.98690068]])
         tolerance = 0.000001
-        status, output = check_output("MsdYi", input, expected, tolerance, Height=1.0, MSD=0.05, Sigma=1.0)
+        status, output = check_output(
+            "MsdYi", input, expected, tolerance, Height=1.0, Msd=0.05, Sigma=1.0)
 
         if not status:
             self.fail("Computed output " + str(output) + " from input " + str(input) +
                       " is not equal to the expected output: " + str(expected))
 
     def test_use_in_fit(self):
-        workspace = create_test_workspace(create_model("MsdYi", Height=1.0, MSD=0.05, Sigma=1.0), 1000)
-        function_string = create_function_string("MsdYi", Height=1.0, MSD=0.05, Sigma=1.0)
-        Fit(Function=function_string, InputWorkspace=workspace, StartX=1.2, EndX=1200)
+        workspace = create_test_workspace(
+            create_model(
+                "MsdYi",
+                Height=1.0,
+                Msd=0.05,
+                Sigma=1.0),
+            1000)
+        function_string = create_function_string(
+            "MsdYi", Height=1.0, Msd=0.05, Sigma=1.0)
+        Fit(Function=function_string,
+            InputWorkspace=workspace,
+            StartX=1.2,
+            EndX=1200)
 
 
 if __name__ == '__main__':

docs/source/algorithms/MSDFit-v1.rst

@@ -61,8 +61,8 @@ Output:
     A0: [ 0.87079958]
     A1: [ 0.03278263]
     Using Yi Model
-    A0: [ 0.75677983]
-    A1: [ 1.76943372]
+    A0: [ 0.95770532]
+    A1: [ 0.58819225]
 
 .. categories::
 

docs/source/release/v3.14.0/indirect_inelastic.rst

@@ -49,6 +49,7 @@ Bugfixes
   message is now displayed.
 - The Probability Density Functions (PDF) workspaces for the FABADA minimiser in ConvFit no longer overwrite each other. 
   Various other improvements in the display of the FABADA PDF's have also been finished.
+- The expression for the Fit type Yi in MSDFit was incorrect and has now been corrected.
 
 
 Data Corrections Interface