Merge remote-tracking branch 'origin/feature/8373_stretchExpFT'

Code/Mantid/Framework/PythonInterface/plugins/functions/StretchedExpFT.py

@@ -0,0 +1,139 @@
+'''*WIKI* 
+
+Provides the Fourier Transform of the Symmetrized Stretched Exponential Function
+<math> S(Q,E) = Fourier{ height(Q) \cdot e^{-|\frac{x}{tau(Q)}|^{beta(Q)} }</math>
+
+If the energy units of energy are micro-eV, then tau is expressed in pico-seconds. If E-units are micro-eV then
+tau is expressed in nano-seconds.
+*WIKI*
+    
+@author Jose Borreguero, NScD
+@date October 06, 2013
+
+Copyright &copy; 2007-8 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
+
+This file is part of Mantid.
+
+Mantid is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 3 of the License, or
+(at your option) any later version.
+
+Mantid is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+File change history is stored at: <https://github.com/mantidproject/mantid>
+Code Documentation is available at: <http://doxygen.mantidproject.org>
+'''
+
+from mantid.api import IFunction1D, FunctionFactory #, BoundaryConstraint
+from mantid import logger
+import numpy as np
+import copy
+
+from pdb import set_trace as tr
+
+class StretchedExpFT(IFunction1D):
+
+    def __init__(self):
+        '''declare some constants'''
+        super(StretchedExpFT, self).__init__()
+        self._meV2ps = 4.136
+        self._parmset = set(['height','tau','beta']) #valid syntaxfor python >= 2.6
+        self._parm2index = {'height':0,'tau':1,'beta':2} #order in which they were defined
+
+    def category(self):
+        return 'QENS'
+
+    def init(self):
+        '''Declare parameters that participate in the fitting'''
+        # Active fitting parameters
+        self.declareParameter('height', 0.1, 'Intensity at the origin')
+        self.declareParameter('tau', 100.0, 'Relaxation time of the standard exponential')
+        self.declareParameter('beta',1.0, 'Stretching exponent')
+        #self.addConstraints() # constraints not yet exposed to python
+        self.validateParams()
+
+    def validateParams(self):
+        '''Check parameters are positive'''
+        height = self.getParameterValue('height')
+        tau = self.getParameterValue('tau')
+        beta = self.getParameterValue('beta')
+        for name,value in {'height':height, 'tau':tau, 'beta':beta}.items():
+            if value <=0:
+                message = 'Parameter {} in StretchedExpFT must be positive. Got {} instead'.format(name, str(value))
+                logger.error(message)
+                #raise ValueError(message)
+                return None
+        return {'height':height, 'tau':tau, 'beta':beta}
+
+    def function1D(self, xvals, **optparms):
+        '''
+        Computes the Fourier transform of the Symmetrized Stretched Exponential
+        
+        The Symmetrized Stretched Exponential:
+                Fourier{ height * exp( - |t/tau|**beta ) }
+        
+        Given a time step dt and M=2*N+1 time points (N negative, one at zero, N positive),
+        then fft will sample frequencies [0, 1/(M*dt), N/(M*dt), -N/(M*dt), (-N+1)/(M*dt),..,1/(M*dt)].
+        
+        Given xvals, let:
+            1/(M*dt) = xvals[1]-xvals[0]
+            N/(M*dt) = max(abs(xvals))
+        Thus:
+            dt = 1/[M*(xvals[1]-xvals[0])]  # M=2*N+1
+            N = max(abs(xvals)) / (xvals[1]-xvals[0])
+
+        Its Fourier transform is real by definition, thus we return the real part
+        of the Fast Fourier Transform (FFT). The FFT step is meant to produce
+        some residual imaginary part due to numerical inaccuracies which we ignore.
+        
+        We take the absolute value of the real part. The Fourier transform introduces
+        an extra factor exp(i*pi*E/de) which amounts to alternating sign every
+        time E increases by de, the energy bin width
+        '''
+        import scipy.fftpack
+        import scipy.interpolate
+
+        p=self.validateParams()
+        if not p:
+            return np.zeros(len(xvals), dtype=float) # return zeros if parameters not valid
+        # override parameter values with optparms (used for the numerical derivative)
+        if optparms:
+            if self._parmset.issubset( set(optparms.keys()) ):
+                for name in self._parmset: p[name] = optparms[name]
+        de = xvals[1]-xvals[0] # meV (or ueV) , energy step
+        # make sure M > len(xvals) so that we can use interp1d later
+        N = 1+ int( max(np.abs(xvals)) / de )
+        M = 2*N+1
+        dt = self._meV2ps / (M*de) # ps ( or ns), time step
+        sampled_times = dt * np.arange(-N, N+1)
+        exponent = -(np.abs(sampled_times)/p['tau'])**p['beta']
+        freqs = de * np.arange(-N, N+1)
+        fourier = p['height']*np.abs( scipy.fftpack.fft( np.exp(exponent) ).real )
+        fourier = np.concatenate( (fourier[N+1:],fourier[0:N+1]) )
+        interpolator = scipy.interpolate.interp1d(freqs, fourier)
+        fourier = interpolator(xvals)
+        return fourier
+
+    def functionDeriv1D(self, xvals, jacobian):
+        '''Numerical derivative'''
+        p = self.validateParams()
+        f0 = self.function1D(xvals)
+        dp = {}
+        for (key,val) in p.items(): dp[key] = 0.1 * val #modify by ten percent
+        for name in self._parmset:
+            pp = copy.copy(p)
+            pp[name] += dp[name]
+            df = (self.function1D(xvals, **pp) - f0) / dp[name]
+            ip = self._parm2index[name]
+            for ix in range(len(xvals)):
+                jacobian.set(ix, ip, df[ix])  
+
+# Required to have Mantid recognise the new function
+FunctionFactory.subscribe(StretchedExpFT)

Code/Mantid/Framework/PythonInterface/test/python/plugins/functions/CMakeLists.txt

@@ -5,6 +5,7 @@
 set ( TEST_PY_FILES
   Example1DFunctionTest.py
   ExamplePeakFunctionTest.py
+  StretchedExpFTTest.py  # comment until scipy installed on the MAC
 )
 
 check_tests_valid ( ${CMAKE_CURRENT_SOURCE_DIR} ${TEST_PY_FILES} )

Code/Mantid/Framework/PythonInterface/test/python/plugins/functions/StretchedExpFTTest.py

@@ -0,0 +1,125 @@
+import unittest
+from mantid.kernel import *
+from mantid.api import *
+from mantid.simpleapi import Fit
+import testhelpers
+
+from pdb import set_trace as tr
+
+class _InternalMakeSEFTData(PythonAlgorithm):
+
+    def PyInit(self):
+        self.declareProperty('height', 1.0, validator=FloatBoundedValidator(lower=0))
+        self.declareProperty('tau',    1.0, validator=FloatBoundedValidator(lower=0))
+        self.declareProperty('beta',   1.0, validator=FloatBoundedValidator(lower=0))
+        self.declareProperty('nhalfbins', 100)
+        self.declareProperty('de', 0.0004) # energy step in meV
+        self.declareProperty(MatrixWorkspaceProperty('OutputWorkspace', '', direction = Direction.Output))
+
+    def PyExec(self):
+        '''Create the data workspace containing one spectrum resulting from the
+        Fourier transform of an stretched exponential, plus some noise'''
+        import numpy as np
+        from scipy.fftpack import fft
+        from scipy.interpolate import interp1d
+
+        meV2ps = 4.136 # from meV (or ueV) to ps (or ns)
+
+        nhalfbins = self.getProperty('nhalfbins').value
+        de = self.getProperty('de').value #bin width, in meV or ueV
+        nbins = 2*nhalfbins
+        xvals = de*np.arange(-nhalfbins, nhalfbins) + de/2  # energy values
+
+        #Find the time domain
+        N = 1+nhalfbins
+        M = 2*N+1
+        dt = meV2ps / (M*de) # ps ( or ns), time step
+        sampled_times = dt * np.arange(-N, N+1)
+
+        height = self.getProperty('height').value
+        tau = self.getProperty('tau').value
+        beta = self.getProperty('beta').value
+
+        #Define the stretched exponential on the time domain, then do its Fourier transform
+        exponent = -(np.abs(sampled_times)/tau)**beta
+        freqs = de * np.arange(-N, N+1)
+        fourier = height*np.abs( fft( np.exp(exponent) ).real )
+        fourier = np.concatenate( (fourier[N+1:],fourier[0:N+1]) )
+        interpolator = interp1d(freqs, fourier)
+        ynominal = interpolator(xvals)
+
+        wspace = WorkspaceFactory.create('Workspace2D', NVectors=1, XLength=nbins, YLength=nbins)
+        # white noise in the [0.2, 1) range, average is 0.6
+        noise = 0.2 + 0.8*np.random.random(nbins)
+        sign = np.random.random(nbins)
+        sign = np.where(sign>0.5, 1, -1) # random sign
+        # average noise is 6% of the signal local intensity
+
+        yvals = (1+0.1*noise*sign) * ynominal
+        error = 0.2*noise*ynominal
+        wspace.dataX(0)[:] = xvals
+        wspace.dataY(0)[:] = yvals
+        # average error is 12% of the signal local intensity
+        wspace.dataE(0)[:] = error
+
+        self.setProperty('OutputWorkspace', wspace) # Stores the workspace as the given name
+
+class StretchedExpFTTest(unittest.TestCase):
+
+    def skipTest(self):
+        try:
+            import scipy.fftpack # Scipy not available in windows debug
+            return False
+        except ImportError:
+            return True
+
+    def test_registered(self):
+        try:
+            FunctionFactory.createFunction('StretchedExpFT')
+        except RuntimeError, exc:
+            self.fail('Could not create StretchedExpFT function: %s' % str(exc))
+
+    def test_fit(self):
+        if self.skipTest(): # python2.6 doesn't have skipping decorators
+            return
+
+        from random import random
+        variation = lambda x: x*( 1+(random()-0.5)/5. )  # range [x*0.9, x*1.1] should be bigger but not until parameter constraints have been exposed to python 
+        #Generate a data workspace using random parameters around {'height':0.1,'tau':100,'beta':1}
+        parms={'height':variation(0.1), 'tau':variation(100), 'beta':variation(1)}
+
+        Nh = 4000
+        de = 0.0004
+        AlgorithmFactory.subscribe(_InternalMakeSEFTData)
+        alg = testhelpers.run_algorithm('_InternalMakeSEFTData', nhalfbins=Nh, de=de,
+                                         OutputWorkspace='_test_seft_data', **parms)
+        input_ws = alg.getProperty('OutputWorkspace').value
+
+        sx = -Nh*de + de/2
+        ex = (Nh-1)*de + de/2
+        func_string = 'name=StretchedExpFT,height=0.1,tau=100,beta=1'
+        Fit(Function=func_string, InputWorkspace='_test_seft_data',
+            StartX=sx, EndX=ex, CreateOutput=1, MaxIterations=20)
+
+        ws = mtd['_test_seft_data_Parameters']
+        fitted=''
+        for irow in range( ws.rowCount() ):
+            row = ws.row(irow)
+            name = row['Name']
+            if name == 'Cost function value':
+                value = row['Value']
+            elif name in parms.keys():
+                fitted += '{0}={1}, '.format(name, row['Value'])
+
+        target = ', '.join( '{0}={1}'.format(key,val) for key,val in parms.items() )       
+        msg='Cost function {0} too high\nTargets were {1},\nbut obtained {2}'.format(value,target,fitted)
+        self.assertTrue(value < 5.0, msg)
+        msg='Cost function {0}\nStarted with height=0.1, tau=100, beta=1\nTargets were {1},\nobtained {2}'.format(value,target,fitted)
+
+        mtd.remove('_test_seft_data')
+        mtd.remove('_test_seft_data_NormalisedCovarianceMatrix')
+        mtd.remove('_test_seft_data_Parameters')
+        mtd.remove('_test_seft_data_Workspace')
+
+if __name__ == '__main__':
+    unittest.main()