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