Skip to content

Commit 63a0230

Browse files
committed
Adding code supporting Pavol's recipe interface. It does pretty much what is
shown in the new pinterface.py example. Have fun!
1 parent 3663e49 commit 63a0230

3 files changed

Lines changed: 318 additions & 11 deletions

File tree

‎diffpy/srfit/fit/recipe.py‎

Lines changed: 240 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,240 @@
1+
#!/usr/bin/env python
2+
########################################################################
3+
#
4+
# diffpy.srfit by DANSE Diffraction group
5+
# Simon J. L. Billinge
6+
# (c) 2011 Trustees of the Columbia University
7+
# in the City of New York. All rights reserved.
8+
#
9+
# File coded by: Chris Farrow
10+
#
11+
# See AUTHORS.txt for a list of people who contributed.
12+
# See LICENSE.txt for license information.
13+
#
14+
########################################################################
15+
"""Recipe class for creating and organizing fit objects."""
16+
17+
from diffpy.srfit.adapters import adapt
18+
from diffpy.srfit.adapters.nodes import Node
19+
from diffpy.srfit.fit import Parameter, Var
20+
from diffpy.srfit.util import getVaried, isAdapter
21+
from diffpy.srfit.util.cachemanager import CacheManager
22+
23+
__all__ = ["Recipe"]
24+
25+
class Recipe(object):
26+
"""Object for generating equations and managing their parameters.
27+
28+
The purpose of this class is to generate equations that can be used as part
29+
of a fit. Various parts of the interface accepts string equations, which
30+
can simplify the creation of fit complexes. See the methods and examples
31+
for use.
32+
33+
Attributes
34+
Each top-level object passed or extracted during creation is accessible as
35+
an attribute of the equation class. These objects are adapted, and so have
36+
the parameter interface.
37+
38+
"""
39+
40+
def __init__(self, **kw):
41+
"""Initialize the equation with a string.
42+
43+
kw -- Named arguments that define the local namespace.
44+
45+
"""
46+
# Adapt and store keywords
47+
self._namespace = {}
48+
# FIXME - rebuild this when pickled.
49+
self._adapters = {}
50+
self.register(**kw)
51+
# FIXME - need objects for accessing constrained parameters.
52+
return
53+
54+
# FIXME - do __setattr__ to set values directly.
55+
56+
def __getattr__(self, name):
57+
"""Get an attribute from the namespace."""
58+
try:
59+
return self._namespace[name]
60+
except KeyError:
61+
msg = "'%s' has no attribute '%s'" % (self, name)
62+
raise AttributeError(msg)
63+
64+
def register(self, **kw):
65+
"""Register new objects.
66+
67+
Object adapters will be available as attributes.
68+
69+
"""
70+
for name, val in kw.items():
71+
adapter = adapt(val, name)
72+
self._namespace[name] = adapter
73+
self._adapters[id(val)] = adapter
74+
return
75+
76+
def _getAdapter(self, obj):
77+
"""Get an adapter for an object.
78+
79+
This retrieves a cached adapter or creates a new one and caches it.
80+
81+
"""
82+
key = id(obj)
83+
adapter = self._adapters.get(key)
84+
if adapter is None:
85+
adapter = adapt(obj)
86+
self._adapters[key] = adapter
87+
return adapter
88+
89+
# XXX Should this call register?
90+
def vars(self, **kw):
91+
"""Create variables and add them to the namespace."""
92+
for name, val in kw.items():
93+
v = Var(name, val)
94+
self._namespace[name] = v
95+
return
96+
97+
def var(self, name, val = None):
98+
"""Create and return variables. It gets added to the namespace."""
99+
v = Var(name, val)
100+
self._namespace[name] = v
101+
return v
102+
103+
def pars(self, **kw):
104+
"""Create parameters and add them to the namespace."""
105+
for name, val in kw.items():
106+
v = Parameter(name, val)
107+
self._namespace[name] = v
108+
return
109+
110+
def par(self, name, val = None):
111+
"""Create and return parameter. It gets added to the namespace."""
112+
p = Parameter(name, val)
113+
self._namespace[name] = p
114+
return p
115+
116+
def constrain(self, target, *args):
117+
"""Get a constraint creator.
118+
119+
target -- Target of the constraint.
120+
args -- Modfying objects.
121+
122+
The args and target define a dependency relationship. When the args are
123+
modified they indicate that the constraint needs to be reevaluated.
124+
This reevaluation occurs whenever the target is retrieved from the
125+
recipe (e.g. recipe.target.value evaluates the constraint).
126+
127+
Returns a callable with which to define the constraint. The constraint
128+
callable accepts a function of no arguments or a string. If a function
129+
is passed, it relies on its globals for evaluation. A string constraint
130+
is evaluated within the namespace defined by the target and the args.
131+
132+
"""
133+
# Make sure everything is registered
134+
tadapter = self._getAdapter(target)
135+
# Check the arguments
136+
targs = [self._getAdapter(arg) for arg in args]
137+
138+
nsset = set(self._namespace.values())
139+
msg = ""
140+
if tadapter not in nsset:
141+
msg = "'%s' is not registered" % target
142+
for targ in targs:
143+
if targ not in nsset:
144+
msg = "'%s' is not registered" % target
145+
if msg:
146+
raise ValueError(msg)
147+
148+
cgen = ConstraintBuilder(self, tadapter, targs)
149+
return cgen
150+
151+
def build(self, streq):
152+
"""Build an equation from a string.
153+
154+
The string will be evaluated in the namespace of registered objects,
155+
parameters and variables.
156+
157+
"""
158+
ns = dict(self._namespace)
159+
eq = eval(streq, ns)
160+
eq._makeFunction([], self._namespace)
161+
return eq
162+
163+
# End class Recipe
164+
165+
class ConstraintBuilder(object):
166+
"""Builder used by Recipes to build constraints."""
167+
168+
def __init__(self, recipe, target, args):
169+
self._recipe = recipe
170+
self._target = target
171+
self._args = args
172+
return
173+
174+
def __call__(self, con):
175+
176+
if isinstance(con, str):
177+
# Create a string getter and adapt it
178+
ns = dict( (arg.name, arg) for arg in self._args )
179+
ns[self._target.name] = self._target
180+
getter = StringGetter(con)
181+
adapter = adapt(ns, "stringpar", getter)
182+
else:
183+
getter = ConGetter()
184+
adapter = adapt(con, "funcpar", getter)
185+
186+
# Make the adapter a function. This sets up the functional relationship
187+
# between the adapter and the arguments.
188+
adapter._makeFunction(self._args, {})
189+
190+
# Create a proxy parameter for the target. This will hold the
191+
# constraint and notify the target whenever its value is invalid.
192+
proxy = Parameter(self._target.name + "_proxy")
193+
proxy._cache = self._target._cache
194+
proxy._cache.addNode(proxy)
195+
proxy.constrain(adapter)
196+
197+
# FIXME - associate the constraint with the target in the recipe so it
198+
# can be looked up, etc.
199+
200+
return
201+
202+
# End class ConstraintBuilder
203+
204+
class ConGetter(object):
205+
"""Getter for function constraints."""
206+
207+
def __call__(self, func):
208+
"""Get the function value.
209+
210+
func -- The function to call. This takes no arguments.
211+
212+
Returns the value of the function.
213+
214+
"""
215+
return func()
216+
217+
# End class ConGetter
218+
219+
class StringGetter(object):
220+
"""Getter for string equations."""
221+
222+
def __init__(self, eqstr):
223+
self.eqstr = eqstr
224+
return
225+
226+
def __call__(self, ns):
227+
"""Get the value.
228+
229+
ns -- Dictionary holding the namespace of reference objects.
230+
231+
Returns None, since we don't know what eqstr does.
232+
233+
"""
234+
valspace = dict((name, ref.value) for name, ref in ns.items())
235+
exec self.eqstr in valspace
236+
return None
237+
238+
# End class StringGetter
239+
240+
__id__ = "$Id$"

‎doc/examples/crystalpdf.py‎

Lines changed: 15 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,8 @@
2424
2525
"""
2626

27+
from numpy import pi
28+
2729
from diffpy.Structure import Structure
2830
from diffpy.srreal.pdfcalculator import PDFCalculator
2931
# This import statement registers pdf-related adapters
@@ -39,20 +41,19 @@ def main(ciffile, datname):
3941
# just for this. We tell 'loadProfile' to use this parser by specifying the
4042
# 'PDF' format.
4143
profile = loadProfile(datname, "PDF")
42-
# Now we can make good use of the profile and set the fit range.
43-
profile.setRange(xmax = 20)
44-
r, gr, dgr = profile
4544

4645
# Create a PDFCalculator. We can use information from the metadata loaded
4746
# from file to initialize it. We can't pass the whole meta dictionary since
4847
# it contains entries that PDFCalculator can't use. In addition, passing a
4948
# non-zero qmin value to PDFCalculator creates a "bad" PDF.
50-
calc = PDFCalculator(qmax = profile.meta["qmax"])
49+
calc = PDFCalculator(qmax = profile.meta["qmax"],
50+
rmin=profile.x[0], rmax=20)
51+
calc.rstep = pi / calc.qmax
5152
calc.setScatteringFactorTableByType(profile.meta["stype"])
52-
rdat = r.get()
53-
calc.rmin = rdat[0]
54-
calc.rstep = rdat[1] - rdat[0]
55-
calc.rmax = rdat[-1] + 0.5 * calc.rstep
53+
54+
# Use the same grid as the profile
55+
profile.setPoints(calc.rgrid)
56+
r, gr, dgr = profile.pars
5657

5758
# Now we adapt the PDFCalculator so it can be used as a symbolic
5859
# calculation object. We must pass our PDFCalculator instance, and
@@ -92,7 +93,7 @@ def main(ciffile, datname):
9293
# Because our data is on the grid defined by 'r', and the fit is on the
9394
# grid defined by 'rcalc', we have to interpolate in order to compare them.
9495
# The numpy function interp has been adapted for this purpose.
95-
fiteq = interp_(r, rcalc, gcalc)
96+
fiteq = gcalc#interp_(r, rcalc, gcalc)
9697

9798
# Create the residual equation. Note that 'chi' creates a vector residual
9899
# that can be dotted into itself to generate 'chi^2'.
@@ -107,8 +108,10 @@ def main(ciffile, datname):
107108
res = residual(reseq, rest1)
108109

109110
# Optimize.
110-
from scipy.optimize import leastsq
111-
leastsq(res.vec, res.values)
111+
#from scipy.optimize import leastsq
112+
#leastsq(res.vec, res.values)
113+
from scipy.optimize import fmin
114+
fmin(res, res.values)
112115

113116
# Get the results
114117
results = FitResults(res)
@@ -125,6 +128,7 @@ def main(ciffile, datname):
125128
# Make the data and the recipe
126129
ciffile = "data/ni.cif"
127130
data = "data/ni-q27r100-neutron.gr"
131+
#data = "data/ni-q27r60-xray.gr"
128132
main(ciffile, data)
129133

130134
# End of file

‎doc/examples/pinterface.py‎

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
#!/usr/bin/env python
2+
3+
from diffpy.srreal.pdfcalculator import PDFCalculator
4+
from diffpy.Structure import Structure
5+
from diffpy.srfit.fit.recipe import Recipe
6+
import diffpy.srfit.pdf
7+
from diffpy.srfit import loadProfile, residual, FitResults
8+
9+
def main():
10+
data = "data/ni-q27r100-neutron.gr"
11+
profile = loadProfile(data, "PDF")
12+
pdfc = PDFCalculator()
13+
pdfc.qmax = profile.meta["qmax"]
14+
pdfc.setScatteringFactorTableByType(profile.meta["stype"])
15+
16+
stru = Structure(filename="./data/ni.cif")
17+
# XXX - made change here. Adapted objects must be registered with the
18+
# recipe. The reason for this is because these objects are held as adapters
19+
# inside the recipe. Changing this would require an overhaul of this
20+
# branch, which is premature at this point.
21+
recipe = Recipe(pdfc=pdfc,stru=stru,profile=profile)
22+
recipe.vars(a=stru.lattice.a, Uiso=0.003, scale=1, qdamp=0.01, delta2=5)
23+
24+
# Link the fit range between the profile and the calculator
25+
recipe.pars(rmin = 0.1, rmax = 20, rstep=0.05)
26+
recipe.constrain(pdfc, recipe.rmin)('pdfc.rmin = rmin')
27+
recipe.constrain(pdfc, recipe.rmax)('pdfc.rmax = rmax')
28+
recipe.constrain(pdfc, recipe.rstep)('pdfc.rstep = rstep')
29+
# XXX pdfc is not a parameter, but it works very well here
30+
recipe.constrain(profile, pdfc)('profile.setPoints(pdfc.rgrid)')
31+
32+
# Set of other constraints
33+
recipe.constrain(stru, recipe.a)('stru.lattice.setLatPar(a=a, b=a, c=a)')
34+
recipe.constrain(stru, recipe.Uiso)('stru.Uisoequiv = Uiso')
35+
recipe.constrain(pdfc, recipe.scale)('pdfc.scale = scale')
36+
recipe.constrain(pdfc, recipe.qdamp)('pdfc.qdamp = qdamp')
37+
recipe.constrain(pdfc, recipe.delta2)('pdfc.delta2 = delta2')
38+
39+
# Make the equation to be fit
40+
gdiff = recipe.build('(profile.y - pdfc(stru)[1]) / profile.dy')
41+
res = residual(gdiff)
42+
43+
from scipy.optimize import leastsq
44+
recipe.rstep.value=0.1
45+
leastsq(res.vec, res.values)
46+
recipe.rstep.value=0.01
47+
leastsq(res.vec, res.values)
48+
49+
results = FitResults(res)
50+
results.show()
51+
52+
# and plot
53+
from pylab import plot, show
54+
r, gr, dgr = profile
55+
r, rcalc = pdfc(stru)
56+
plot(r, gr, 'bo')
57+
plot(r, rcalc, color='r', linestyle='-')
58+
diff = gr - rcalc
59+
plot(r, diff + 1.15*min(gr), color='g', linestyle='-')
60+
show()
61+
62+
if __name__ == "__main__":
63+
main()

0 commit comments

Comments
 (0)