Commit a7fd8abb authored by Shahroz Ahmed's avatar Shahroz Ahmed
Browse files

Merge branch 'master' of github.com:mantidproject/mantid into...

Merge branch 'master' of github.com:mantidproject/mantid into 14923_unit_tests_for_cropped_calib_engg_gui
parents 3a28e436 f6583045
......@@ -114,8 +114,9 @@ private:
/// Calculates the distance of the current solution
double distance(const Kernel::DblMatrix &s2, const std::vector<double> &beta);
/// Populates the output workspaces
void populateOutputWS(const API::MatrixWorkspace_sptr &inWS, size_t spec,
size_t nspec, const std::vector<double> &data,
void populateOutputWS(const API::MatrixWorkspace_sptr &inWS, bool complex,
size_t spec, size_t nspec,
const std::vector<double> &data,
const std::vector<double> &image,
API::MatrixWorkspace_sptr &outData,
API::MatrixWorkspace_sptr &outImage);
......
......@@ -57,6 +57,12 @@ void MaxEnt::init() {
new WorkspaceProperty<>("InputWorkspace", "", Direction::Input),
"An input workspace.");
declareProperty("ComplexData", false,
"Whether or not the input data are complex. If true, the "
"input workspace is expected to have an even number of "
"histograms, with real and imaginary parts arranged in "
"consecutive workspaces");
auto mustBeNonNegative = boost::make_shared<BoundedValidator<double>>();
mustBeNonNegative->setLower(1E-12);
declareProperty(new PropertyWithValue<double>("A", 0.4, mustBeNonNegative,
......@@ -122,6 +128,12 @@ std::map<std::string, std::string> MaxEnt::validateInputs() {
}
}
size_t nhistograms = inWS->getNumberHistograms();
bool complex = getProperty("ComplexData");
if (complex && (nhistograms % 2))
result["InputWorkspace"] = "The number of histograms in the input "
"workspace must be even for complex data";
return result;
}
......@@ -132,6 +144,8 @@ void MaxEnt::exec() {
// Read input workspace
MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");
// Complex data?
bool complex = getProperty("ComplexData");
// Background (default level, sky background, etc)
double background = getProperty("A");
......@@ -154,25 +168,39 @@ void MaxEnt::exec() {
size_t npoints = inWS->blocksize();
// Number of X bins
size_t npointsX = inWS->isHistogramData() ? npoints + 1 : npoints;
// We need to handle complex data
npoints *= 2;
// Create output workspaces
MatrixWorkspace_sptr outImageWS =
WorkspaceFactory::Instance().create(inWS, 2 * nspec, npointsX, npoints);
MatrixWorkspace_sptr outDataWS =
WorkspaceFactory::Instance().create(inWS, nspec, npointsX, npoints);
// Create evol workspaces
MatrixWorkspace_sptr outEvolChi =
WorkspaceFactory::Instance().create(inWS, nspec, niter, niter);
MatrixWorkspace_sptr outEvolTest =
WorkspaceFactory::Instance().create(inWS, nspec, niter, niter);
// Start distribution (flat background)
std::vector<double> image(npoints, background);
// Output workspaces
MatrixWorkspace_sptr outImageWS;
MatrixWorkspace_sptr outDataWS;
MatrixWorkspace_sptr outEvolChi;
MatrixWorkspace_sptr outEvolTest;
if (complex) {
outImageWS =
WorkspaceFactory::Instance().create(inWS, nspec, npointsX, npoints);
outDataWS =
WorkspaceFactory::Instance().create(inWS, nspec, npointsX, npoints);
outEvolChi =
WorkspaceFactory::Instance().create(inWS, nspec / 2, niter, niter);
outEvolTest =
WorkspaceFactory::Instance().create(inWS, nspec / 2, niter, niter);
nspec = nspec / 2;
} else {
outImageWS =
WorkspaceFactory::Instance().create(inWS, 2 * nspec, npointsX, npoints);
outDataWS =
WorkspaceFactory::Instance().create(inWS, nspec, npointsX, npoints);
outEvolChi = WorkspaceFactory::Instance().create(inWS, nspec, niter, niter);
outEvolTest =
WorkspaceFactory::Instance().create(inWS, nspec, niter, niter);
}
// We need to handle complex data
npoints *= 2;
for (size_t s = 0; s < nspec; s++) {
// Start distribution (flat background)
std::vector<double> image(npoints, background);
// Read data from the input workspace
// Only real part, complex part is zero
std::vector<double> data(npoints, 0.);
......@@ -181,6 +209,12 @@ void MaxEnt::exec() {
data[2 * i] = inWS->readY(s)[i];
error[2 * i] = inWS->readE(s)[i];
}
if (complex) {
for (size_t i = 0; i < npoints / 2; i++) {
data[2 * i + 1] = inWS->readY(2 * s + 1)[i];
error[2 * i + 1] = inWS->readE(2 * s + 1)[i];
}
}
// To record the algorithm's progress
std::vector<double> evolChi(niter, 0.);
......@@ -248,7 +282,7 @@ void MaxEnt::exec() {
std::vector<double> solutionData = transformImageToData(newImage);
// Populate the output workspaces
populateOutputWS(inWS, s, nspec, solutionData, newImage, outDataWS,
populateOutputWS(inWS, complex, s, nspec, solutionData, newImage, outDataWS,
outImageWS);
// Populate workspaces recording the evolution of Chi and Test
......@@ -756,8 +790,21 @@ double MaxEnt::distance(const DblMatrix &s2, const std::vector<double> &beta) {
return dist;
}
void MaxEnt::populateOutputWS(const MatrixWorkspace_sptr &inWS, size_t spec,
size_t nspec, const std::vector<double> &data,
/** Populates the output workspaces
* @param inWS :: [input] The input workspace
* @param complex :: [input] Whether or not the input is complex
* @param spec :: [input] The current spectrum being analyzed
* @param nspec :: [input] The total number of histograms in the input workspace
* @param data :: [input] The reconstructed data
* @param image :: [input] The reconstructed image
* @param outData :: [output] The output workspace containing the reconstructed
* data
* @param outImage :: [output] The output workspace containing the reconstructed
* image
*/
void MaxEnt::populateOutputWS(const MatrixWorkspace_sptr &inWS, bool complex,
size_t spec, size_t nspec,
const std::vector<double> &data,
const std::vector<double> &image,
MatrixWorkspace_sptr &outData,
MatrixWorkspace_sptr &outImage) {
......@@ -774,10 +821,18 @@ void MaxEnt::populateOutputWS(const MatrixWorkspace_sptr &inWS, size_t spec,
// Reconstructed data
for (int i = 0; i < npoints; i++)
for (int i = 0; i < npoints; i++) {
YR[i] = data[2 * i];
YI[i] = data[2 * i + 1];
}
outData->dataX(spec) = inWS->readX(spec);
outData->dataY(spec).assign(YR.begin(), YR.end());
outData->dataE(spec).assign(E.begin(), E.end());
if (complex) {
outData->dataX(nspec + spec) = inWS->readX(spec);
outData->dataY(nspec + spec).assign(YI.begin(), YI.end());
outData->dataE(nspec + spec).assign(E.begin(), E.end());
}
// Reconstructed image
......@@ -794,15 +849,14 @@ void MaxEnt::populateOutputWS(const MatrixWorkspace_sptr &inWS, size_t spec,
if (npointsX == npoints + 1)
X[npoints] = X[npoints - 1] + delta;
// Real part
outImage->dataX(spec).assign(X.begin(), X.end());
outImage->dataY(spec).assign(YR.begin(), YR.end());
outImage->dataE(spec).assign(E.begin(), E.end());
// Imaginary part
outImage->dataX(nspec + spec).assign(X.begin(), X.end());
outImage->dataY(nspec + spec).assign(YI.begin(), YI.end());
// No errors
outData->dataE(spec).assign(E.begin(), E.end());
outImage->dataE(spec).assign(E.begin(), E.end());
outImage->dataE(spec + nspec).assign(E.begin(), E.end());
outImage->dataE(nspec + spec).assign(E.begin(), E.end());
}
} // namespace Algorithms
......
......@@ -5,6 +5,7 @@
#include "MantidAPI/FrameworkManager.h"
#include "MantidAPI/AlgorithmManager.h"
#include "MantidTestHelpers/WorkspaceCreationHelper.h"
using namespace Mantid::API;
using Mantid::MantidVec;
......@@ -22,9 +23,96 @@ public:
TS_ASSERT(alg->isInitialized())
}
void test_real_data() {
// Run one iteration, we just want to test the output workspaces' dimensions
int nHist = 5;
int nBins = 10;
auto ws = WorkspaceCreationHelper::Create2DWorkspace(nHist, nBins);
IAlgorithm_sptr alg = AlgorithmManager::Instance().create("MaxEnt");
alg->initialize();
alg->setChild(true);
alg->setProperty("InputWorkspace", ws);
alg->setPropertyValue("MaxIterations", "1");
alg->setPropertyValue("ReconstructedImage", "image");
alg->setPropertyValue("ReconstructedData", "data");
alg->setPropertyValue("EvolChi", "evolChi");
alg->setPropertyValue("EvolAngle", "evolAngle");
TS_ASSERT_THROWS_NOTHING(alg->execute());
MatrixWorkspace_sptr data = alg->getProperty("ReconstructedData");
MatrixWorkspace_sptr image = alg->getProperty("ReconstructedImage");
MatrixWorkspace_sptr chi = alg->getProperty("EvolChi");
MatrixWorkspace_sptr angle = alg->getProperty("EvolAngle");
TS_ASSERT_EQUALS(data->getNumberHistograms(), nHist);
TS_ASSERT_EQUALS(image->getNumberHistograms(), nHist * 2);
TS_ASSERT_EQUALS(chi->getNumberHistograms(), nHist);
TS_ASSERT_EQUALS(angle->getNumberHistograms(), nHist);
TS_ASSERT_EQUALS(data->blocksize(), nBins);
TS_ASSERT_EQUALS(image->blocksize(), nBins);
TS_ASSERT_EQUALS(chi->blocksize(), 1);
TS_ASSERT_EQUALS(angle->blocksize(), 1);
}
void test_complex_data() {
// Run one iteration, we just want to test the output workspaces' dimensions
int nHist = 6;
int nBins = 10;
auto ws = WorkspaceCreationHelper::Create2DWorkspace(nHist, nBins);
IAlgorithm_sptr alg = AlgorithmManager::Instance().create("MaxEnt");
alg->initialize();
alg->setChild(true);
alg->setProperty("InputWorkspace", ws);
alg->setProperty("ComplexData", true);
alg->setPropertyValue("MaxIterations", "1");
alg->setPropertyValue("ReconstructedImage", "image");
alg->setPropertyValue("ReconstructedData", "data");
alg->setPropertyValue("EvolChi", "evolChi");
alg->setPropertyValue("EvolAngle", "evolAngle");
TS_ASSERT_THROWS_NOTHING(alg->execute());
MatrixWorkspace_sptr data = alg->getProperty("ReconstructedData");
MatrixWorkspace_sptr image = alg->getProperty("ReconstructedImage");
MatrixWorkspace_sptr chi = alg->getProperty("EvolChi");
MatrixWorkspace_sptr angle = alg->getProperty("EvolAngle");
TS_ASSERT_EQUALS(data->getNumberHistograms(), nHist);
TS_ASSERT_EQUALS(image->getNumberHistograms(), nHist);
TS_ASSERT_EQUALS(chi->getNumberHistograms(), nHist / 2);
TS_ASSERT_EQUALS(angle->getNumberHistograms(), nHist / 2);
TS_ASSERT_EQUALS(data->blocksize(), nBins);
TS_ASSERT_EQUALS(image->blocksize(), nBins);
TS_ASSERT_EQUALS(chi->blocksize(), 1);
TS_ASSERT_EQUALS(angle->blocksize(), 1);
}
void test_bad_complex_data() {
auto ws = WorkspaceCreationHelper::Create2DWorkspace(5, 10);
IAlgorithm_sptr alg = AlgorithmManager::Instance().create("MaxEnt");
alg->initialize();
alg->setChild(true);
alg->setProperty("InputWorkspace", ws);
alg->setProperty("ComplexData", true);
alg->setPropertyValue("MaxIterations", "1");
alg->setPropertyValue("ReconstructedImage", "image");
alg->setPropertyValue("ReconstructedData", "data");
alg->setPropertyValue("EvolChi", "evolChi");
alg->setPropertyValue("EvolAngle", "evolAngle");
TS_ASSERT_THROWS_ANYTHING(alg->execute());
}
void test_cosine() {
auto ws = createWorkspace(0.0);
auto ws = createWorkspace(50, 0.0);
IAlgorithm_sptr alg = AlgorithmManager::Instance().create("MaxEnt");
alg->initialize();
......@@ -57,7 +145,7 @@ public:
void test_sine() {
auto ws = createWorkspace(M_PI / 2.);
auto ws = createWorkspace(50, M_PI / 2.);
IAlgorithm_sptr alg = AlgorithmManager::Instance().create("MaxEnt");
alg->initialize();
......@@ -88,14 +176,12 @@ public:
TS_ASSERT_DELTA(data->readY(0)[27], 0.7218, 0.0001);
}
MatrixWorkspace_sptr createWorkspace(double phase) {
MatrixWorkspace_sptr createWorkspace(size_t maxt, double phase) {
// Create cosine with phase 'phase'
// Frequency of the oscillations
double w = 1.6;
// Maximum time
size_t maxt = 50;
MantidVec X;
MantidVec Y;
......@@ -106,7 +192,6 @@ public:
Y.push_back(cos(w * x + phase));
E.push_back(0.1);
}
auto createWS = AlgorithmManager::Instance().create("CreateWorkspace");
createWS->initialize();
createWS->setChild(true);
......
......@@ -20,7 +20,7 @@ default.instrument =
Q.convention = Inelastic
# Set of PyQt interfaces to add to the Interfaces menu, separated by a space. Interfaces are seperated from their respective categories by a "/".
mantidqt.python_interfaces = Direct/DGS_Reduction.py Direct/DGSPlanner.py SANS/ORNL_SANS.py Reflectometry/REFL_Reduction.py Reflectometry/REFL_SF_Calculator.py Reflectometry/REFM_Reduction.py Utility/TofConverter.py Reflectometry/ISIS_Reflectometry.py Diffraction/Powder_Diffraction_Reduction.py Utility/FilterEvents.py Diffraction/HFIR_Powder_Diffraction_Reduction.py Diffraction/HFIR_4Circle_Reduction.py
mantidqt.python_interfaces = Direct/DGS_Reduction.py Direct/DGSPlanner.py SANS/ORNL_SANS.py Reflectometry/REFL_Reduction.py Reflectometry/REFL_SF_Calculator.py Reflectometry/REFM_Reduction.py Utility/TofConverter.py Reflectometry/ISIS_Reflectometry.py Diffraction/Powder_Diffraction_Reduction.py Utility/FilterEvents.py Diffraction/HFIR_Powder_Diffraction_Reduction.py Diffraction/HFIR_4Circle_Reduction.py Utility/QECoverage.py
mantidqt.python_interfaces_directory = @MANTID_ROOT@/scripts
......
#pylint: disable=invalid-name
# pylint: disable=invalid-name
'''
@author Jose Borreguero, NScD
@date October 06, 2013
......@@ -24,49 +24,49 @@ 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.api import IFunction1D, FunctionFactory
from mantid import logger
import numpy as np
import copy
class StretchedExpFT(IFunction1D):
#pylint: disable=super-on-old-class
# pylint: disable=super-on-old-class
def __init__(self):
'''declare some constants'''
"""declare some constants"""
super(StretchedExpFT, self).__init__()
self._h = 4.135665616 #meV*Thz
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
self._h = 4.135665616 # meV*Thz
self._parmset = set(['height', 'tau', 'beta', 'Origin']) # valid syntax for python >= 2.6
self._parm2index = {'height': 0, 'tau': 1, 'beta': 2, 'Origin': 3} # order in which they were defined
def category(self):
return 'QuasiElastic'
def init(self):
'''Declare parameters that participate in the fitting'''
"""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.declareParameter('beta', 1.0, 'Stretching exponent')
self.declareParameter('Origin', 0.0, 'Origin of the peak')
self.validateParams()
def validateParams(self):
'''Check parameters are positive'''
"""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:
origin = self.getParameterValue('Origin')
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)
# raise ValueError(message)
return None
return {'height':height, 'tau':tau, 'beta':beta}
return {'height': height, 'tau': tau, 'beta': beta, 'Origin': origin}
def function1D(self, xvals, **optparms):
'''
Computes the Fourier transform of the Symmetrized Stretched Exponential
""" Computes the Fourier transform of the Symmetrized Stretched Exponential
The Symmetrized Stretched Exponential:
Fourier{ height * exp( - |t/tau|**beta ) }
......@@ -88,42 +88,44 @@ class StretchedExpFT(IFunction1D):
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()
p = self.validateParams()
if not p:
return np.zeros(len(xvals), dtype=float) # return zeros if parameters not valid
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()) ):
if self._parmset.issubset(set(optparms.keys())):
for name in self._parmset:
p[name] = optparms[name]
de = (xvals[1]-xvals[0]) / 2 # increase the long-time range, increase the low-frequency resolution
emax = max( abs(xvals) )
Nmax = 16000 # increase short-times resolution, increase the resolution of the structure factor tail
while ( emax / de ) < Nmax:
emax = 2 * emax
N = int( emax / de )
dt = ( float(N) / ( 2 * N + 1 ) ) * (self._h / emax) # extent to negative times and t==0
sampled_times = dt * np.arange(-N, N+1) # len( sampled_times ) < 64000
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]) )
de = (xvals[1] - xvals[0]) / 2 # increase the long-time range, increase the low-frequency resolution
emax = max(abs(xvals))
nmax = 16000 # increase short-times resolution, increase the resolution of the structure factor tail
while (emax / de) < nmax:
emax *= 2
N = int(emax / de)
dt = (float(N) / (2 * N + 1)) * (self._h / emax) # extent to negative times and t==0
sampled_times = dt * np.arange(-N, N + 1) # len( sampled_times ) < 64000
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)
fourier = interpolator(xvals - p['Origin'])
return fourier
def functionDeriv1D(self, xvals, jacobian):
'''Numerical derivative'''
"""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 (key, val) in p.items():
dp[key] = 0.1 * val # modify by ten percent
if val == 0.0:
dp[key] = 1E-06 # some small, supposedly, increment
for name in self._parmset:
pp = copy.copy(p)
pp[name] += dp[name]
......@@ -132,5 +134,5 @@ class StretchedExpFT(IFunction1D):
for ix in range(len(xvals)):
jacobian.set(ix, ip, df[ix])
# Required to have Mantid recognise the new function
FunctionFactory.subscribe(StretchedExpFT)
FunctionFactory.subscribe(StretchedExpFT) # Required to have Mantid recognise the new function
......@@ -3,126 +3,140 @@ from mantid.kernel import *
from mantid.api import *
from mantid.simpleapi import Fit
import testhelpers
#from pdb import set_trace as tr
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
plank_constant = 4.135665616 # in units of meV*THz
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 + 1) + de/2 #
#Find the time domain
M = 2 * nhalfbins + 1
dt = plank_constant / ( M * de ) # ps ( or ns), time step
sampled_times = dt * np.arange( -nhalfbins, nhalfbins + 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( -nhalfbins, nhalfbins + 1 )
fourier = height * np.abs( fft( np.exp( exponent ) ).real )
fourier = np.concatenate( ( fourier[ nhalfbins + 1 : ],fourier[ 0 : nhalfbins + 1 ] ) )
interpolator = interp1d( freqs, fourier )
N = nhalfbins / 2
ynominal = interpolator( xvals[ N : -N-1 ] ) # high energies are not well represented, avoid them
wspace = WorkspaceFactory.create( 'Workspace2D' , NVectors = 1, XLength = 2 * N, YLength = 2 * N )
# white noise in the [0.2, 1) range, average is 0.6
noise = 0.2 + 0.8 * np.random.random( 2 * N )
sign = np.random.random( 2 * N )
sign = np.where( sign > 0.5, 1, -1 ) # random sign
# average noise is 6% of the signal local intensity
yvals = ( 1 + 0.05 * noise * sign ) * ynominal
error = 0.1 * noise * ynominal
wspace.dataX(0)[:] = xvals[ N : -N-1 ]
wspace.dataY(0)[:] = yvals
# average error is 12% of the signal local intensity
wspace.dataE(0)[:] = error
""" # Only for debugging purposes
buffer = '#h = {0}, tau = {1}, beta = {2}\n'.format( height, tau, beta )
for iy in range( len(yvals) ):
buffer += '{0} {1} {2}\n'.format( xvals[N+iy], yvals[iy], error[iy] )
open( '/tmp/junk.dat', 'w' ).write( buffer )
"""
self.setProperty( 'OutputWorkspace', wspace ) # Stores the workspace as the given name
def PyInit(self):
self.declareProperty('height', 1.0, validator=FloatBoundedValidator(lower=0))
self.declareProperty('tau', 100.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
plank_constant = 4.135665616 # in units of meV*THz
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+1)+de/2
# Find the time domain
M = 2*nhalfbins+1
dt = plank_constant/(M*de) # ps ( or ns), time step
sampled_times = dt*np.arange(-nhalfbins, nhalfbins+1)
height = self.getProperty('height').value
tau = self.getProperty('tau').value