Refs #8373 Finished writing the test

modified:   StretchedExpFT.py
modified:   StretchedExpFTTest.py

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

@@ -28,7 +28,7 @@
 from mantid.api import IFunction1D, FunctionFactory #, BoundaryConstraint
 from mantid import logger
 import numpy as np
-from scipy.fftpack import (fft, fftfreq, fftshift)
+from scipy.fftpack import fft
 from scipy.interpolate import interp1d
 import copy
 
@@ -75,7 +75,7 @@ def function1D(self, xvals, **optparms):
                 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 fftfreq will sample frequencies [0, 1/(M*dt), N/(M*dt), -N/(M*dt), (-N+1)/(M*dt),..,1/(M*dt)].
+        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]

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

@@ -4,6 +4,7 @@
 from mantid.simpleapi import Fit
 import testhelpers
 
+from pdb import set_trace as tr
 
 class _InternalMakeSEFTData(PythonAlgorithm):
 
@@ -12,35 +13,55 @@ def PyInit(self):
         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
-        nbins = 1 + 2*nhalfbins
-        wspace = WorkspaceFactory.create('Workspace2D', NVectors=1, XLength=nbins, YLength=nbins)
+        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
-
-        xvals = np.arange(-nhalfbins, 1+nhalfbins) * (1./nhalfbins) # from -1 to 1
-        exponent = -(np.abs(xvals)/tau)**beta
-        ynominal = height*fft( np.exp(exponent) ).real
+        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
+        #tr()
         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)[:] = 0.2*noise*ynominal
-
+        wspace.dataE(0)[:] = error
+        
         self.setProperty('OutputWorkspace', wspace) # Stores the workspace as the given name
 
 class StretchedExpFTTest(unittest.TestCase):
@@ -53,15 +74,21 @@ def test_registered(self):
 
     def test_fit(self):
         from random import random
-        parms={'height':10*random(), 'tau':0.1+random(), 'beta':0.1+2*random()}
+        variation = lambda x: x+x*(random()-0.5)  # range [-x/2, x/2]
+        #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=4000,
+        alg = testhelpers.run_algorithm('_InternalMakeSEFTData', nhalfbins=Nh, de=de,
                                          OutputWorkspace='_test_seft_data', **parms)
         input_ws = alg.getProperty('OutputWorkspace').value
 
-        func_string = 'name=StretechedExpFT,' + ','.join( '{}={}'.format(key,val) for key,val in parms.items() )
-        Fit(Function=func_string, InputWorkspace='_test_data',
-            StartX=-1, EndX=1, CreateOutput=1, MaxIterations=20)
+        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=''
@@ -71,12 +98,13 @@ def test_fit(self):
             if name == 'Cost function value':
                 value = row['Value']
             elif name in parms.keys():
-                fitted += '{}={}'.format(name, row['Value'])
+                fitted += '{}={}, '.format(name, row['Value'])
 
-        target = ','.join( '{}={}'.format(key,val) for key,val in parms.items() )       
-        msg='Cost function {} too high. Targets were {} but obtained {}'.format(value,target,fitted)
+        target = ', '.join( '{}={}'.format(key,val) for key,val in parms.items() )       
+        msg='Cost function {} too high\nTargets were {},\nbut obtained {}'.format(value,target,fitted)
         self.assertLess(value, 5, msg)
-
+        msg='Cost function {}\nStarted with height=0.1, tau=100, beta=1\nTargets were {},\nobtained {}'.format(value,target,fitted)
+
         mtd.remove('_test_seft_data')
         mtd.remove('_test_seft_data_NormalisedCovarianceMatrix')
         mtd.remove('_test_seft_data_Parameters')