source: sasview/park-1.2.1/park/peaks.py @ 0353540

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 0353540 was 3570545, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Adding park Part 2

  • Property mode set to 100644
File size: 9.2 KB
Line 
1# This program is public domain
2""" 
3A park model implementing a multipeak fitter.
4
5*** WARNING *** this example was used to inform the design process,
6and has not yet been updated to correspond to the current implementation.
7Do not use this as a tutorial.
8
9This is an example model showing how to put together a multi-part
10fit objective function.
11
12:group main: Peaks
13:group peaks: Gaussian, Lorentzian, Voigt
14:group backgrounds: Constant, Quadratic, Linear
15"""
16import numpy
17import park
18
19class Peaks(park.Model): 
20    """Peak fitter"""
21    # Increment the minor version when adding new functions
22    # Increment the major version if refactoring the class
23    __version__ = "0.9"
24    peak_types = {}
25
26    def __init__(self, datafile=None):
27        """
28        Create a new multipeak model.  If datafile is present it will
29        be loaded directly.  If not, then it can be set later with
30        self.data.load(datafile), or replaced with a specialized
31        data loader for the particular instrument.
32        """
33        # The initializer creates self.parameterset which will hold
34        # the set of fitting parameters for all models.  It will also
35        # define self.data = park.data.Data1D() which computes the
36        # residuals given an eval function.
37        park.Model.__init__(self,filename=datafile)
38        self.parameterset = park.ParameterSet()
39        self.peaks = {}
40        self.peaknum = 0
41
42    def eval(self, x):
43        """
44        Returns the raw theory value at x, unconvoluted and unweighted.
45       
46        The residuals will be calculated by park.Model.residuals by
47        calling
48        """
49        x = self.data[0]
50        y = numpy.empty(x.shape)
51        for name,peak in self.peaks.itervalues():
52            y += peak(x)
53        return y
54       
55    def add_peak(self, type, name=None, **kw):
56        """
57        Add a peak to the model.
58       
59        The name of the peak is used to distinguish it from other peaks
60        in the model, and appears as part of the parameter name in
61        constraint expressions.  For example, if name is 'P3' and this
62        is part of model 'M1' then the sigma parameter of a gaussian
63        peak would be 'M1.P3.sigma'.
64
65        The peak can be specified either by peak type and initial arguments
66        or by the peak itself, which is a function of x with a parameters
67        attribute.
68       
69        Available peak types include:
70            gaussian(scale=1, center=0, sigma=1)
71                = scale exp ( -0.5 (x-center)**2 / sigma**2 )
72            lorentzian(scale=1, center=0, gamma=1)
73                = scale/pi gamma / ((x-center)**2 + gamma**2))
74            voigt(scale=1, center=0, sigma=1, gamma=1)
75                = scale (lorentzian(mu,gamma) * gaussian(mu,sigma))(x-center)
76                where * represents convolution
77       
78        Available background functions include:
79            constant(C=0)
80                = C
81            slope(C=0, B=1)
82                = B x + C
83            quadratic(A=1, B=0, C=0)
84                = A x**2 + B x + C
85
86        Additional peak types can be registered.  The complete dictionary
87        of available types is in Peaks.peak_types.
88        """
89        self.peaknum
90        if name == None: name = 'P'+str(self.peaknum)
91        if isinstance(type,basestring):
92            peak = self.peak_types[type](**kw)
93        else:
94            type = peak
95        self.peaks[name] = peak
96        self.parameters[name] = peak.parameters
97        setattr(self,name,peak)
98
99    def remove_peak(self, name):
100        """
101        Remove the named peak from the model.
102        """
103        del self.peaks[name]
104        del self.parameterset[name]
105
106    @classmethod
107    def register_peak_generator(cls, type, factory):
108        """
109        Register a new peak generator.  This will return a callable
110        object with a parameter set that can be used in the peak
111        fitting class.
112
113        The individual peak functions a number of attributes:
114       
115        __name__
116            the name to be displayed in the list of peaks; the
117            default is __class__.__name__
118        __factory__
119            the name of a factory function in the python
120            namespace which returns a new peak.  These parameters
121            will be set from the stored parameter list when the
122            peak model is reloaded.  This can be omitted if it
123            is simply the module and name of the class.
124        __doc__
125            the tool tip to be displayed with the peak
126        __call__(x)
127            evaluate the peak function at x
128        """
129        cls.peak_types[type]=factory
130
131
132
133# ================== Peaks and backgrounds ===================
134class Gaussian(object):
135    """Gaussian peak: scale exp ( -0.5 (x-center)**2 / sigma**2 )"""
136    def __init__(self, name, scale=1., center=0, sigma=1):
137        park.define_parameters(name, ['scale','center','sigma'],
138                               scale=scale,center=center,sigma=sigma)
139        self.parameters['scale'].deriv = self.dscale
140        self.parameters['center'].deriv = self.dcenter
141        self.parameters['sigma'].deriv = self.dsigma
142    def __call__(self, x):
143        return self.scale * numpy.exp(-((x-self.center)/self.sigma)**2)
144    def dscale(self, x):
145        return self(x)/self.scale
146    def dcenter(self, x):
147        return 2*(x-self.center)/self.sigma*self(x)
148    def dsigma(self, x):
149        return 2*(x-self.center)/self.sigma**3*self(x)
150
151class Lorentzian(object):
152    """Lorentzian peak (HWHM): scale/pi  gamma/((x-center)**2 + gamma**2)"""
153    def __init__(self, name, scale=1., center=0, gamma=1):
154        park.define_parameters(name, ['scale','center','gamma'],
155                               scale=scale,center=center,gamma=gamma)
156        self.parameters['scale'].deriv = self.dscale
157        self.parameters['center'].deriv = self.dcenter
158        self.parameters['gamma'].deriv = self.dgamma
159    def __call__(self, x):
160        return self.scale/numpy.pi * self.gamma/((x-self.center)**2 + self.gamma**2)
161    def dscale(self, x):
162        return self(x)/self.scale
163    def dcenter(self, x):
164        return 2*(x-self.center)/((x-self.center)**2+self.gamma**2)*self(x)
165    def dgamma(self, x):
166        return self(x)*(1/self.gamma - self.gamma/((x-self.center)**2 + self.gamma**2))
167
168class Voigt(object):
169    """Voigt peak (HWHM,sigma): A [G(sigma) * L(gamma)](x-center)"""
170    def __init__(self, name, scale=1, center=0, sigma=1, gamma=1):
171        park.define_parameters(name, ['scale','center','sigma','gamma'],
172                               scale=scale, center=center, sigma=sigma,
173                               gamma=gamma)
174    def __call__(self, x):
175        return self.scale*voigt(x-self.center, 
176                                sigma=self.sigma, gamma=self.gamma)
177
178class Quadratic(object):
179    """Quadratic background: A x**2 + B x + C"""
180    name = "Background: quadratic"
181    def __init__(self, name, A=1, B=0, C=0):
182        park.define_parameters(name, ['A', 'B', 'C'], A=A, B=B, C=C)
183        self.parameters['A'].deriv = self.dA
184        self.parameters['B'].deriv = self.dB
185        self.parameters['C'].deriv = self.dC
186    def __call__(self, x):
187        return numpy.polyval([self.a.value, self.b.value, self.c.value],x)
188    def dA(self, x): return x**2
189    def dB(self, x): return x
190    def dC(self, x): return numpy.ones(x.shape)
191
192class Linear(object):
193    """Linear background: B x**2 + C"""
194    name = "Background: linear"
195    def __init__(self, name, B=0, C=0):
196        park.define_parameters(name, ['B', 'C'], B=B, C=C)
197        self.parameters['B'].deriv = self.dB
198        self.parameters['C'].deriv = self.dC
199    def __call__(self, x):
200        return numpy.polyval([self.a.value, self.b.value, self.c.value],x)
201    def dB(self, x): return x
202    def dC(self, x): return numpy.ones(x.shape)
203
204class Constant(object):
205    """Constant background: C"""
206    name = "Background: constant"
207    def __init__(self, name, C=0):
208        park.define_parameters(name, ['C'], C=C)
209        self.parameters['C'].deriv = self.dC
210    def __call__(self, x):
211        return numpy.polyval([self.a.value, self.b.value, self.c.value],x)
212    def dC(self, x): return numpy.ones(x.shape)
213
214
215def voigt(x, sigma, gamma):
216    """
217    Return the voigt function, which is the convolution of a Lorentz
218    function with a Gaussian.
219   
220    :Parameters:
221     gamma : real
222      The half-width half-maximum of the Lorentzian
223     sigma : real
224      The 1-sigma width of the Gaussian, which is one standard deviation.
225
226    Ref: W.I.F. David, J. Appl. Cryst. (1986). 19, 63-64
227   
228    Note: adjusted to use stddev and HWHM rather than FWHM parameters
229    """
230    # wofz function = w(z) = Fad[d][e][y]eva function = exp(-z**2)erfc(-iz)
231    from scipy.special import wofz
232    z = (x+1j*gamma)/(sigma*sqrt(2))
233    V = wofz(z)/(sqrt(2*pi)*sigma)
234    return V.real
235
236
237def init():
238    """
239    Register peak types with the Peaks model.
240    """
241    for cls in [Gaussian, Lorentzian, Voigt, Quadratic, Linear, Constant]:
242        Peaks.register_peak_generator(cls.__name__,cls)
243init()
244
245# ========================================
246
247def test():
248    x = numpy.linspace(-5,5,10)
249   
250    pass
251
252if __name__ == "__main__": test()
Note: See TracBrowser for help on using the repository browser.