forked from AmbaPant/mantid
-
Notifications
You must be signed in to change notification settings - Fork 1
/
DSFinterp1DFit.py
194 lines (182 loc) · 8.82 KB
/
DSFinterp1DFit.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
# 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 & CSNS, Institute of High Energy Physics, CAS
# SPDX - License - Identifier: GPL - 3.0 +
# pylint: disable=no-init,invalid-name
'''
@author Jose Borreguero, NScD
@date October 06, 2013
'''
import numpy
import scipy.interpolate
from mantid.api import IFunction1D, FunctionFactory
from mantid.simpleapi import mtd
from mantid import logger
# pylint: disable=too-many-instance-attributes
class DSFinterp1DFit(IFunction1D):
_RegressionTypes = None
_minWindow = None
_InputWorkspaces = None
_LoadErrors = None
_WorkspaceIndex = None
_ParameterValues = None
_fmin = None
_fmax = None
_LocalRegression = None
_RegressionType = None
_RegressionWindow = None
_xvalues = None
_channelgroup = None
def category(self):
return 'QuasiElastic'
def init(self):
'''Declare parameters and attributes that participate in the fitting'''
# Active fitting parameters
self.declareParameter('Intensity', 1.0, 'Intensity')
self.declareParameter('TargetParameter', 1.0,
'Target value of the structure factor parameter')
self.declareAttribute('InputWorkspaces', '')
self.declareAttribute('LoadErrors', False)
self.declareAttribute('WorkspaceIndex', 0)
self.declareAttribute('ParameterValues', '')
self.declareAttribute('LocalRegression', True)
self.declareAttribute('RegressionType', 'quadratic')
self.declareAttribute('RegressionWindow', 6)
# "private" attributes associated to the declare function attributes
self._InputWorkspaces = None
self._LoadErrors = None
self._WorkspaceIndex = None
self._ParameterValues = None
self._fmin = None
self._fmax = None
self._LocalRegression = None
self._RegressionType = None
# valid syntaxfor python >= 2.6
self._RegressionTypes = set(['linear', 'quadratic'])
self._RegressionWindow = None
self._minWindow = {'linear': 3, 'quadratic': 4}
# channelgroup to interpolate values
self._channelgroup = None
self._xvalues = None # energies of the channels
def setAttributeValue(self, name, value):
if name == "InputWorkspaces":
self._InputWorkspaces = value.split()
if ',' in value:
self._InputWorkspaces = [x.strip() for x in value.split(',')]
elif name == 'LoadErrors':
self._LoadErrors = bool(value)
elif name == 'WorkspaceIndex':
self._WorkspaceIndex = int(value)
elif name == 'ParameterValues':
self._ParameterValues = []
self._fmin = 0.0
self._fmax = 0.0
if value:
self._ParameterValues = [float(f) for f in value.split()]
self._fmin = min(self._ParameterValues)
self._fmax = max(self._ParameterValues)
elif name == 'LocalRegression':
self._LocalRegression = bool(value)
elif name == 'RegressionType':
self._RegressionType = value.lower()
elif name == 'RegressionWindow':
self._RegressionWindow = value
def validateParams(self):
'''Check parameters within expected range'''
intensity = self.getParameterValue('Intensity')
if intensity <= 0:
message = 'Parameter Intensity in DSFinterp1DFit must be positive. Got {0} instead'.format(
intensity)
logger.error(message)
return None
f = self.getParameterValue('TargetParameter')
if f < self._fmin or f > self._fmax:
message = 'TargetParameter {0} is out of bounds [{1}, {2}]. Applying penalty...'.format(
f, self._fmin, self._fmax)
logger.error(message)
return None
return {'Intensity': intensity, 'TargetParameter': f}
def function1D(self, xvals):
''' Fit using the interpolated structure factor '''
p = self.validateParams()
if not p:
# return zeros if parameters not valid
return numpy.zeros(len(xvals), dtype=float)
# The first time the function is called requires initialization of the interpolator
if self._channelgroup is None:
# Check consistency of the input
# check InputWorkspaces have at least the workspace index
for w in self._InputWorkspaces:
if mtd[w].getNumberHistograms() <= self._WorkspaceIndex:
message = 'Numer of histograms in Workspace {0} does not allow for workspace index {1}'.format(
w, self._WorkspaceIndex)
logger.error(message)
raise IndexError(message)
# check number of input workspaces and parameters is the same
if len(self._ParameterValues) != len(self._InputWorkspaces):
message = 'Number of InputWorkspaces and ParameterValues should be the same.' +\
' Found {0} and {1}, respectively'.format(
len(self._InputWorkspaces), len(self._ParameterValues))
logger.error(message)
raise ValueError(message)
# check the regression type is valid
if self._RegressionType not in self._RegressionTypes:
message = 'Regression type {0} not implemented. choose one of {1}'.format(self._RegressionType,
', '.join(self._RegressionTypes))
logger.error(message)
raise NotImplementedError(message)
# check the regression window is appropriate for the regression type selected
if self._RegressionWindow < self._minWindow[self._RegressionType]:
message = 'RegressionWindow must be equal or bigger than ' +\
'{0} for regression type {1}'.format(
self._minWindow[self._RegressionType], self._RegressionType)
logger.error(message)
raise ValueError(message)
# Initialize the energies of the channels with the first of the input workspaces
self._xvalues = numpy.copy(
mtd[self._InputWorkspaces[0]].dataX(self._WorkspaceIndex))
if len(self._xvalues) == 1 + len(mtd[self._InputWorkspaces[0]].dataY(self._WorkspaceIndex)):
# Deal with histogram data
self._xvalues = (self._xvalues[1:]+self._xvalues[:-1])/2.0
# Initialize the channel group
nf = len(self._ParameterValues)
# Load the InputWorkspaces into a group of dynamic structure factors
from dsfinterp.dsf import Dsf
from dsfinterp.dsfgroup import DsfGroup
dsfgroup = DsfGroup()
for idsf in range(nf):
dsf = Dsf()
dsf.SetIntensities(
mtd[self._InputWorkspaces[idsf]].dataY(self._WorkspaceIndex))
dsf.errors = None # do not incorporate error data
if self._LoadErrors:
dsf.SetErrors(mtd[self._InputWorkspaces[idsf]].dataE(
self._WorkspaceIndex))
dsf.SetFvalue(self._ParameterValues[idsf])
dsfgroup.InsertDsf(dsf)
# Create the interpolator
from dsfinterp.channelgroup import ChannelGroup
self._channelgroup = ChannelGroup()
self._channelgroup.InitFromDsfGroup(dsfgroup)
if self._LocalRegression:
self._channelgroup.InitializeInterpolator(
running_regr_type=self._RegressionType, windowlength=self._RegressionWindow)
else:
self._channelgroup.InitializeInterpolator(windowlength=0)
# channel group has been initialized, so evaluate the interpolator
dsf = self._channelgroup(p['TargetParameter'])
# Linear interpolation between the energies of the channels and the xvalues we require
# NOTE: interpolator evaluates to zero for any of the xvals outside of the domain defined by self._xvalues
intensities_interpolator = scipy.interpolate.interp1d(
self._xvalues, p['Intensity']*dsf.intensities, kind='linear')
return intensities_interpolator(xvals) # can we pass by reference?
# Required to have Mantid recognize the new function
# pylint: disable=unused-import
try:
import dsfinterp # noqa
FunctionFactory.subscribe(DSFinterp1DFit)
except ImportError:
logger.debug('Failed to subscribe fit function DSFinterp1DFit. '
+ 'Python package dsfinterp may be missing (https://pypi.python.org/pypi/dsfinterp)')