1 | # -*- coding: utf-8 -*- |
---|
2 | ########### SVN repository information ################### |
---|
3 | # $Date: 2015-04-28 19:49:55 +0000 (Tue, 28 Apr 2015) $ |
---|
4 | # $Author: toby $ |
---|
5 | # $Revision: 1818 $ |
---|
6 | # $URL: trunk/GSASIImapvars.py $ |
---|
7 | # $Id: GSASIImapvars.py 1818 2015-04-28 19:49:55Z toby $ |
---|
8 | ########### SVN repository information ################### |
---|
9 | """ |
---|
10 | *GSASIImapvars: Parameter constraints* |
---|
11 | ====================================== |
---|
12 | |
---|
13 | Module to implements algebraic contraints, parameter redefinition |
---|
14 | and parameter simplification contraints. |
---|
15 | |
---|
16 | Parameter redefinition (new vars) is done by creating one or more relationships |
---|
17 | between a set of parameters |
---|
18 | |
---|
19 | :: |
---|
20 | |
---|
21 | Mx1 * Px + My1 * Py +... |
---|
22 | Mx2 * Px + Mz2 * Pz + ... |
---|
23 | |
---|
24 | where Pj is a parameter name and Mjk is a constant. |
---|
25 | |
---|
26 | Constant constraint Relations can also be supplied in the form of an equation: |
---|
27 | |
---|
28 | :: |
---|
29 | |
---|
30 | nx1 * Px + ny1 * Py +... = C1 |
---|
31 | |
---|
32 | where Cn is a constant. These equations define an algebraic |
---|
33 | constrant. |
---|
34 | |
---|
35 | Parameters can also be "fixed" (held), which prevents them from being refined. |
---|
36 | |
---|
37 | All of the above three cases are input using routines |
---|
38 | GroupConstraints and GenerateConstraints. The input consists of a list of |
---|
39 | relationship dictionaries: |
---|
40 | |
---|
41 | .. code-block:: python |
---|
42 | |
---|
43 | constrDict = [ |
---|
44 | {'0:12:Scale': 2.0, '0:14:Scale': 4.0, '0:13:Scale': 3.0, '0:0:Scale': 0.5}, |
---|
45 | {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0}, |
---|
46 | {'0::A0': 0.0}] |
---|
47 | fixedList = ['5.0', None, '0'] |
---|
48 | |
---|
49 | Where the dictionary defines the first part of an expression and the corresponding fixedList |
---|
50 | item is either None (for parameter redefinition) or the constant value, for a constant |
---|
51 | constraint equation. A dictionary that contains a single term defines a variable that |
---|
52 | will be fixed (held). The multiplier and the fixedList value in this case are ignored. |
---|
53 | |
---|
54 | Parameters can also be equivalenced or "slaved" to another parameter, such that one |
---|
55 | (independent) parameter is equated to several (now dependent) parameters. In |
---|
56 | algebraic form this is: |
---|
57 | |
---|
58 | :: |
---|
59 | |
---|
60 | P0 = M1 * P1 = M2 * P2 = ... |
---|
61 | |
---|
62 | Thus parameters P0, P1 and P2,... are linearly equivalent. Routine StoreEquivalence is |
---|
63 | used to specify these equivalences. |
---|
64 | |
---|
65 | Parameter redefinition (new vars) describes a new, independent, parameter, which |
---|
66 | is defined in terms of dependent parameters that are defined in the |
---|
67 | Model, while fixed constrained relations effectively reduce the complexity |
---|
68 | of the Model by removing a degree of freedom. It is possible for a parameter to appear |
---|
69 | in both a parameter redefinition expression and a fixed constraint equation, but a |
---|
70 | parameter cannot be used a parameter equivalance cannot be used elsewhere (not fixed, |
---|
71 | constrained or redefined). Likewise a fixed parameter cannot be used elsewhere (not |
---|
72 | equivalanced, constrained or redefined). |
---|
73 | |
---|
74 | Relationships are grouped so that a set of dependent parameters appear |
---|
75 | in only one group (done in routine GroupConstraints.) Note that if a |
---|
76 | group contains relationships/equations that involve N dependent |
---|
77 | parameters, there must exist N-C new parameters, where C is the number |
---|
78 | of contraint equations in the group. Routine GenerateConstraints takes |
---|
79 | the output from GroupConstraints and generates the |
---|
80 | "missing" relationships and saves that information in the module's |
---|
81 | global variables. Each generated parameter is named sequentially using paramPrefix. |
---|
82 | |
---|
83 | A list of parameters that will be varied is specified as input to GenerateConstraints |
---|
84 | (varyList). A fixed parameter will simply be removed from this list preventing that |
---|
85 | parameter from being varied. Note that all parameters in a constraint relationship |
---|
86 | must specified as varied (appear in varyList) or none can be varied. This is |
---|
87 | checked in GenerateConstraints. Likewise, if all parameters in a constraint are |
---|
88 | not referenced in a refinement, the constraint is ignored, but if some parameters |
---|
89 | in a constraint group are not referenced in a refinement, but others are this |
---|
90 | constitutes and error. |
---|
91 | |
---|
92 | * When a new variable is created, the variable is assigned the name associated |
---|
93 | in the constraint definition or it is assigned a default name of form |
---|
94 | ``::constr<n>`` (see paramPrefix). The vary setting for variables used in the |
---|
95 | constraint are ignored. |
---|
96 | Note that any generated "missing" relations are not varied. Only |
---|
97 | the input relations can be are varied. |
---|
98 | |
---|
99 | * If all parameters in a fixed constraint equation are varied, the generated "missing" |
---|
100 | relations in the group are all varied. This provides the N-C degrees of freedom. |
---|
101 | |
---|
102 | *External Routines* |
---|
103 | ------------------- |
---|
104 | |
---|
105 | To define a set of constrained and unconstrained relations, one |
---|
106 | defines a list of dictionary defining constraint parameters and their |
---|
107 | values, a list of fixed values for each constraint and a list of |
---|
108 | parameters to be varied. In addition, one uses |
---|
109 | :func:`StoreEquivalence` to define parameters that are equivalent. One |
---|
110 | can then use :func:`CheckConstraints` to check that the input is |
---|
111 | internally consistent and finally :func:`GroupConstraints` and |
---|
112 | :func:`GenerateConstraints` to generate the internally used |
---|
113 | tables. Routines :func:`Map2Dict` is used to initialize the parameter |
---|
114 | dictionary and :func:`Dict2Map`, :func:`Dict2Deriv`, and |
---|
115 | :func:`ComputeDepESD` are used to apply constraints. Routine |
---|
116 | :func:`VarRemapShow` is used to print out the constraint information, |
---|
117 | as stored by :func:`GenerateConstraints`. |
---|
118 | |
---|
119 | :func:`InitVars` |
---|
120 | This is optionally used to clear out all defined previously defined constraint information |
---|
121 | |
---|
122 | :func:`StoreEquivalence` |
---|
123 | To implement parameter redefinition, one calls StoreEquivalence. This should be called for every set of |
---|
124 | equivalence relationships. There is no harm in using StoreEquivalence with the same independent variable: |
---|
125 | |
---|
126 | .. code-block:: python |
---|
127 | |
---|
128 | StoreEquivalence('x',('y',)) |
---|
129 | StoreEquivalence('x',('z',)) |
---|
130 | |
---|
131 | or equivalently |
---|
132 | |
---|
133 | .. code-block:: python |
---|
134 | |
---|
135 | StoreEquivalence('x',('y','z')) |
---|
136 | |
---|
137 | The latter will run more efficiently. Note that mixing independent and dependent variables is |
---|
138 | problematic. This is not allowed: |
---|
139 | |
---|
140 | .. code-block:: python |
---|
141 | |
---|
142 | StoreEquivalence('x',('y',)) |
---|
143 | StoreEquivalence('y',('z',)) |
---|
144 | |
---|
145 | Use StoreEquivalence before calling GenerateConstraints or CheckConstraints |
---|
146 | |
---|
147 | :func:`CheckConstraints` |
---|
148 | To check that input in internally consistent, use CheckConstraints |
---|
149 | |
---|
150 | :func:`Map2Dict` |
---|
151 | To determine values for the parameters created in this module, one |
---|
152 | calls Map2Dict. This will not apply contraints. |
---|
153 | |
---|
154 | :func:`Dict2Map` |
---|
155 | To take values from the new independent parameters and constraints, |
---|
156 | one calls Dict2Map. This will apply contraints. |
---|
157 | |
---|
158 | :func:`Dict2Deriv` |
---|
159 | Use Dict2Deriv to determine derivatives on independent parameters |
---|
160 | from those on dependent ones |
---|
161 | |
---|
162 | :func:`ComputeDepESD` |
---|
163 | Use ComputeDepESD to compute uncertainties on dependent variables |
---|
164 | |
---|
165 | :func:`VarRemapShow` |
---|
166 | To show a summary of the parameter remapping, one calls VarRemapShow |
---|
167 | |
---|
168 | *Global Variables* |
---|
169 | ------------------ |
---|
170 | |
---|
171 | dependentParmList: |
---|
172 | contains a list by group of lists of |
---|
173 | parameters used in the group. Note that parameters listed in |
---|
174 | dependentParmList should not be refined as they will not affect |
---|
175 | the model |
---|
176 | |
---|
177 | indParmList: |
---|
178 | a list of groups of Independent parameters defined in |
---|
179 | each group. This contains both parameters used in parameter |
---|
180 | redefinitions as well as names of generated new parameters. |
---|
181 | |
---|
182 | fixedVarList: |
---|
183 | a list of variables that have been 'fixed' |
---|
184 | by defining them as equal to a constant (::var: = 0). Note that |
---|
185 | the constant value is ignored at present. These variables are |
---|
186 | later removed from varyList which prevents them from being refined. |
---|
187 | Unlikely to be used externally. |
---|
188 | |
---|
189 | arrayList: |
---|
190 | a list by group of relationship matrices to relate |
---|
191 | parameters in dependentParmList to those in indParmList. Unlikely |
---|
192 | to be used externally. |
---|
193 | |
---|
194 | invarrayList: |
---|
195 | a list by group of relationship matrices to relate |
---|
196 | parameters in indParmList to those in dependentParmList. Unlikely |
---|
197 | to be used externally. |
---|
198 | |
---|
199 | fixedDict: |
---|
200 | a dictionary containing the fixed values corresponding |
---|
201 | to parameter equations. The dict key is an ascii string, but the |
---|
202 | dict value is a float. Unlikely to be used externally. |
---|
203 | |
---|
204 | *Routines* |
---|
205 | ---------- |
---|
206 | |
---|
207 | Note that parameter names in GSAS-II are strings of form ``<ph>:<hst>:<nam>`` |
---|
208 | |
---|
209 | """ |
---|
210 | |
---|
211 | import numpy as np |
---|
212 | import GSASIIpath |
---|
213 | GSASIIpath.SetVersionNumber("$Revision: 1818 $") |
---|
214 | # data used for constraints; |
---|
215 | debug = False # turns on printing as constraint input is processed |
---|
216 | # note that constraints are stored listed by contraint groups, where each constraint |
---|
217 | # group contains those parameters that must be handled together |
---|
218 | dependentParmList = [] # contains a list of parameters in each group |
---|
219 | # note that parameters listed in dependentParmList should not be refined |
---|
220 | arrayList = [] # a list of of relationship matrices |
---|
221 | invarrayList = [] # a list of inverse relationship matrices |
---|
222 | indParmList = [] # a list of names for the new parameters |
---|
223 | fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations |
---|
224 | # key is original ascii string, value is float |
---|
225 | fixedVarList = [] # List of variables that should not be refined |
---|
226 | |
---|
227 | # prefix for parameter names |
---|
228 | paramPrefix = "::constr" |
---|
229 | consNum = 0 # number of the next constraint to be created |
---|
230 | |
---|
231 | def InitVars(): |
---|
232 | '''Initializes all constraint information''' |
---|
233 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,consNum |
---|
234 | dependentParmList = [] # contains a list of parameters in each group |
---|
235 | arrayList = [] # a list of of relationship matrices |
---|
236 | invarrayList = [] # a list of inverse relationship matrices |
---|
237 | indParmList = [] # a list of names for the new parameters |
---|
238 | fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations |
---|
239 | consNum = 0 # number of the next constraint to be created |
---|
240 | fixedVarList = [] |
---|
241 | |
---|
242 | def VarKeys(constr): |
---|
243 | """Finds the keys in a constraint that represent variables |
---|
244 | e.g. eliminates any that start with '_' |
---|
245 | |
---|
246 | :param dict constr: a single constraint entry of form:: |
---|
247 | |
---|
248 | {'var1': mult1, 'var2': mult2,... '_notVar': val,...} |
---|
249 | |
---|
250 | (see :func:`GroupConstraints`) |
---|
251 | :returns: a list of keys where any keys beginning with '_' are |
---|
252 | removed. |
---|
253 | """ |
---|
254 | return [i for i in constr.keys() if not i.startswith('_')] |
---|
255 | |
---|
256 | |
---|
257 | def GroupConstraints(constrDict): |
---|
258 | """divide the constraints into groups that share no parameters. |
---|
259 | |
---|
260 | :param dict constrDict: a list of dicts defining relationships/constraints |
---|
261 | |
---|
262 | :: |
---|
263 | |
---|
264 | constrDict = [{<constr1>}, {<constr2>}, ...] |
---|
265 | |
---|
266 | where {<constr1>} is {'var1': mult1, 'var2': mult2,... } |
---|
267 | |
---|
268 | :returns: two lists of lists: |
---|
269 | |
---|
270 | * a list of grouped contraints where each constraint grouped containts a list |
---|
271 | of indices for constraint constrDict entries |
---|
272 | * a list containing lists of parameter names contained in each group |
---|
273 | |
---|
274 | """ |
---|
275 | assignedlist = [] # relationships that have been used |
---|
276 | groups = [] # contains a list of grouplists |
---|
277 | ParmList = [] |
---|
278 | for i,consi in enumerate(constrDict): |
---|
279 | if i in assignedlist: continue # already in a group, skip |
---|
280 | # starting a new group |
---|
281 | grouplist = [i,] |
---|
282 | assignedlist.append(i) |
---|
283 | groupset = set(VarKeys(consi)) |
---|
284 | changes = True # always loop at least once |
---|
285 | while(changes): # loop until we can't find anything to add to the current group |
---|
286 | changes = False # but don't loop again unless we find something |
---|
287 | for j,consj in enumerate(constrDict): |
---|
288 | if j in assignedlist: continue # already in a group, skip |
---|
289 | if len(set(VarKeys(consj)) & groupset) > 0: # true if this needs to be added |
---|
290 | changes = True |
---|
291 | grouplist.append(j) |
---|
292 | assignedlist.append(j) |
---|
293 | groupset = groupset | set(VarKeys(consj)) |
---|
294 | group = sorted(grouplist) |
---|
295 | varlist = sorted(list(groupset)) |
---|
296 | groups.append(group) |
---|
297 | ParmList.append(varlist) |
---|
298 | return groups,ParmList |
---|
299 | |
---|
300 | def CheckConstraints(varyList,constrDict,fixedList): |
---|
301 | '''Takes a list of relationship entries comprising a group of |
---|
302 | constraints and checks for inconsistencies such as conflicts in |
---|
303 | parameter/variable definitions and or inconsistently varied parameters. |
---|
304 | |
---|
305 | :param list varyList: a list of parameters names that will be varied |
---|
306 | |
---|
307 | :param dict constrDict: a list of dicts defining relationships/constraints |
---|
308 | (as created in :func:`GSASIIstrIO.ProcessConstraints` and |
---|
309 | documented in :func:`GroupConstraints`) |
---|
310 | |
---|
311 | :param list fixedList: a list of values specifying a fixed value for each |
---|
312 | dict in constrDict. Values are either strings that can be converted to |
---|
313 | floats or ``None`` if the constraint defines a new parameter rather |
---|
314 | than a constant. |
---|
315 | |
---|
316 | :returns: two strings: |
---|
317 | |
---|
318 | * the first lists conflicts internal to the specified constraints |
---|
319 | * the second lists conflicts where the varyList specifies some |
---|
320 | parameters in a constraint, but not all |
---|
321 | |
---|
322 | If there are no errors, both strings will be empty |
---|
323 | ''' |
---|
324 | import re |
---|
325 | global dependentParmList,arrayList,invarrayList,indParmList,consNum |
---|
326 | errmsg = '' |
---|
327 | warnmsg = '' |
---|
328 | fixVlist = [] |
---|
329 | # process fixed variables (holds) |
---|
330 | for cdict in constrDict: |
---|
331 | # N.B. No "_" names in holds |
---|
332 | if len(cdict) == 1: |
---|
333 | fixVlist.append(cdict.keys()[0]) |
---|
334 | |
---|
335 | # process equivalences: make a list of dependent and independent vars |
---|
336 | # and check for repeated uses (repetition of a parameter as an |
---|
337 | # independent var is OK) |
---|
338 | indepVarList = [] |
---|
339 | depVarList = [] |
---|
340 | multdepVarList = [] |
---|
341 | for varlist,mapvars,multarr,invmultarr in zip( |
---|
342 | dependentParmList,indParmList,arrayList,invarrayList): |
---|
343 | if multarr is None: # an equivalence |
---|
344 | zeromult = False |
---|
345 | for mv in mapvars: |
---|
346 | varied = 0 |
---|
347 | notvaried = '' |
---|
348 | if mv in varyList: |
---|
349 | varied += 1 |
---|
350 | else: |
---|
351 | if notvaried: notvaried += ', ' |
---|
352 | notvaried += mv |
---|
353 | if mv not in indepVarList: indepVarList.append(mv) |
---|
354 | for v,m in zip(varlist,invmultarr): |
---|
355 | if v in indepVarList: |
---|
356 | errmsg += '\nVariable '+v+' is used to set values in a constraint before its value is set in another constraint\n' |
---|
357 | if m == 0: zeromult = True |
---|
358 | if v in varyList: |
---|
359 | varied += 1 |
---|
360 | else: |
---|
361 | if notvaried: notvaried += ', ' |
---|
362 | notvaried += v |
---|
363 | if v in depVarList: |
---|
364 | multdepVarList.append(v) |
---|
365 | else: |
---|
366 | depVarList.append(v) |
---|
367 | if varied > 0 and varied != len(varlist)+1: |
---|
368 | warnmsg += "\nNot all variables refined in equivalence:\n\t" |
---|
369 | s = "" |
---|
370 | for v in varlist: |
---|
371 | if s != "": s+= " & " |
---|
372 | s += str(v) |
---|
373 | warnmsg += str(mv) + " => " + s |
---|
374 | warnmsg += '\nNot refined: ' + notvaried + '\n' |
---|
375 | if zeromult: |
---|
376 | errmsg += "\nZero multiplier is invalid in equivalence:\n\t" |
---|
377 | s = "" |
---|
378 | for v in varlist: |
---|
379 | if s != "": s+= " & " |
---|
380 | s += str(v) |
---|
381 | errmsg += str(mv) + " => " + s + '\n' |
---|
382 | |
---|
383 | # check for errors: |
---|
384 | if len(multdepVarList) > 0: |
---|
385 | errmsg += "\nThe following parameters(s) are used in conflicting Equivalence relations as dependent variables:\n" |
---|
386 | s = '' |
---|
387 | for var in sorted(set(multdepVarList)): |
---|
388 | if s != "": s+= ", " |
---|
389 | s += str(var) |
---|
390 | errmsg += '\t'+ s + '\n' |
---|
391 | equivVarList = list(set(indepVarList).union(set(depVarList))) |
---|
392 | if debug: print 'equivVarList',equivVarList |
---|
393 | inboth = set(fixVlist).intersection(set(equivVarList)) |
---|
394 | if len(inboth) > 0: |
---|
395 | errmsg += "\nThe following parameter(s) are used in both Equivalence and Fixed constraints:\n" |
---|
396 | s = '' |
---|
397 | for var in sorted(inboth): |
---|
398 | if s != "": s+= ", " |
---|
399 | s += str(var) |
---|
400 | errmsg += '\t'+ s + '\n' |
---|
401 | |
---|
402 | groups,parmlist = GroupConstraints(constrDict) |
---|
403 | # scan through parameters in each relationship. Are all varied? If only some are |
---|
404 | # varied, create a warning message. |
---|
405 | for group,varlist in zip(groups,parmlist): |
---|
406 | if len(varlist) == 1: continue |
---|
407 | for rel in group: |
---|
408 | varied = 0 |
---|
409 | notvaried = '' |
---|
410 | for var in constrDict[rel]: |
---|
411 | if var.startswith('_'): continue |
---|
412 | if not re.match('[0-9]*:[0-9\*]*:',var): |
---|
413 | warnmsg += "\nVariable "+str(var)+" does not begin with a ':'" |
---|
414 | if var in varyList: |
---|
415 | varied += 1 |
---|
416 | else: |
---|
417 | if notvaried: notvaried += ', ' |
---|
418 | notvaried += var |
---|
419 | if var in fixVlist: |
---|
420 | errmsg += '\nParameter '+var+" is Fixed and used in a constraint:\n\t" |
---|
421 | errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])+"\n" |
---|
422 | if varied > 0 and varied != len(VarKeys(constrDict[rel])): |
---|
423 | warnmsg += "\nNot all variables refined in constraint:\n\t" |
---|
424 | warnmsg += _FormatConstraint(constrDict[rel],fixedList[rel]) |
---|
425 | warnmsg += '\nNot refined: ' + notvaried + '\n' |
---|
426 | if errmsg or warnmsg: |
---|
427 | return errmsg,warnmsg |
---|
428 | |
---|
429 | # now look for process each group and create the relations that are needed to form |
---|
430 | # non-singular square matrix |
---|
431 | for group,varlist in zip(groups,parmlist): |
---|
432 | if len(varlist) == 1: continue # a constraint group with a single variable can be ignored |
---|
433 | if len(varlist) < len(group): # too many relationships -- no can do |
---|
434 | errmsg += "\nOver-constrained input. " |
---|
435 | errmsg += "There are more constraints " + str(len(group)) |
---|
436 | errmsg += "\n\tthan variables " + str(len(varlist)) + "\n" |
---|
437 | for rel in group: |
---|
438 | errmsg += _FormatConstraint(constrDict[rel],fixedList[rel]) |
---|
439 | errmsg += "\n" |
---|
440 | continue |
---|
441 | try: |
---|
442 | multarr = _FillArray(group,constrDict,varlist) |
---|
443 | _RowEchelon(len(group),multarr,varlist) |
---|
444 | except: |
---|
445 | errmsg += "\nSingular input. " |
---|
446 | errmsg += "There are internal inconsistencies in these constraints\n" |
---|
447 | for rel in group: |
---|
448 | errmsg += _FormatConstraint(constrDict[rel],fixedList[rel]) |
---|
449 | errmsg += "\n" |
---|
450 | continue |
---|
451 | try: |
---|
452 | multarr = _FillArray(group,constrDict,varlist,FillDiagonals=True) |
---|
453 | GramSchmidtOrtho(multarr,len(group)) |
---|
454 | except: |
---|
455 | errmsg += "\nUnexpected singularity with constraints (in Gram-Schmidt)\n" |
---|
456 | for rel in group: |
---|
457 | errmsg += _FormatConstraint(constrDict[rel],fixedList[rel]) |
---|
458 | errmsg += "\n" |
---|
459 | continue |
---|
460 | mapvar = [] |
---|
461 | group = group[:] |
---|
462 | # scan through all generated and input variables |
---|
463 | # Check again for inconsistent variable use |
---|
464 | # for new variables -- where varied and unvaried parameters get grouped |
---|
465 | # together. I don't think this can happen when not flagged before, but |
---|
466 | # it does not hurt to check again. |
---|
467 | for i in range(len(varlist)): |
---|
468 | varied = 0 |
---|
469 | notvaried = '' |
---|
470 | if len(group) > 0: |
---|
471 | rel = group.pop(0) |
---|
472 | fixedval = fixedList[rel] |
---|
473 | for var in VarKeys(constrDict[rel]): |
---|
474 | if var in varyList: |
---|
475 | varied += 1 |
---|
476 | else: |
---|
477 | if notvaried: notvaried += ', ' |
---|
478 | notvaried += var |
---|
479 | else: |
---|
480 | fixedval = None |
---|
481 | if fixedval is None: |
---|
482 | varname = paramPrefix + str(consNum) # assign a name to a variable |
---|
483 | mapvar.append(varname) |
---|
484 | consNum += 1 |
---|
485 | else: |
---|
486 | mapvar.append(fixedval) |
---|
487 | if varied > 0 and notvaried != '': |
---|
488 | warnmsg += "\nNot all variables refined in generated constraint" |
---|
489 | warnmsg += '\nPlease report this unexpected error\n' |
---|
490 | for rel in group: |
---|
491 | warnmsg += _FormatConstraint(constrDict[rel],fixedList[rel]) |
---|
492 | warnmsg += "\n" |
---|
493 | warnmsg += '\n\tNot refined: ' + notvaried + '\n' |
---|
494 | try: |
---|
495 | np.linalg.inv(multarr) |
---|
496 | except: |
---|
497 | errmsg += "\nSingular input. " |
---|
498 | errmsg += "The following constraints are not " |
---|
499 | errmsg += "linearly independent\n\tor do not " |
---|
500 | errmsg += "allow for generation of a non-singular set\n" |
---|
501 | errmsg += 'Please report this unexpected error\n' |
---|
502 | for rel in group: |
---|
503 | errmsg += _FormatConstraint(constrDict[rel],fixedList[rel]) |
---|
504 | errmsg += "\n" |
---|
505 | return errmsg,warnmsg |
---|
506 | |
---|
507 | def GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict=None,SeqHist=None): |
---|
508 | '''Takes a list of relationship entries comprising a group of |
---|
509 | constraints and builds the relationship lists and their inverse |
---|
510 | and stores them in global variables Also checks for internal |
---|
511 | conflicts or inconsistencies in parameter/variable definitions. |
---|
512 | |
---|
513 | :param list groups: a list of grouped contraints where each constraint |
---|
514 | grouped containts a list of indices for constraint constrDict entries, |
---|
515 | created in :func:`GroupConstraints` (returned as 1st value) |
---|
516 | |
---|
517 | :param list parmlist: a list containing lists of parameter names |
---|
518 | contained in each group, created in :func:`GroupConstraints` |
---|
519 | (returned as 2nd value) |
---|
520 | |
---|
521 | :param list varyList: a list of parameters names (strings of form |
---|
522 | ``<ph>:<hst>:<nam>``) that will be varied. Note that this is changed here. |
---|
523 | |
---|
524 | :param dict constrDict: a list of dicts defining relationships/constraints |
---|
525 | (as defined in :func:`GroupConstraints`) |
---|
526 | |
---|
527 | :param list fixedList: a list of values specifying a fixed value for each |
---|
528 | dict in constrDict. Values are either strings that can be converted to |
---|
529 | floats, float values or None if the constraint defines a new parameter. |
---|
530 | |
---|
531 | :param dict parmDict: a dict containing all parameters defined in current |
---|
532 | refinement. |
---|
533 | |
---|
534 | :param int SeqHist: number of current histogram, when used in a sequential |
---|
535 | refinement. None (default) otherwise. Wildcard variable names are |
---|
536 | set to the current histogram, when found if not None. |
---|
537 | ''' |
---|
538 | global dependentParmList,arrayList,invarrayList,indParmList,consNum |
---|
539 | msg = '' |
---|
540 | |
---|
541 | # process fixed (held) variables |
---|
542 | for cdict in constrDict: |
---|
543 | if len(cdict) == 1: |
---|
544 | fixedVarList.append(cdict.keys()[0]) |
---|
545 | |
---|
546 | # process equivalences: make a list of dependent and independent vars |
---|
547 | # and check for repeated uses (repetition of a parameter as an |
---|
548 | # independent var is OK [A=B; A=C], but chaining: [A=B; B=C] is not good) |
---|
549 | indepVarList = [] |
---|
550 | depVarList = [] |
---|
551 | multdepVarList = [] |
---|
552 | translateTable = {} # lookup table for wildcard referenced variables |
---|
553 | for varlist,mapvars,multarr,invmultarr in zip( # process equivalences |
---|
554 | dependentParmList,indParmList,arrayList,invarrayList): |
---|
555 | if multarr is None: # true only if an equivalence |
---|
556 | zeromult = False |
---|
557 | for mv in mapvars: |
---|
558 | #s = '' |
---|
559 | varied = 0 |
---|
560 | notvaried = '' |
---|
561 | if mv in varyList: |
---|
562 | varied += 1 |
---|
563 | else: |
---|
564 | if notvaried: notvaried += ', ' |
---|
565 | notvaried += mv |
---|
566 | if parmDict is not None and mv not in parmDict: |
---|
567 | print "Dropping equivalence for variable "+str(mv)+". Not defined in this refinement" |
---|
568 | #msg += "\nCannot equivalence to variable "+str(mv)+". Not defined in this refinement" |
---|
569 | #continue |
---|
570 | else: |
---|
571 | if mv not in indepVarList: indepVarList.append(mv) |
---|
572 | for v,m in zip(varlist,invmultarr): |
---|
573 | if parmDict is not None and v not in parmDict: |
---|
574 | print "Dropping equivalence for dep. variable "+str(v)+". Not defined in this refinement" |
---|
575 | continue |
---|
576 | if m == 0: zeromult = True |
---|
577 | if v in varyList: |
---|
578 | varied += 1 |
---|
579 | else: |
---|
580 | if notvaried: notvaried += ', ' |
---|
581 | notvaried += v |
---|
582 | if v in depVarList: |
---|
583 | multdepVarList.append(v) |
---|
584 | else: |
---|
585 | depVarList.append(v) |
---|
586 | if varied > 0 and varied != len(varlist)+1: |
---|
587 | msg += "\nNot all variables refined in equivalence:\n\t" |
---|
588 | s = "" |
---|
589 | for v in varlist: |
---|
590 | if s != "": s+= " & " |
---|
591 | s += str(v) |
---|
592 | msg += str(mv) + " => " + s |
---|
593 | msg += '\nNot refined: ' + notvaried + '\n' |
---|
594 | if zeromult: |
---|
595 | msg += "\nZero multiplier is invalid in equivalence:\n\t" |
---|
596 | s = "" |
---|
597 | for v in varlist: |
---|
598 | if s != "": s+= " & " |
---|
599 | s += str(v) |
---|
600 | msg += str(mv) + " => " + s + '\n' |
---|
601 | # save the lists of dep. and indep. vars (after dropping unused) |
---|
602 | global independentVars |
---|
603 | independentVars = indepVarList |
---|
604 | #print 'independentVars=',independentVars |
---|
605 | equivVarList = list(set(indepVarList).union(set(depVarList))) |
---|
606 | |
---|
607 | # scan through parameters in each relationship. Are all varied? If only some are |
---|
608 | # varied, create an error message. |
---|
609 | for group,varlist in zip(groups,parmlist): |
---|
610 | if len(varlist) == 1: continue |
---|
611 | for rel in group: |
---|
612 | varied = 0 |
---|
613 | notvaried = '' |
---|
614 | unused = 0 |
---|
615 | notused = '' |
---|
616 | for var in constrDict[rel]: |
---|
617 | if var.startswith('_'): continue |
---|
618 | if var.split(':')[1] == '*' and SeqHist is not None: |
---|
619 | # convert wildcard var to reference current histogram; save translation in table |
---|
620 | sv = var.split(':') |
---|
621 | sv[1] = str(SeqHist) |
---|
622 | translateTable[var] = ':'.join(sv) |
---|
623 | var = translateTable[var] |
---|
624 | if parmDict is not None and var not in parmDict: |
---|
625 | unused += 1 |
---|
626 | if notvaried: notused += ', ' |
---|
627 | notused += var |
---|
628 | if var in varyList: |
---|
629 | varied += 1 |
---|
630 | else: |
---|
631 | if notvaried: notvaried += ', ' |
---|
632 | notvaried += var |
---|
633 | if var in fixedVarList: |
---|
634 | msg += '\nError: parameter '+var+" is Fixed and used in a constraint:\n\t" |
---|
635 | msg += _FormatConstraint(constrDict[rel],fixedList[rel])+"\n" |
---|
636 | #if unused > 0:# and unused != len(VarKeys(constrDict[rel])): |
---|
637 | if unused > 0 and unused != len(VarKeys(constrDict[rel])): |
---|
638 | msg += "\nSome (but not all) variables in constraint are not defined:\n\t" |
---|
639 | msg += _FormatConstraint(constrDict[rel],fixedList[rel]) |
---|
640 | msg += '\nNot used: ' + notused + '\n' |
---|
641 | if varied > 0 and varied != len(VarKeys(constrDict[rel])): |
---|
642 | msg += "\nNot all variables refined in constraint:\n\t" |
---|
643 | msg += _FormatConstraint(constrDict[rel],fixedList[rel]) |
---|
644 | msg += '\nNot refined: ' + notvaried + '\n' |
---|
645 | # if there were errors found, go no farther |
---|
646 | if msg: |
---|
647 | print ' *** ERROR in constraint definitions! ***' |
---|
648 | print msg |
---|
649 | raise Exception |
---|
650 | |
---|
651 | # now process each group and create the relations that are needed to form |
---|
652 | # a non-singular square matrix |
---|
653 | # If all are varied and this is a constraint equation, then set VaryFree flag |
---|
654 | # so that the newly created relationships will be varied |
---|
655 | for group,varlist in zip(groups,parmlist): |
---|
656 | if len(varlist) == 1: continue |
---|
657 | # for constraints, if all included variables are refined, |
---|
658 | # set the VaryFree flag, and remaining degrees of freedom will be |
---|
659 | # varied (since consistency was checked, if any one variable is |
---|
660 | # refined, then assume that all are) |
---|
661 | varsList = [] # make a list of all the referenced variables as well |
---|
662 | VaryFree = False |
---|
663 | for rel in group: |
---|
664 | varied = 0 |
---|
665 | unused = 0 |
---|
666 | for var in VarKeys(constrDict[rel]): |
---|
667 | var = translateTable.get(var,var) # replace wildcards |
---|
668 | if parmDict is not None and var not in parmDict: |
---|
669 | unused += 1 |
---|
670 | if var not in varsList: varsList.append(var) |
---|
671 | if var in varyList: varied += 1 |
---|
672 | if fixedList[rel] is not None and varied > 0: |
---|
673 | VaryFree = True |
---|
674 | if len(varlist) < len(group): # too many relationships -- no can do |
---|
675 | msg = 'too many relationships' |
---|
676 | break |
---|
677 | # Since we checked before, if any variables are unused, then all must be. |
---|
678 | # If so, this set of relationships can be ignored |
---|
679 | if unused: |
---|
680 | if debug: print('Constraint ignored (all variables undefined)') |
---|
681 | if debug: print (' '+_FormatConstraint(constrDict[rel],fixedList[rel])) |
---|
682 | continue |
---|
683 | # fill in additional degrees of freedom |
---|
684 | try: |
---|
685 | arr = _FillArray(group,constrDict,varlist) |
---|
686 | _RowEchelon(len(group),arr,varlist) |
---|
687 | constrArr = _FillArray(group,constrDict,varlist,FillDiagonals=True) |
---|
688 | GramSchmidtOrtho(constrArr,len(group)) |
---|
689 | except: |
---|
690 | msg = 'Singular relationships' |
---|
691 | break |
---|
692 | mapvar = [] |
---|
693 | group = group[:] |
---|
694 | # scan through all generated and input relationships, we need to add to the varied list |
---|
695 | # all the new parameters where VaryFree has been set or where a New Var is varied. |
---|
696 | # |
---|
697 | # If a group does not contain any fixed values (constraint equations) |
---|
698 | # and nothing in the group is varied, drop this group, so that the |
---|
699 | # dependent parameters can be refined individually. |
---|
700 | unused = True |
---|
701 | for i in range(len(varlist)): |
---|
702 | if len(group) > 0: # get the original equation reference |
---|
703 | rel = group.pop(0) |
---|
704 | fixedval = fixedList[rel] |
---|
705 | varyflag = constrDict[rel].get('_vary',False) |
---|
706 | varname = constrDict[rel].get('_name','') |
---|
707 | else: # this relationship has been generated |
---|
708 | varyflag = False |
---|
709 | varname = '' |
---|
710 | fixedval = None |
---|
711 | if fixedval is None: # this is a new variable, not a constraint |
---|
712 | if not varname: |
---|
713 | varname = paramPrefix + str(consNum) # no assigned name, create one |
---|
714 | consNum += 1 |
---|
715 | mapvar.append(varname) |
---|
716 | # vary the new relationship if it is a degree of freedom in |
---|
717 | # a set of contraint equations or if a New Var is flagged to be varied. |
---|
718 | if VaryFree or varyflag: |
---|
719 | unused = False |
---|
720 | varyList.append(varname) |
---|
721 | # fix (prevent varying) of all the variables inside the constraint group |
---|
722 | # (dependent vars) |
---|
723 | for var in varsList: |
---|
724 | if var in varyList: varyList.remove(var) |
---|
725 | else: |
---|
726 | unused = False |
---|
727 | mapvar.append(fixedval) |
---|
728 | if unused: # skip over constraints that don't matter (w/o fixed value or any refined variables) |
---|
729 | if debug: print('Constraint ignored (all variables unrefined)') |
---|
730 | if debug: print (' '+_FormatConstraint(constrDict[rel],fixedList[rel])) |
---|
731 | continue |
---|
732 | dependentParmList.append([translateTable.get(var,var) for var in varlist]) |
---|
733 | arrayList.append(constrArr) |
---|
734 | invarrayList.append(np.linalg.inv(constrArr)) |
---|
735 | indParmList.append(mapvar) |
---|
736 | if msg: |
---|
737 | print ' *** ERROR in constraint definitions! ***' |
---|
738 | print msg |
---|
739 | print VarRemapShow(varyList) |
---|
740 | raise Exception |
---|
741 | # setup dictionary containing the fixed values |
---|
742 | global fixedDict |
---|
743 | # key is original ascii string, value is float |
---|
744 | for fixedval in fixedList: |
---|
745 | if fixedval: |
---|
746 | fixedDict[fixedval] = float(fixedval) |
---|
747 | |
---|
748 | # make list of dependent variables |
---|
749 | global dependentVars |
---|
750 | depVarList = [] |
---|
751 | for varlist,mapvars,invmultarr in zip( # process equivalences |
---|
752 | dependentParmList,indParmList,invarrayList): |
---|
753 | for i,mv in enumerate(varlist): |
---|
754 | if mv not in depVarList: depVarList.append(mv) |
---|
755 | dependentVars = depVarList |
---|
756 | if debug: # on debug, show what is parsed & generated, semi-readable |
---|
757 | print 50*'-' |
---|
758 | print VarRemapShow(varyList) |
---|
759 | print 'Varied: ',varyList |
---|
760 | print 'Not Varied: ',fixedVarList |
---|
761 | |
---|
762 | def StoreEquivalence(independentVar,dependentList): |
---|
763 | '''Takes a list of dependent parameter(s) and stores their |
---|
764 | relationship to a single independent parameter (independentVar) |
---|
765 | |
---|
766 | :param str independentVar: name of master parameter that will be used to determine the value |
---|
767 | to set the dependent variables |
---|
768 | |
---|
769 | :param list dependentList: a list of parameters that will set from |
---|
770 | independentVar. Each item in the list can be a string with the parameter |
---|
771 | name or a tuple containing a name and multiplier: |
---|
772 | ``['parm1',('parm2',.5),]`` |
---|
773 | |
---|
774 | ''' |
---|
775 | |
---|
776 | global dependentParmList,arrayList,invarrayList,indParmList |
---|
777 | mapList = [] |
---|
778 | multlist = [] |
---|
779 | for var in dependentList: |
---|
780 | if isinstance(var, basestring): |
---|
781 | mult = 1.0 |
---|
782 | elif len(var) == 2: |
---|
783 | var,mult = var |
---|
784 | else: |
---|
785 | raise Exception("Cannot parse "+repr(var) + " as var or (var,multiplier)") |
---|
786 | mapList.append(var) |
---|
787 | multlist.append(tuple((mult,))) |
---|
788 | # added relationships to stored values |
---|
789 | arrayList.append(None) |
---|
790 | invarrayList.append(np.array(multlist)) |
---|
791 | indParmList.append(tuple((independentVar,))) |
---|
792 | dependentParmList.append(mapList) |
---|
793 | return |
---|
794 | |
---|
795 | def GetDependentVars(): |
---|
796 | '''Return a list of dependent variables: e.g. variables that are |
---|
797 | constrained in terms of other variables |
---|
798 | |
---|
799 | :returns: a list of variable names |
---|
800 | |
---|
801 | ''' |
---|
802 | return dependentVars |
---|
803 | |
---|
804 | def GetIndependentVars(): |
---|
805 | '''Return a list of independent variables: e.g. variables that are |
---|
806 | created by constraints of other variables |
---|
807 | |
---|
808 | :returns: a list of variable names |
---|
809 | |
---|
810 | ''' |
---|
811 | return independentVars |
---|
812 | |
---|
813 | def PrintIndependentVars(parmDict,varyList,sigDict,PrintAll=False,pFile=None): |
---|
814 | '''Print the values and uncertainties on the independent variables''' |
---|
815 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict |
---|
816 | printlist = [] |
---|
817 | mapvars = GetIndependentVars() |
---|
818 | for i,name in enumerate(mapvars): |
---|
819 | if name in fixedDict: continue |
---|
820 | if PrintAll or name in varyList: |
---|
821 | sig = sigDict.get(name) |
---|
822 | printlist.append([name,parmDict[name],sig]) |
---|
823 | if len(printlist) == 0: return |
---|
824 | s1 = '' |
---|
825 | print >>pFile,130*'-' |
---|
826 | print >>pFile,"Variables generated by constraints" |
---|
827 | printlist.append(3*[None]) |
---|
828 | for name,val,esd in printlist: |
---|
829 | if len(s1) > 120 or name is None: |
---|
830 | print >>pFile,'' |
---|
831 | print >>pFile,s1 |
---|
832 | print >>pFile,s2 |
---|
833 | print >>pFile,s3 |
---|
834 | s1 = '' |
---|
835 | if name is None: break |
---|
836 | if s1 == "": |
---|
837 | s1 = ' name :' |
---|
838 | s2 = ' value :' |
---|
839 | s3 = ' sig :' |
---|
840 | s1 += '%15s' % (name) |
---|
841 | s2 += '%15.4f' % (val) |
---|
842 | if esd is None: |
---|
843 | s3 += '%15s' % ('n/a') |
---|
844 | else: |
---|
845 | s3 += '%15.4f' % (esd) |
---|
846 | |
---|
847 | def ComputeDepESD(covMatrix,varyList,parmDict): |
---|
848 | '''Compute uncertainties for dependent parameters from independent ones |
---|
849 | returns a dictionary containing the esd values for dependent parameters |
---|
850 | ''' |
---|
851 | sigmaDict = {} |
---|
852 | for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList): |
---|
853 | #if invmultarr is None: continue # probably not needed |
---|
854 | try: |
---|
855 | valuelist = [parmDict[var] for var in mapvars] |
---|
856 | except KeyError: |
---|
857 | continue |
---|
858 | # get the v-covar matrix for independent parameters |
---|
859 | vcov = np.zeros((len(mapvars),len(mapvars))) |
---|
860 | for i1,name1 in enumerate(mapvars): |
---|
861 | if name1 not in varyList: continue |
---|
862 | iv1 = varyList.index(name1) |
---|
863 | for i2,name2 in enumerate(mapvars): |
---|
864 | if name2 not in varyList: continue |
---|
865 | iv2 = varyList.index(name2) |
---|
866 | vcov[i1][i2] = covMatrix[iv1][iv2] |
---|
867 | # vec is the vector that multiplies each of the independent values |
---|
868 | for v,vec in zip(varlist,invmultarr): |
---|
869 | sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec))) |
---|
870 | return sigmaDict |
---|
871 | |
---|
872 | def _FormatConstraint(RelDict,RelVal): |
---|
873 | '''Formats a Constraint or Function for use in a convenient way''' |
---|
874 | linelen = 45 |
---|
875 | s = [""] |
---|
876 | for var,val in RelDict.items(): |
---|
877 | if var.startswith('_'): continue |
---|
878 | if len(s[-1]) > linelen: s.append(' ') |
---|
879 | m = val |
---|
880 | if s[-1] != "" and m >= 0: |
---|
881 | s[-1] += ' + ' |
---|
882 | elif s[-1] != "": |
---|
883 | s[-1] += ' - ' |
---|
884 | m = abs(m) |
---|
885 | s[-1] += '%.3f*%s '%(m,var) |
---|
886 | if len(s[-1]) > linelen: s.append(' ') |
---|
887 | if RelVal is None: |
---|
888 | s[-1] += ' = New variable' |
---|
889 | else: |
---|
890 | s[-1] += ' = ' + RelVal |
---|
891 | s1 = '' |
---|
892 | for s2 in s: |
---|
893 | if s1 != '': s1 += '\n\t' |
---|
894 | s1 += s2 |
---|
895 | return s1 |
---|
896 | |
---|
897 | def VarRemapShow(varyList,inputOnly=False): |
---|
898 | '''List out the saved relationships. This should be done after the constraints have been |
---|
899 | defined using :func:`StoreEquivalence`, :func:`GroupConstraints` and :func:`GenerateConstraints`. |
---|
900 | |
---|
901 | :returns: a string containing the details of the contraint relationships |
---|
902 | ''' |
---|
903 | s = '' |
---|
904 | if len(fixedVarList) > 0: |
---|
905 | s += 'Fixed Variables:\n' |
---|
906 | for v in fixedVarList: |
---|
907 | s += ' ' + v + '\n' |
---|
908 | s += 'Variable mapping relations:\n' |
---|
909 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict |
---|
910 | for varlist,mapvars,multarr,invmultarr in zip( |
---|
911 | dependentParmList,indParmList,arrayList,invarrayList): |
---|
912 | for i,mv in enumerate(mapvars): |
---|
913 | if multarr is None: |
---|
914 | s += ' ' + str(mv) + ' is equivalent to parameter(s): ' |
---|
915 | j = 0 |
---|
916 | for v,m in zip(varlist,invmultarr): |
---|
917 | if debug: print 'v,m[0]: ',v,m[0] |
---|
918 | if j > 0: s += ' & ' |
---|
919 | j += 1 |
---|
920 | s += str(v) |
---|
921 | if m != 1: |
---|
922 | s += " / " + str(m[0]) |
---|
923 | s += '\n' |
---|
924 | continue |
---|
925 | s += ' %s = ' % mv |
---|
926 | j = 0 |
---|
927 | for m,v in zip(multarr[i,:],varlist): |
---|
928 | if m == 0: continue |
---|
929 | if j > 0: s += ' + ' |
---|
930 | j += 1 |
---|
931 | s += '(%s * %s)' % (m,v) |
---|
932 | if mv in varyList: s += ' VARY' |
---|
933 | s += '\n' |
---|
934 | if inputOnly: return s |
---|
935 | s += 'Inverse variable mapping relations:\n' |
---|
936 | for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList): |
---|
937 | for i,mv in enumerate(varlist): |
---|
938 | s += ' %s = ' % mv |
---|
939 | j = 0 |
---|
940 | for m,v in zip(invmultarr[i,:],mapvars): |
---|
941 | if m == 0: continue |
---|
942 | if j > 0: s += ' + ' |
---|
943 | j += 1 |
---|
944 | s += '(%s * %s)' % (m,v) |
---|
945 | s += '\n' |
---|
946 | return s |
---|
947 | |
---|
948 | def Dict2Deriv(varyList,derivDict,dMdv): |
---|
949 | '''Compute derivatives for Independent Parameters from the |
---|
950 | derivatives for the original parameters |
---|
951 | |
---|
952 | :param list varyList: a list of parameters names that will be varied |
---|
953 | |
---|
954 | :param dict derivDict: a dict containing derivatives for parameter values keyed by the |
---|
955 | parameter names. |
---|
956 | |
---|
957 | :param list dMdv: a Jacobian, as a list of np.array containing derivatives for dependent |
---|
958 | parameter computed from derivDict |
---|
959 | |
---|
960 | ''' |
---|
961 | global dependentParmList,arrayList,invarrayList,indParmList,invarrayList |
---|
962 | for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList): |
---|
963 | for i,name in enumerate(mapvars): |
---|
964 | # grouped variables: need to add in the derv. w/r |
---|
965 | # dependent variables to the independent ones |
---|
966 | if name not in varyList: continue # skip if independent var not varied |
---|
967 | if multarr is None: |
---|
968 | for v,m in zip(varlist,invmultarr): |
---|
969 | if debug: print 'start dMdv',dMdv[varyList.index(name)] |
---|
970 | if debug: print 'add derv',v,'/',m[0],'to derv',name,'add=',derivDict[v] / m[0] |
---|
971 | if m == 0: continue |
---|
972 | dMdv[varyList.index(name)] += derivDict[v] / m[0] |
---|
973 | else: |
---|
974 | for v,m in zip(varlist,multarr[i,:]): |
---|
975 | if debug: print 'start dMdv',dMdv[varyList.index(name)] |
---|
976 | if debug: print 'add derv',v,'*',m,'to derv',name,'add=',m * derivDict[v] |
---|
977 | if m == 0: continue |
---|
978 | dMdv[varyList.index(name)] += m * derivDict[v] |
---|
979 | |
---|
980 | def Map2Dict(parmDict,varyList): |
---|
981 | '''Create (or update) the Independent Parameters from the original |
---|
982 | set of Parameters |
---|
983 | |
---|
984 | Removes dependent variables from the varyList |
---|
985 | |
---|
986 | This should be done once, after the constraints have been |
---|
987 | defined using :func:`StoreEquivalence`, |
---|
988 | :func:`GroupConstraints` and :func:`GenerateConstraints` and |
---|
989 | before any variable refinement is done. This completes the parameter |
---|
990 | dictionary by defining independent parameters and it satisfies the |
---|
991 | constraint equations in the initial parameters |
---|
992 | |
---|
993 | :param dict parmDict: a dict containing parameter values keyed by the |
---|
994 | parameter names. |
---|
995 | This will contain updated values for both dependent and independent |
---|
996 | parameters after Dict2Map is called. It will also contain some |
---|
997 | unexpected entries of every constant value {'0':0.0} & {'1.0':1.0}, |
---|
998 | which do not cause any problems. |
---|
999 | |
---|
1000 | :param list varyList: a list of parameters names that will be varied |
---|
1001 | |
---|
1002 | |
---|
1003 | ''' |
---|
1004 | # process the Independent vars: remove dependent ones from varylist |
---|
1005 | # and then compute values for the Independent ones from their dependents |
---|
1006 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict |
---|
1007 | for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList): |
---|
1008 | for item in varlist: |
---|
1009 | try: |
---|
1010 | varyList.remove(item) |
---|
1011 | except ValueError: |
---|
1012 | pass |
---|
1013 | if multarr is None: continue |
---|
1014 | valuelist = [parmDict[var] for var in varlist] |
---|
1015 | parmDict.update(zip(mapvars, |
---|
1016 | np.dot(multarr,np.array(valuelist))) |
---|
1017 | ) |
---|
1018 | # now remove fixed variables from the varyList |
---|
1019 | global fixedVarList |
---|
1020 | for item in fixedVarList: |
---|
1021 | try: |
---|
1022 | varyList.remove(item) |
---|
1023 | except ValueError: |
---|
1024 | pass |
---|
1025 | # Set constrained parameters to their fixed values |
---|
1026 | parmDict.update(fixedDict) |
---|
1027 | |
---|
1028 | def Dict2Map(parmDict,varyList): |
---|
1029 | '''Applies the constraints defined using :func:`StoreEquivalence`, |
---|
1030 | :func:`GroupConstraints` and :func:`GenerateConstraints` by changing |
---|
1031 | values in a dict containing the parameters. This should be |
---|
1032 | done before the parameters are used for any computations |
---|
1033 | |
---|
1034 | :param dict parmDict: a dict containing parameter values keyed by the |
---|
1035 | parameter names. |
---|
1036 | This will contain updated values for both dependent and independent |
---|
1037 | parameters after Dict2Map is called. It will also contain some |
---|
1038 | unexpected entries of every constant value {'0':0.0} & {'1.0':1.0}, |
---|
1039 | which do not cause any problems. |
---|
1040 | |
---|
1041 | :param list varyList: a list of parameters names that will be varied |
---|
1042 | |
---|
1043 | ''' |
---|
1044 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict |
---|
1045 | # reset fixed values (should not be needed, but very quick) |
---|
1046 | # - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended |
---|
1047 | # not needed, but also not a problem - BHT |
---|
1048 | parmDict.update(fixedDict) |
---|
1049 | for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList): |
---|
1050 | #if invmultarr is None: continue |
---|
1051 | try: |
---|
1052 | valuelist = [parmDict[var] for var in mapvars] |
---|
1053 | except KeyError: |
---|
1054 | continue |
---|
1055 | parmDict.update(zip(varlist, |
---|
1056 | np.dot(invmultarr,np.array(valuelist))) |
---|
1057 | ) |
---|
1058 | |
---|
1059 | #====================================================================== |
---|
1060 | # internal routines follow (these routines are unlikely to be called |
---|
1061 | # from outside the module) |
---|
1062 | |
---|
1063 | def GramSchmidtOrtho(a,nkeep=0): |
---|
1064 | '''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to |
---|
1065 | find orthonormal unit vectors relative to first row. |
---|
1066 | |
---|
1067 | If nkeep is non-zero, the first nkeep rows in the array are not changed |
---|
1068 | |
---|
1069 | input: |
---|
1070 | arrayin: a 2-D non-singular square array |
---|
1071 | returns: |
---|
1072 | a orthonormal set of unit vectors as a square array |
---|
1073 | ''' |
---|
1074 | def proj(a,b): |
---|
1075 | 'Projection operator' |
---|
1076 | return a*(np.dot(a,b)/np.dot(a,a)) |
---|
1077 | for j in range(nkeep,len(a)): |
---|
1078 | for i in range(j): |
---|
1079 | a[j] -= proj(a[i],a[j]) |
---|
1080 | if np.allclose(np.linalg.norm(a[j]),0.0): |
---|
1081 | raise Exception,"Singular input to GramSchmidtOrtho" |
---|
1082 | a[j] /= np.linalg.norm(a[j]) |
---|
1083 | return a |
---|
1084 | |
---|
1085 | def _FillArray(sel,dict,collist,FillDiagonals=False): |
---|
1086 | '''Construct a n by n matrix (n = len(collist) |
---|
1087 | filling in the rows using the relationships defined |
---|
1088 | in the dictionaries selected by sel |
---|
1089 | |
---|
1090 | If FillDiagonals is True, diagonal elements in the |
---|
1091 | array are set to 1.0 |
---|
1092 | ''' |
---|
1093 | n = len(collist) |
---|
1094 | if FillDiagonals: |
---|
1095 | arr = np.eye(n) |
---|
1096 | else: |
---|
1097 | arr = np.zeros(2*[n,]) |
---|
1098 | # fill the top rows |
---|
1099 | for i,cnum in enumerate(sel): |
---|
1100 | for j,var in enumerate(collist): |
---|
1101 | arr[i,j] = dict[cnum].get(var,0) |
---|
1102 | return arr |
---|
1103 | |
---|
1104 | def _SwapColumns(i,m,v): |
---|
1105 | '''Swap columns in matrix m as well as the labels in v |
---|
1106 | so that element (i,i) is replaced by the first non-zero element in row i after that element |
---|
1107 | |
---|
1108 | Throws an exception if there are no non-zero elements in that row |
---|
1109 | ''' |
---|
1110 | for j in range(i+1,len(v)): |
---|
1111 | if not np.allclose(m[i,j],0): |
---|
1112 | m[:,(i,j)] = m[:,(j,i)] |
---|
1113 | v[i],v[j] = v[j],v[i] |
---|
1114 | return |
---|
1115 | else: |
---|
1116 | raise Exception,'Singular input' |
---|
1117 | |
---|
1118 | def _RowEchelon(m,arr,collist): |
---|
1119 | '''Convert the first m rows in Matrix arr to row-echelon form |
---|
1120 | exchanging columns in the matrix and collist as needed. |
---|
1121 | |
---|
1122 | throws an exception if the matrix is singular because |
---|
1123 | the first m rows are not linearly independent |
---|
1124 | ''' |
---|
1125 | n = len(collist) |
---|
1126 | for i in range(m): |
---|
1127 | if np.allclose(arr[i,i],0): |
---|
1128 | _SwapColumns(i,arr,collist) |
---|
1129 | arr[i,:] /= arr[i,i] # normalize row |
---|
1130 | # subtract current row from subsequent rows to set values to left of diagonal to 0 |
---|
1131 | for j in range(i+1,m): |
---|
1132 | arr[j,:] -= arr[i,:] * arr[j,i] |
---|
1133 | |
---|
1134 | if __name__ == "__main__": |
---|
1135 | parmdict = {} |
---|
1136 | constrDict = [ |
---|
1137 | {'0:12:Scale': 2.0, '0:11:Scale': 1.0, '0:14:Scale': 4.0, '0:13:Scale': 3.0, '0:0:Scale': 0.5}, |
---|
1138 | {'0:0:eA': 0.0}, |
---|
1139 | {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0}, |
---|
1140 | {'1::C(10,0,1)': 1.0, '2::C(10,0,1)': 1.0}, |
---|
1141 | {'1::AUiso:0': 1.0, '0::AUiso:0': 1.0}, |
---|
1142 | {'0::A0': 0.0} |
---|
1143 | ] |
---|
1144 | fixedList = ['5.0', '0', None, None, '1.0', '0'] |
---|
1145 | StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) |
---|
1146 | #StoreEquivalence('1::atomx:3',('2::atomx:3', ('2::atomz:3',2,), )) # error: dependent & independent vars mixed |
---|
1147 | #StoreEquivalence('1::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) # error: dependent vars repeated |
---|
1148 | #StoreEquivalence('0:1:eA',('0:0:eA',)) # error: equiv & fixed |
---|
1149 | #StoreEquivalence('0:99:Scale',('0:12:Scale',)) # error: equiv & constrained |
---|
1150 | #StoreEquivalence('0:12:Scale',('0:99:Scale',)) # error: equiv & constrained |
---|
1151 | varylist = ['2::atomx:3', |
---|
1152 | '2::C(10,6,1)', '1::C(10,6,1)', |
---|
1153 | '2::atomy:3', '2::atomz:3', |
---|
1154 | '0:12:Scale', '0:11:Scale', '0:14:Scale', '0:13:Scale', '0:0:Scale'] |
---|
1155 | # e,w = CheckConstraints([, |
---|
1156 | # [{'2:0:Scale': 1.0, '5:0:Scale': 1.0, '10:0:Scale': 1.0, '6:0:Scale': 1.0, '9:0:Scale': 1.0, '8:0:Scale': 1.0,# '3:0:Scale': 1.0, '4:0:Scale': 1.0, '7:0:Scale': 1.0, '1:0:Scale': 1.0, '0:0:Scale': 1.0}], |
---|
1157 | # ['1.0']) |
---|
1158 | # if e: print 'error=',e |
---|
1159 | # if w: print 'error=',w |
---|
1160 | # varyList = ['0::A0', '0::AUiso:0', '0::Afrac:1', '0::Afrac:2', '0::Afrac:3', '0::Afrac:4', '0::dAx:5', '0::dAy:5', '0::dAz:5', '0::AUiso:5', ':0:Back;0', ':0:Back;1', ':0:Back;2', ':0:Back;3', ':0:Back;4', ':0:Back;5', ':0:Back;6', ':0:Back;7', ':0:Back;8', ':0:Back;9', ':0:Back;10', ':0:Back;11', ':0:U', ':0:V', ':0:W', ':0:X', ':0:Y', ':0:Scale', ':0:DisplaceX', ':0:DisplaceY'] |
---|
1161 | # constrDict = [ |
---|
1162 | # {'0::Afrac:4': 24.0, '0::Afrac:1': 16.0, '0::Afrac:3': 24.0, '0::Afrac:2': 16.0}, |
---|
1163 | # {'0::Afrac:1': 1.0, '0::Afrac:2': 1.0}, |
---|
1164 | # {'0::Afrac:4': 1.0, '0::Afrac:3': 1.0}] |
---|
1165 | # fixedList = ['40.0', '1.0', '1.0'] |
---|
1166 | |
---|
1167 | errmsg, warnmsg = CheckConstraints(varylist,constrDict,fixedList) |
---|
1168 | if errmsg: |
---|
1169 | print "*** Error ********************" |
---|
1170 | print errmsg |
---|
1171 | if warnmsg: |
---|
1172 | print "*** Warning ********************" |
---|
1173 | print warnmsg |
---|
1174 | if errmsg or warnmsg: |
---|
1175 | sys.exit() |
---|
1176 | groups,parmlist = GroupConstraints(constrDict) |
---|
1177 | GenerateConstraints(groups,parmlist,varylist,constrDict,fixedList) |
---|
1178 | print VarRemapShow(varylist) |
---|
1179 | parmdict.update( { |
---|
1180 | '0:12:Scale': 1.0, '0:11:Scale': 1.0, '0:14:Scale': 1.0, '0:13:Scale': 1.0, '0:0:Scale': 2.0, |
---|
1181 | '0:0:eA': 0.0, |
---|
1182 | '2::C(10,6,1)': 0.2, '1::C(10,6,1)': 0.3, |
---|
1183 | '1::C(10,0,1)': 0.2, '2::C(10,0,1)': 0.3, |
---|
1184 | '1::AUiso:0': 0.02, '0::AUiso:0': 0.03, |
---|
1185 | '0::A0': 0.0, |
---|
1186 | '2::atomx:3':0.23,'2::atomy:3':-.23, '2::atomz:3':-0.11, |
---|
1187 | }) |
---|
1188 | print 'parmdict start',parmdict |
---|
1189 | print 'varylist start',varylist |
---|
1190 | before = parmdict.copy() |
---|
1191 | Map2Dict(parmdict,varylist) |
---|
1192 | print 'parmdict before and after Map2Dict' |
---|
1193 | print ' key / before / after' |
---|
1194 | for key in sorted(parmdict.keys()): |
---|
1195 | print ' '+key,'\t',before.get(key),'\t',parmdict[key] |
---|
1196 | print 'varylist after',varylist |
---|
1197 | before = parmdict.copy() |
---|
1198 | Dict2Map(parmdict,varylist) |
---|
1199 | print 'after Dict2Map' |
---|
1200 | print ' key / before / after' |
---|
1201 | for key in sorted(parmdict.keys()): |
---|
1202 | print ' '+key,'\t',before.get(key),'\t',parmdict[key] |
---|
1203 | # dMdv = len(varylist)*[0] |
---|
1204 | # deriv = {} |
---|
1205 | # for i,v in enumerate(parmdict.keys()): deriv[v]=i |
---|
1206 | # Dict2Deriv(varylist,deriv,dMdv) |
---|