1 | # This program is public domain |
---|
2 | """ |
---|
3 | Parameter expression evaluator. |
---|
4 | |
---|
5 | For systems in which constraints are expressed as string expressions rather |
---|
6 | than python code, :func:`compile_constraints` can construct an expression |
---|
7 | evaluator that substitutes the computed values of the expressions into the |
---|
8 | parameters. |
---|
9 | |
---|
10 | The compiler requires a symbol table, an expression set and a context. |
---|
11 | The symbol table maps strings containing fully qualified names such as |
---|
12 | 'M1.c[3].full_width' to parameter objects with a 'value' property that |
---|
13 | can be queried and set. The expression set maps symbol names from the |
---|
14 | symbol table to string expressions. The context provides additional symbols |
---|
15 | for the expressions in addition to the usual mathematical functions and |
---|
16 | constants. |
---|
17 | |
---|
18 | The expressions are compiled and interpreted by python, with only minimal |
---|
19 | effort to make sure that they don't contain bad code. The resulting |
---|
20 | constraints function returns 0 so it can be used directly in a fit problem |
---|
21 | definition. |
---|
22 | |
---|
23 | Extracting the symbol table from the model depends on the structure of the |
---|
24 | model. If fitness.parameters() is set correctly, then this should simply |
---|
25 | be a matter of walking the parameter data, remembering the path to each |
---|
26 | parameter in the symbol table. For compactness, dictionary elements should |
---|
27 | be referenced by .name rather than ["name"]. Model name can be used as the |
---|
28 | top level. |
---|
29 | |
---|
30 | Getting the parameter expressions applied correctly is challenging. |
---|
31 | The following monkey patch works by overriding model_update in FitProblem |
---|
32 | so that after setp(p) is called and, the constraints expression can be |
---|
33 | applied before telling the underlying fitness function that the model |
---|
34 | is out of date:: |
---|
35 | |
---|
36 | # Override model update so that parameter constraints are applied |
---|
37 | problem._model_update = problem.model_update |
---|
38 | def model_update(): |
---|
39 | constraints() |
---|
40 | problem._model_update() |
---|
41 | problem.model_update = model_update |
---|
42 | |
---|
43 | Ideally, this interface will change |
---|
44 | """ |
---|
45 | from __future__ import print_function |
---|
46 | |
---|
47 | import math |
---|
48 | import re |
---|
49 | |
---|
50 | # simple pattern which matches symbols. Note that it will also match |
---|
51 | # invalid substrings such as a3...9, but given syntactically correct |
---|
52 | # input it will only match symbols. |
---|
53 | _symbol_pattern = re.compile('([a-zA-Z_][a-zA-Z_0-9.]*)') |
---|
54 | |
---|
55 | def _symbols(expr,symtab): |
---|
56 | """ |
---|
57 | Given an expression string and a symbol table, return the set of symbols |
---|
58 | used in the expression. Symbols are only returned once even if they |
---|
59 | occur multiple times. The return value is a set with the elements in |
---|
60 | no particular order. |
---|
61 | |
---|
62 | This is the first step in computing a dependency graph. |
---|
63 | """ |
---|
64 | matches = [m.group(0) for m in _symbol_pattern.finditer(expr)] |
---|
65 | return set([symtab[m] for m in matches if m in symtab]) |
---|
66 | |
---|
67 | def _substitute(expr,mapping): |
---|
68 | """ |
---|
69 | Replace all occurrences of symbol s with mapping[s] for s in mapping. |
---|
70 | """ |
---|
71 | # Find the symbols and the mapping |
---|
72 | matches = [(m.start(),m.end(),mapping[m.group(1)]) |
---|
73 | for m in _symbol_pattern.finditer(expr) |
---|
74 | if m.group(1) in mapping] |
---|
75 | |
---|
76 | # Split the expression in to pieces, with new symbols replacing old |
---|
77 | pieces = [] |
---|
78 | offset = 0 |
---|
79 | for start,end,text in matches: |
---|
80 | pieces += [expr[offset:start],text] |
---|
81 | offset = end |
---|
82 | pieces.append(expr[offset:]) |
---|
83 | |
---|
84 | # Join the pieces and return them |
---|
85 | return "".join(pieces) |
---|
86 | |
---|
87 | def _find_dependencies(symtab, exprs): |
---|
88 | """ |
---|
89 | Returns a list of pair-wise dependencies from the parameter expressions. |
---|
90 | |
---|
91 | For example, if p3 = p1+p2, then find_dependencies([p1,p2,p3]) will |
---|
92 | return [(p3,p1),(p3,p2)]. For base expressions without dependencies, |
---|
93 | such as p4 = 2*pi, this should return [(p4, None)] |
---|
94 | """ |
---|
95 | deps = [(target,source) |
---|
96 | for target,expr in exprs.items() |
---|
97 | for source in _symbols_or_none(expr,symtab)] |
---|
98 | return deps |
---|
99 | |
---|
100 | # Hack to deal with expressions without dependencies --- return a fake |
---|
101 | # dependency of None. |
---|
102 | # The better solution is fix order_dependencies so that it takes a |
---|
103 | # dictionary of {symbol: dependency_list}, for which no dependencies |
---|
104 | # is simply []; fix in parameter_mapping as well |
---|
105 | def _symbols_or_none(expr,symtab): |
---|
106 | syms = _symbols(expr,symtab) |
---|
107 | return syms if len(syms) else [None] |
---|
108 | |
---|
109 | def _parameter_mapping(pairs): |
---|
110 | """ |
---|
111 | Find the parameter substitution we need so that expressions can |
---|
112 | be evaluated without having to traverse a chain of |
---|
113 | model.layer.parameter.value |
---|
114 | """ |
---|
115 | left,right = zip(*pairs) |
---|
116 | pars = list(sorted(p for p in set(left+right) if p is not None)) |
---|
117 | definition = dict( ('P%d'%i,p) for i,p in enumerate(pars) ) |
---|
118 | # p is None when there is an expression with no dependencies |
---|
119 | substitution = dict( (p,'P%d.value'%i) |
---|
120 | for i,p in enumerate(sorted(pars)) |
---|
121 | if p is not None) |
---|
122 | return definition, substitution |
---|
123 | |
---|
124 | def no_constraints(): |
---|
125 | """ |
---|
126 | This parameter set has no constraints between the parameters. |
---|
127 | """ |
---|
128 | pass |
---|
129 | |
---|
130 | def compile_constraints(symtab, exprs, context={}): |
---|
131 | """ |
---|
132 | Build and return a function to evaluate all parameter expressions in |
---|
133 | the proper order. |
---|
134 | |
---|
135 | Input: |
---|
136 | |
---|
137 | *symtab* is the symbol table for the model: { 'name': parameter } |
---|
138 | |
---|
139 | *exprs* is the set of computed symbols: { 'name': 'expression' } |
---|
140 | |
---|
141 | *context* is any additional context needed to evaluate the expression |
---|
142 | |
---|
143 | Return: |
---|
144 | |
---|
145 | updater function which sets parameter.value for each expression |
---|
146 | |
---|
147 | Raises: |
---|
148 | |
---|
149 | AssertionError - model, parameter or function is missing |
---|
150 | |
---|
151 | SyntaxError - improper expression syntax |
---|
152 | |
---|
153 | ValueError - expressions have circular dependencies |
---|
154 | |
---|
155 | This function is not terribly sophisticated, and it would be easy to |
---|
156 | trick. However it handles the common cases cleanly and generates |
---|
157 | reasonable messages for the common errors. |
---|
158 | |
---|
159 | This code has not been fully audited for security. While we have |
---|
160 | removed the builtins and the ability to import modules, there may |
---|
161 | be other vectors for users to perform more than simple function |
---|
162 | evaluations. Unauthenticated users should not be running this code. |
---|
163 | |
---|
164 | Parameter names are assumed to contain only _.a-zA-Z0-9#[] |
---|
165 | |
---|
166 | Both names are provided for inverse functions, e.g., acos and arccos. |
---|
167 | |
---|
168 | Should try running the function to identify syntax errors before |
---|
169 | running it in a fit. |
---|
170 | |
---|
171 | Use help(fn) to see the code generated for the returned function fn. |
---|
172 | dis.dis(fn) will show the corresponding python vm instructions. |
---|
173 | """ |
---|
174 | |
---|
175 | # Sort the parameters in the order they need to be evaluated |
---|
176 | deps = _find_dependencies(symtab, exprs) |
---|
177 | if deps == []: return no_constraints |
---|
178 | order = order_dependencies(deps) |
---|
179 | |
---|
180 | |
---|
181 | # Rather than using the full path to the parameters in the parameter |
---|
182 | # expressions, instead use Pn, and substitute Pn.value for each occurrence |
---|
183 | # of the parameter in the expression. |
---|
184 | names = list(sorted(symtab.keys())) |
---|
185 | parameters = dict(('P%d'%i, symtab[k]) for i,k in enumerate(names)) |
---|
186 | mapping = dict((k, 'P%d.value'%i) for i,k in enumerate(names)) |
---|
187 | |
---|
188 | |
---|
189 | # Initialize dictionary with available functions |
---|
190 | global_context = {} |
---|
191 | global_context.update(math.__dict__) |
---|
192 | global_context.update(dict(arcsin=math.asin,arccos=math.acos, |
---|
193 | arctan=math.atan,arctan2=math.atan2)) |
---|
194 | global_context.update(context) |
---|
195 | global_context.update(parameters) |
---|
196 | global_context['id'] = id |
---|
197 | local_context = {} |
---|
198 | |
---|
199 | # Define the constraints function |
---|
200 | assignments = ["=".join((p,exprs[p])) for p in order] |
---|
201 | code = [_substitute(s, mapping) for s in assignments] |
---|
202 | functiondef = """ |
---|
203 | def eval_expressions(): |
---|
204 | ''' |
---|
205 | %s |
---|
206 | ''' |
---|
207 | %s |
---|
208 | return 0 |
---|
209 | """%("\n ".join(assignments),"\n ".join(code)) |
---|
210 | |
---|
211 | #print("Function: "+functiondef) |
---|
212 | # CRUFT: python < 3.0; doc builder isn't allowing the following exec |
---|
213 | # https://stackoverflow.com/questions/4484872/why-doesnt-exec-work-in-a-function-with-a-subfunction/41368813#comment73790496_41368813 |
---|
214 | #exec(functiondef, global_context, local_context) |
---|
215 | eval(compile(functiondef, '<string>', 'exec'), global_context, local_context) |
---|
216 | retfn = local_context['eval_expressions'] |
---|
217 | |
---|
218 | # Remove garbage added to globals by exec |
---|
219 | global_context.pop('__doc__', None) |
---|
220 | global_context.pop('__name__', None) |
---|
221 | global_context.pop('__file__', None) |
---|
222 | global_context.pop('__builtins__') |
---|
223 | #print globals.keys() |
---|
224 | |
---|
225 | return retfn |
---|
226 | |
---|
227 | def order_dependencies(pairs): |
---|
228 | """ |
---|
229 | Order elements from pairs so that b comes before a in the |
---|
230 | ordered list for all pairs (a,b). |
---|
231 | """ |
---|
232 | #print "order_dependencies",pairs |
---|
233 | emptyset = set() |
---|
234 | order = [] |
---|
235 | |
---|
236 | # Break pairs into left set and right set |
---|
237 | # Note: pairs is array or list, so use "len(pairs) > 0" to check for empty. |
---|
238 | left,right = [set(s) for s in zip(*pairs)] if len(pairs) > 0 else ([],[]) |
---|
239 | while len(pairs) > 0: |
---|
240 | #print "within",pairs |
---|
241 | # Find which items only occur on the right |
---|
242 | independent = right - left |
---|
243 | if independent == emptyset: |
---|
244 | cycleset = ", ".join(str(s) for s in left) |
---|
245 | raise ValueError("Cyclic dependencies amongst %s"%cycleset) |
---|
246 | |
---|
247 | # The possibly resolvable items are those that depend on the independents |
---|
248 | dependent = set([a for a,b in pairs if b in independent]) |
---|
249 | pairs = [(a,b) for a,b in pairs if b not in independent] |
---|
250 | if pairs == []: |
---|
251 | resolved = dependent |
---|
252 | else: |
---|
253 | left,right = [set(s) for s in zip(*pairs)] |
---|
254 | resolved = dependent - left |
---|
255 | #print "independent",independent,"dependent",dependent,"resolvable",resolved |
---|
256 | order += resolved |
---|
257 | #print "new order",order |
---|
258 | order.reverse() |
---|
259 | return order |
---|
260 | |
---|
261 | # ========= Test code ======== |
---|
262 | def _check(msg,pairs): |
---|
263 | """ |
---|
264 | Verify that the list n contains the given items, and that the list |
---|
265 | satisfies the partial ordering given by the pairs in partial order. |
---|
266 | """ |
---|
267 | # Note: pairs is array or list, so use "len(pairs) > 0" to check for empty. |
---|
268 | left,right = zip(*pairs) if len(pairs) > 0 else ([],[]) |
---|
269 | items = set(left) |
---|
270 | n = order_dependencies(pairs) |
---|
271 | if set(n) != items or len(n) != len(items): |
---|
272 | n.sort() |
---|
273 | items = list(items); items.sort() |
---|
274 | raise ValueError("%s expect %s to contain %s for %s"%(msg,n,items,pairs)) |
---|
275 | for lo,hi in pairs: |
---|
276 | if lo in n and hi in n and n.index(lo) >= n.index(hi): |
---|
277 | raise ValueError("%s expect %s before %s in %s for %s"%(msg,lo,hi,n,pairs)) |
---|
278 | |
---|
279 | def test_deps(): |
---|
280 | import numpy as np |
---|
281 | |
---|
282 | # Null case |
---|
283 | _check("test empty",[]) |
---|
284 | |
---|
285 | # Some dependencies |
---|
286 | _check("test1",[(2,7),(1,5),(1,4),(2,1),(3,1),(5,6)]) |
---|
287 | _check("test1 renumbered",[(6,1),(7,3),(7,4),(6,7),(5,7),(3,2)]) |
---|
288 | _check("test1 numpy",np.array([(2,7),(1,5),(1,4),(2,1),(3,1),(5,6)])) |
---|
289 | |
---|
290 | # No dependencies |
---|
291 | _check("test2",[(4,1),(3,2),(8,4)]) |
---|
292 | |
---|
293 | # Cycle test |
---|
294 | pairs = [(1,4),(4,3),(4,5),(5,1)] |
---|
295 | try: |
---|
296 | n = order_dependencies(pairs) |
---|
297 | except ValueError: |
---|
298 | pass |
---|
299 | else: |
---|
300 | raise ValueError("test3 expect ValueError exception for %s"%(pairs,)) |
---|
301 | |
---|
302 | # large test for gross speed check |
---|
303 | A = np.random.randint(4000,size=(1000,2)) |
---|
304 | A[:,1] += 4000 # Avoid cycles |
---|
305 | _check("test-large",A) |
---|
306 | |
---|
307 | # depth tests |
---|
308 | k = 200 |
---|
309 | A = np.array([range(0,k),range(1,k+1)]).T |
---|
310 | _check("depth-1",A) |
---|
311 | |
---|
312 | A = np.array([range(1,k+1),range(0,k)]).T |
---|
313 | _check("depth-2",A) |
---|
314 | |
---|
315 | def test_expr(): |
---|
316 | import inspect, dis |
---|
317 | import math |
---|
318 | |
---|
319 | symtab = {'a.b.x':1, 'a.c':2, 'a.b':3, 'b.x':4} |
---|
320 | expr = 'a.b.x + sin(4*pi*a.c) + a.b.x/a.b' |
---|
321 | |
---|
322 | # Check symbol lookup |
---|
323 | assert _symbols(expr, symtab) == set([1,2,3]) |
---|
324 | |
---|
325 | # Check symbol rename |
---|
326 | assert _substitute(expr,{'a.b.x':'Q'}) == 'Q + sin(4*pi*a.c) + Q/a.b' |
---|
327 | assert _substitute(expr,{'a.b':'Q'}) == 'a.b.x + sin(4*pi*a.c) + a.b.x/Q' |
---|
328 | |
---|
329 | |
---|
330 | # Check dependency builder |
---|
331 | # Fake parameter class |
---|
332 | class Parameter: |
---|
333 | def __init__(self, name, value=0, expression=''): |
---|
334 | self.path = name |
---|
335 | self.value = value |
---|
336 | self.expression = expression |
---|
337 | def iscomputed(self): return (self.expression != '') |
---|
338 | def __repr__(self): return self.path |
---|
339 | def world(*pars): |
---|
340 | symtab = dict((p.path,p) for p in pars) |
---|
341 | exprs = dict((p.path,p.expression) for p in pars if p.iscomputed()) |
---|
342 | return symtab, exprs |
---|
343 | p1 = Parameter('G0.sigma',5) |
---|
344 | p2 = Parameter('other',expression='2*pi*sin(G0.sigma/.1875) + M1.G1') |
---|
345 | p3 = Parameter('M1.G1',6) |
---|
346 | p4 = Parameter('constant',expression='2*pi*35') |
---|
347 | # Simple chain |
---|
348 | assert set(_find_dependencies(*world(p1,p2,p3))) == set([(p2.path,p1),(p2.path,p3)]) |
---|
349 | # Constant expression |
---|
350 | assert set(_find_dependencies(*world(p1,p4))) == set([(p4.path,None)]) |
---|
351 | # No dependencies |
---|
352 | assert set(_find_dependencies(*world(p1,p3))) == set([]) |
---|
353 | |
---|
354 | # Check function builder |
---|
355 | fn = compile_constraints(*world(p1,p2,p3)) |
---|
356 | |
---|
357 | # Inspect the resulting function |
---|
358 | if 0: |
---|
359 | print(inspect.getdoc(fn)) |
---|
360 | print(dis.dis(fn)) |
---|
361 | |
---|
362 | # Evaluate the function and see if it updates the |
---|
363 | # target value as expected |
---|
364 | fn() |
---|
365 | expected = 2*math.pi*math.sin(5/.1875) + 6 |
---|
366 | assert p2.value == expected,"Value was %s, not %s"%(p2.value,expected) |
---|
367 | |
---|
368 | # Check empty dependency set doesn't crash |
---|
369 | fn = compile_constraints(*world(p1,p3)) |
---|
370 | fn() |
---|
371 | |
---|
372 | # Check that constants are evaluated properly |
---|
373 | fn = compile_constraints(*world(p4)) |
---|
374 | fn() |
---|
375 | assert p4.value == 2*math.pi*35 |
---|
376 | |
---|
377 | # Check additional context example; this also tests multiple |
---|
378 | # expressions |
---|
379 | class Table: |
---|
380 | Si = 2.09 |
---|
381 | values = {'Si': 2.07} |
---|
382 | tbl = Table() |
---|
383 | p5 = Parameter('lookup',expression="tbl.Si") |
---|
384 | fn = compile_constraints(*world(p1,p2,p3,p5),context=dict(tbl=tbl)) |
---|
385 | fn() |
---|
386 | assert p5.value == 2.09,"Value for %s was %s"%(p5.expression,p5.value) |
---|
387 | p5.expression = "tbl.values['Si']" |
---|
388 | fn = compile_constraints(*world(p1,p2,p3,p5),context=dict(tbl=tbl)) |
---|
389 | fn() |
---|
390 | assert p5.value == 2.07,"Value for %s was %s"%(p5.expression,p5.value) |
---|
391 | |
---|
392 | |
---|
393 | # Verify that we capture invalid expressions |
---|
394 | for expr in ['G4.cage', 'M0.cage', 'M1.G1 + *2', |
---|
395 | 'piddle', |
---|
396 | '5; import sys; print "p0wned"', |
---|
397 | '__import__("sys").argv']: |
---|
398 | try: |
---|
399 | p6 = Parameter('broken',expression=expr) |
---|
400 | fn = compile_constraints(*world(p6)) |
---|
401 | fn() |
---|
402 | except Exception as msg: |
---|
403 | #print(msg) |
---|
404 | pass |
---|
405 | else: |
---|
406 | raise "Failed to raise error for %s"%expr |
---|
407 | |
---|
408 | if __name__ == "__main__": |
---|
409 | test_expr() |
---|
410 | test_deps() |
---|