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 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 | Relations are typically supplied as input to InputParse, where either |
---|
18 | of two keywords can be attached to a relationship: |
---|
19 | * VARY which means the generated parameter name will be included in |
---|
20 | the list of variables to be refined (varyList) |
---|
21 | * VARYFREE which means that all generated parameter names for a |
---|
22 | group, including the ones for generated relationships (see below) |
---|
23 | will be included in the list of variables to be refined (varyList) |
---|
24 | The case of VARY and VARYFREE is immaterial. |
---|
25 | |
---|
26 | Relations can also be supplied in the form of an equation: |
---|
27 | nx1 * Px + ny1 * Py +... = C1 |
---|
28 | where Cn is a constant. These equations define an algebraic |
---|
29 | constrant. The keyword VARY makes no sense when used with a constraint |
---|
30 | equation, but VARYFREE can be used. #RVD??? is VARYFREE required or default?? |
---|
31 | |
---|
32 | Unconstrained relations describe a new, independent, parameter, which |
---|
33 | is defined in terms of dependent parameters that are defined in the |
---|
34 | Model, while constrained relations effectively reduce the complexity |
---|
35 | of the Model by removing a degree of freedom. |
---|
36 | |
---|
37 | Relationships are grouped so that a set of dependent parameters appear |
---|
38 | in only one group (done in routine GroupConstraints.) Note that if a |
---|
39 | group contains relationships/equations that involve N dependent |
---|
40 | parameters, there must exist N-C new parameters, where C is the number |
---|
41 | of contraint equations in the group. Routine GenerateConstraints takes |
---|
42 | the output from InputParse and GroupConstraints generates the |
---|
43 | "missing" relationships and saves that information in the module's |
---|
44 | global variables. Each generated parameter is named |
---|
45 | sequentially using paramPrefix. |
---|
46 | |
---|
47 | Parameter redefinition is done by equating one (independent) parameter to several |
---|
48 | (now dependent) parameters, in algebraic form: |
---|
49 | P1 = n2 * P2 = n3 * P3 ... |
---|
50 | Each equality in the relationship reduces the complexity of the model |
---|
51 | by one potentially variable parameter. Input is provided to routine |
---|
52 | StoreEquivalence in the form of an independent parameter and a list of |
---|
53 | dependent parameters, optionally with a multiplier. |
---|
54 | |
---|
55 | Note that none of the dependent parameters in any constraint or |
---|
56 | reformulation can be refined (see dependentParmList, below). |
---|
57 | |
---|
58 | External Routines: |
---|
59 | To define a set of constrained and unconstrained relations, one |
---|
60 | calls InputParse, GroupConstraints and GenerateConstraints. It |
---|
61 | is best to supply all relations/equations in one call to these |
---|
62 | routines so that grouping is done correctly. |
---|
63 | |
---|
64 | To implement parameter redefinition, one calls |
---|
65 | StoreEquivalence. This should be called for every set of |
---|
66 | equivalence relationships. There is no harm in using |
---|
67 | StoreEquivalence with the same independent variable: |
---|
68 | StoreEquivalence('x',('y',)) |
---|
69 | StoreEquivalence('x',('z',)) |
---|
70 | (though StoreEquivalence('x',('y','z')) will run more |
---|
71 | efficiently) but mixing independent and dependent variables is |
---|
72 | problematic. This may not work properly: |
---|
73 | StoreEquivalence('x',('y',)) |
---|
74 | StoreEquivalence('y',('z',)) |
---|
75 | |
---|
76 | To show a summary of the parameter remapping, one calls |
---|
77 | VarRemapShow |
---|
78 | |
---|
79 | To determine values for the parameters created in this module, one |
---|
80 | calls Map2Dict. This will not apply contraints. |
---|
81 | |
---|
82 | To take values from the new independent parameters and constraints, |
---|
83 | one calls Dict2Map. This will apply contraints. |
---|
84 | |
---|
85 | Global Variables: |
---|
86 | dependentParmList: contains a list by group of lists of |
---|
87 | parameters used in the group. Note that parameters listed in |
---|
88 | dependentParmList should not be refined as they will not affect |
---|
89 | the model |
---|
90 | indParmList: a list of groups of Independent parameters defined in |
---|
91 | each group. This contains both parameters used in parameter |
---|
92 | redefinitions as well as names of generated new parameters. |
---|
93 | |
---|
94 | fixedVarList: a list of variables that have been 'fixed' |
---|
95 | by defining them as equal to a constant (::var: = 0). Note that |
---|
96 | the constant value is ignored at present. These variables are |
---|
97 | later removed from varyList which prevents them from being refined. |
---|
98 | Unlikely to be used externally. |
---|
99 | arrayList: a list by group of relationship matrices to relate |
---|
100 | parameters in dependentParmList to those in indParmList. Unlikely |
---|
101 | to be used externally. |
---|
102 | invarrayList: a list by group of relationship matrices to relate |
---|
103 | parameters in indParmList to those in dependentParmList. Unlikely |
---|
104 | to be used externally. |
---|
105 | fixedDict: a dictionary containing the fixed values corresponding |
---|
106 | to parameter equations. The dict key is an ascii string, but the |
---|
107 | dict value is a float. Unlikely to be used externally. |
---|
108 | """ |
---|
109 | |
---|
110 | import numpy as np |
---|
111 | import re |
---|
112 | # data used for constraints; |
---|
113 | debug = False # turns on printing as constraint input is processed |
---|
114 | # note that constraints are stored listed by contraint groups, where each constraint |
---|
115 | # group contains those parameters that must be handled together |
---|
116 | dependentParmList = [] # contains a list of parameters in each group |
---|
117 | # note that parameters listed in dependentParmList should not be refined |
---|
118 | arrayList = [] # a list of of relationship matrices |
---|
119 | invarrayList = [] # a list of inverse relationship matrices |
---|
120 | indParmList = [] # a list of names for the new parameters |
---|
121 | fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations |
---|
122 | # key is original ascii string, value is float |
---|
123 | fixedVarList = [] # List of variables that should not be refined |
---|
124 | |
---|
125 | # compile regular expressions used for parsing input |
---|
126 | rex_mult = re.compile('[+-]?[0-9.]+[eE]?[+-]?[0-9]*') |
---|
127 | rex_star = re.compile('\s*\*\s*') |
---|
128 | rex_var = re.compile('[a-zA-Z0-9_:]+') |
---|
129 | rex_plusminus = re.compile('\s*[+-=]\s*') |
---|
130 | rex_vary = re.compile('\s*Vary\s*', re.IGNORECASE) |
---|
131 | rex_varyfree = re.compile('(.*)\s*VaryFree\s*', re.IGNORECASE) |
---|
132 | |
---|
133 | # prefix for parameter names |
---|
134 | paramPrefix = "::constr:" |
---|
135 | consNum = 0 # number of the next constraint to be created |
---|
136 | |
---|
137 | def InitVars(): |
---|
138 | '''Initializes all constraint information''' |
---|
139 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,consNum |
---|
140 | dependentParmList = [] # contains a list of parameters in each group |
---|
141 | arrayList = [] # a list of of relationship matrices |
---|
142 | invarrayList = [] # a list of inverse relationship matrices |
---|
143 | indParmList = [] # a list of names for the new parameters |
---|
144 | fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations |
---|
145 | consNum = 0 # number of the next constraint to be created |
---|
146 | fixedVarList = [] |
---|
147 | |
---|
148 | def InputParse(mapList): |
---|
149 | '''Converts a set relationships used to remap parameters into composite |
---|
150 | parameters or to one or more constants. |
---|
151 | |
---|
152 | Input: |
---|
153 | mapList: a list of character strings where each string defines |
---|
154 | a relationship in form: |
---|
155 | ('<const11>*<prm1> [+-] <const12>*<prm2> [+-] ... = <value> [VaryFree]', |
---|
156 | '<const21>*<prm1> [+-] <const22>*<prm2> [+-] ... [Vary/VaryFree]', |
---|
157 | ...) |
---|
158 | these define either a constraint or a new independent parameter, |
---|
159 | where <constXX> is a constant containing an optional sign, digits, |
---|
160 | optionally a decimal point and/or an exponent, prefixed by e or E, |
---|
161 | and <prmN> is the name of a parameter defined in the Model. |
---|
162 | Parameters can be included in any order. Multiplying a parameter |
---|
163 | by 0 causes that parameter to be included in a group of |
---|
164 | relationships. Note that if the same parameter is included twice |
---|
165 | in a relationship/equation, the coefficients are added. |
---|
166 | |
---|
167 | When an equality is defined, a constant, <value>, is |
---|
168 | included. This then describes a constraint equation. |
---|
169 | |
---|
170 | Vary is an optional keyword that indicates the indicated remapped |
---|
171 | parameter should be varied (note, case insensitive). Should be the |
---|
172 | last item on the line. Makes no sense with an equation. |
---|
173 | |
---|
174 | VaryFree is an optional keyword that indicates all possible |
---|
175 | remapped parameters should be varied (note, case |
---|
176 | insensitive). Makes most sense with a constraint equation. Should |
---|
177 | be the last item on the line. |
---|
178 | |
---|
179 | returns |
---|
180 | constrDict: a list with a dict for each item in mapList that |
---|
181 | defines the relationship. The keys are parameter names and the |
---|
182 | values are the multiplier for the parameter name |
---|
183 | constrFlag: a list for each item in mapList with a list contains |
---|
184 | 'Vary' and/or 'VaryFree' |
---|
185 | fixedList: a list for each item in mapList. Contains the value |
---|
186 | (as a string) for each contraint equation or None for a |
---|
187 | constraint relationship. |
---|
188 | ''' |
---|
189 | i = 0 |
---|
190 | constrDict = [] |
---|
191 | fixedList = [] |
---|
192 | constrFlag = [] |
---|
193 | if debug: print 50*'-','\n(Input)' |
---|
194 | for line in mapList: |
---|
195 | inputline = line[:] |
---|
196 | line = line.strip() |
---|
197 | if len(line) == 0: continue # empty lines are ignored |
---|
198 | constrDict.append({}) |
---|
199 | constrFlag.append([]) |
---|
200 | fixedList.append(None) |
---|
201 | i += 1 |
---|
202 | fixedval = None |
---|
203 | if debug: print '\t',line |
---|
204 | try: |
---|
205 | j = 0 |
---|
206 | sign = 1.0 |
---|
207 | while len(line): |
---|
208 | j += 1 |
---|
209 | m = sign * float(rex_mult.match(line).group()) |
---|
210 | line = line[rex_mult.match(line).end():] |
---|
211 | j += 1 |
---|
212 | line = line[rex_star.match(line).end():] |
---|
213 | j += 1 |
---|
214 | v = rex_var.match(line).group() |
---|
215 | line = line[rex_var.match(line).end():] |
---|
216 | #print m,'times',v |
---|
217 | #if v not in varlist: varlist.append(v) |
---|
218 | if constrDict[i-1].get(v) == None: |
---|
219 | constrDict[i-1][v] = m |
---|
220 | else: |
---|
221 | constrDict[i-1][v] += m |
---|
222 | if len(line.strip()) == 0: break |
---|
223 | j += 1 |
---|
224 | # advance to next separator (+-=) |
---|
225 | pm = rex_plusminus.match(line) |
---|
226 | if pm != None: |
---|
227 | line = line[pm.end():] |
---|
228 | pm = pm.group() |
---|
229 | else: |
---|
230 | pm = '' |
---|
231 | if pm.strip() == '+': |
---|
232 | sign = 1.0 |
---|
233 | elif pm.strip() == '-': |
---|
234 | sign = -1.0 |
---|
235 | elif pm.strip() == '=': |
---|
236 | # found a fixed value, also check for a VaryFree flag |
---|
237 | if rex_varyfree.match(line): |
---|
238 | constrFlag[-1].append('VaryFree') |
---|
239 | #fixedval = float(rex_varyfree.split(line)[1]) |
---|
240 | fixedval = rex_varyfree.split(line)[1].strip() |
---|
241 | else: |
---|
242 | #fixedval = float(line.strip()) |
---|
243 | fixedval = line.strip() |
---|
244 | fixedList[-1] = fixedval |
---|
245 | line = "" |
---|
246 | elif rex_varyfree.match(line): |
---|
247 | constrFlag[-1].append('VaryFree') |
---|
248 | line = line[rex_varyfree.match(line).end():] |
---|
249 | elif rex_vary.match(line): |
---|
250 | constrFlag[-1].append('Vary') |
---|
251 | line = line[rex_vary.match(line).end():] |
---|
252 | else: |
---|
253 | # something else is on the line, but not a keyword |
---|
254 | raise Exception,'Error in line '+str(inputline)+'\nat '+str(line) |
---|
255 | except SyntaxError: |
---|
256 | if debug: print 'Error in line',i,'token',j,'@','"'+line+'"' |
---|
257 | raise Exception,'Error in line %d token %d, beginning with %s'% (i,j, line) |
---|
258 | |
---|
259 | if debug: # on debug, show what is parsed in a semi-readable |
---|
260 | print 50*'-','\n(parsed relationship/equation & flag)' |
---|
261 | for i in range(len(constrDict)): |
---|
262 | flags = '' |
---|
263 | for f in constrFlag[i]: |
---|
264 | if flags != '': |
---|
265 | flags = flags + ' + ' + f |
---|
266 | else: |
---|
267 | flags = f |
---|
268 | if fixedList[i] is None: |
---|
269 | print '#',i+1,constrDict[i],'\t',constrFlag[i] |
---|
270 | else: |
---|
271 | print '#',i+1,constrDict[i],'=',fixedList[i],'\t',constrFlag[i] |
---|
272 | |
---|
273 | return constrDict,constrFlag,fixedList |
---|
274 | |
---|
275 | def GroupConstraints(constrDict): |
---|
276 | """divide the constraints into groups that share no parameters. |
---|
277 | Input |
---|
278 | constrDict: a list of dicts defining relationships/constraints |
---|
279 | (see InputParse) |
---|
280 | Returns two lists of lists: |
---|
281 | a list of relationship list indices for each group and |
---|
282 | a list containing the parameters used in each group |
---|
283 | """ |
---|
284 | assignedlist = [] # relationships that have been used |
---|
285 | groups = [] # contains a list of grouplists |
---|
286 | ParmList = [] |
---|
287 | for i,consi in enumerate(constrDict): |
---|
288 | if i in assignedlist: continue # already in a group, skip |
---|
289 | # starting a new group |
---|
290 | grouplist = [i,] |
---|
291 | assignedlist.append(i) |
---|
292 | groupset = set(consi.keys()) |
---|
293 | changes = True # always loop at least once |
---|
294 | while(changes): # loop until we can't find anything to add to the current group |
---|
295 | changes = False # but don't loop again unless we find something |
---|
296 | for j,consj in enumerate(constrDict): |
---|
297 | if j in assignedlist: continue # already in a group, skip |
---|
298 | if len(set(consj.keys()) & groupset) > 0: # true if this needs to be added |
---|
299 | changes = True |
---|
300 | grouplist.append(j) |
---|
301 | assignedlist.append(j) |
---|
302 | groupset = groupset | set(consj.keys()) |
---|
303 | group = sorted(grouplist) |
---|
304 | varlist = sorted(list(groupset)) |
---|
305 | groups.append(group) |
---|
306 | ParmList.append(varlist) |
---|
307 | return groups,ParmList |
---|
308 | |
---|
309 | def GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList): |
---|
310 | '''Takes a list of relationship entries comprising a group of constraints and |
---|
311 | builds the relationship lists and their inverse and stores them in global variables |
---|
312 | ''' |
---|
313 | global dependentParmList,arrayList,invarrayList,indParmList,consNum |
---|
314 | for group,varlist in zip(groups,parmlist): |
---|
315 | VaryFree = False |
---|
316 | for row in group: |
---|
317 | if 'VaryFree' in constrFlag[row]: VaryFree = True |
---|
318 | if len(varlist) == 1: |
---|
319 | #print "Debug: Fixing parameter (%s)" % (varlist[0]) |
---|
320 | fixedVarList.append(varlist[0]) |
---|
321 | continue |
---|
322 | arr = MakeArray(constrDict, group, varlist) |
---|
323 | constrArr = FillInMissingRelations(arr,len(group)) |
---|
324 | mapvar = [] |
---|
325 | group = group[:] |
---|
326 | for i in range(len(varlist)): |
---|
327 | vary = VaryFree |
---|
328 | if len(group) > 0: |
---|
329 | rel = group.pop(0) |
---|
330 | fixedval = fixedList[rel] |
---|
331 | if 'Vary' in constrFlag[rel]: vary = True |
---|
332 | else: |
---|
333 | fixedval = None |
---|
334 | if fixedval is None: |
---|
335 | varname = paramPrefix + str(consNum) |
---|
336 | mapvar.append(varname) |
---|
337 | consNum += 1 |
---|
338 | if vary: varyList.append(varname) |
---|
339 | else: |
---|
340 | mapvar.append(fixedval) |
---|
341 | dependentParmList.append(varlist) |
---|
342 | arrayList.append(constrArr) |
---|
343 | invarrayList.append(np.linalg.inv(constrArr)) |
---|
344 | indParmList.append(mapvar) |
---|
345 | # setup dictionary containing the fixed values |
---|
346 | global fixedDict |
---|
347 | # key is original ascii string, value is float |
---|
348 | for fixedval in fixedList: |
---|
349 | if fixedval: |
---|
350 | fixedDict[fixedval] = float(fixedval) |
---|
351 | |
---|
352 | if debug: # on debug, show what is parsed & generated, semi-readable |
---|
353 | print 50*'-' |
---|
354 | for group,varlist,multarr,inv,mapvar in zip(groups,parmlist,arrayList,invarrayList,indParmList): |
---|
355 | print '\n*** relation(s) in group:',group,'\tvars in group:',varlist |
---|
356 | print 'new parameters:', mapvar |
---|
357 | print 'Input relationship matrix' |
---|
358 | print multarr[:len(group)] |
---|
359 | added = len(group) - len(varlist) |
---|
360 | if added < 0: |
---|
361 | print 'added relationship rows' |
---|
362 | print multarr[added:] |
---|
363 | print 'Inverse relationship matrix' |
---|
364 | print inv |
---|
365 | |
---|
366 | def StoreEquivalence(independentVar,dependentList): |
---|
367 | '''Takes a list of dependent parameter(s) and stores their |
---|
368 | relationship to a single independent parameter (independentVar) |
---|
369 | input: |
---|
370 | independentVar -- name of parameter that will set others (may |
---|
371 | be varied) |
---|
372 | dependentList -- list of parameter that will set from |
---|
373 | independentVar (may not be varied). Each item can be a parameter |
---|
374 | name or a tuple containing a name and multiplier: |
---|
375 | ('parm1',('parm2',.5),) |
---|
376 | ''' |
---|
377 | |
---|
378 | global dependentParmList,arrayList,invarrayList,indParmList |
---|
379 | mapList = [] |
---|
380 | multlist = [] |
---|
381 | for var in dependentList: |
---|
382 | if isinstance(var, basestring): |
---|
383 | mult = 1.0 |
---|
384 | elif len(var) == 2: |
---|
385 | var,mult = var |
---|
386 | else: |
---|
387 | raise Exception, "Cannot parse "+repr(var) + " as var or (var,multiplier)" |
---|
388 | mapList.append(var) |
---|
389 | multlist.append(tuple((mult,))) |
---|
390 | # added relationships to stored values |
---|
391 | arrayList.append(None) |
---|
392 | invarrayList.append(np.array(multlist)) |
---|
393 | indParmList.append(tuple((independentVar,))) |
---|
394 | dependentParmList.append(mapList) |
---|
395 | return |
---|
396 | |
---|
397 | def SetVaryFlags(varyList): |
---|
398 | '''Adds independent variables to the varyList provided that all |
---|
399 | dependent variables are being varied. |
---|
400 | Ignores independent variables where no dependent variables are |
---|
401 | being varied. |
---|
402 | Returns a non-empty text message where some but not all dependent |
---|
403 | variables are being varied. |
---|
404 | ''' |
---|
405 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict |
---|
406 | msg = "" |
---|
407 | for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList): |
---|
408 | for mv in mapvars: |
---|
409 | varied = [] |
---|
410 | notvaried = [] |
---|
411 | if mv in fixedDict: continue |
---|
412 | if multarr is None: |
---|
413 | if mv in varyList: |
---|
414 | varied.append(mv) |
---|
415 | else: |
---|
416 | notvaried.append(mv) |
---|
417 | for v in varlist: |
---|
418 | if v in varyList: |
---|
419 | varied.append(v) |
---|
420 | else: |
---|
421 | notvaried.append(v) |
---|
422 | if len(varied) > 0 and len(notvaried) > 0: |
---|
423 | if msg != "": msg += "\n" |
---|
424 | msg += "Error: inconsistent use of parameter " + mv |
---|
425 | msg += "\n varied: " |
---|
426 | for v in varied: msg += v |
---|
427 | msg += "\n not varied: " |
---|
428 | for v in notvaried: msg += v |
---|
429 | elif len(varied) > 0 and multarr is not None: |
---|
430 | varyList.append(mv) |
---|
431 | return msg |
---|
432 | |
---|
433 | def GetDependentVars(): |
---|
434 | '''Return a list of dependent variables: e.g. variables that are |
---|
435 | constrained in terms of other variables''' |
---|
436 | dependentVars = [] |
---|
437 | global dependentParmList |
---|
438 | for lst in dependentParmList: |
---|
439 | for itm in lst: dependentVars.append(itm) |
---|
440 | return dependentVars |
---|
441 | |
---|
442 | def GetIndependentVars(): |
---|
443 | '''Return a list of independent variables: e.g. variables that are |
---|
444 | created by constrains of other variables''' |
---|
445 | independentVars = [] |
---|
446 | global indParmList,fixedDict |
---|
447 | for lst in indParmList: |
---|
448 | for name in lst: |
---|
449 | if name in fixedDict: continue |
---|
450 | independentVars.append(name) |
---|
451 | return independentVars |
---|
452 | |
---|
453 | def PrintIndependentVars(parmDict,varyList,sigDict,PrintAll=False): |
---|
454 | '''Print the values and uncertainties on the independent variables''' |
---|
455 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict |
---|
456 | printlist = [] |
---|
457 | mapvars = GetIndependentVars() |
---|
458 | for i,name in enumerate(mapvars): |
---|
459 | if name in fixedDict: continue |
---|
460 | if PrintAll or name in varyList: |
---|
461 | sig = sigDict.get(name) |
---|
462 | printlist.append([name,parmDict[name],sig]) |
---|
463 | if len(printlist) == 0: return |
---|
464 | s1 = '' |
---|
465 | print 130*'-' |
---|
466 | print "Variables generated by constraints" |
---|
467 | printlist.append(3*[None]) |
---|
468 | for name,val,esd in printlist: |
---|
469 | if len(s1) > 110 or name is None: |
---|
470 | print |
---|
471 | print s1 |
---|
472 | print s2 |
---|
473 | print s3 |
---|
474 | s1 = '' |
---|
475 | if name is None: break |
---|
476 | if s1 == "": |
---|
477 | s1 = ' name :' |
---|
478 | s2 = ' value :' |
---|
479 | s3 = ' sig :' |
---|
480 | s1 += '%12s' % (name) |
---|
481 | s2 += '%12.6f' % (val) |
---|
482 | if esd is None: |
---|
483 | s3 += '%12s' % ('n/a') |
---|
484 | else: |
---|
485 | s3 += '%12.6f' % (esd) |
---|
486 | |
---|
487 | def ComputeDepESD(covMatrix,varyList,parmDict): |
---|
488 | '''Compute uncertainties for dependent parameters from independent ones |
---|
489 | returns a dictionary containing the esd values for dependent parameters |
---|
490 | ''' |
---|
491 | sigmaDict = {} |
---|
492 | for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList): |
---|
493 | if invmultarr is None: continue # probably not needed |
---|
494 | valuelist = [parmDict[var] for var in mapvars] |
---|
495 | # get the v-covar matrix for independent parameters |
---|
496 | vcov = np.zeros((len(mapvars),len(mapvars))) |
---|
497 | for i1,name1 in enumerate(mapvars): |
---|
498 | if name1 not in varyList: continue |
---|
499 | iv1 = varyList.index(name1) |
---|
500 | for i2,name2 in enumerate(mapvars): |
---|
501 | if name2 not in varyList: continue |
---|
502 | iv2 = varyList.index(name2) |
---|
503 | vcov[i1][i2] = covMatrix[iv1][iv2] |
---|
504 | # vec is the vector that multiplies each of the independent values |
---|
505 | for v,vec in zip(varlist,invmultarr): |
---|
506 | sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec))) |
---|
507 | return sigmaDict |
---|
508 | |
---|
509 | def VarRemapShow(varyList): |
---|
510 | '''List out the saved relationships. |
---|
511 | Returns a string containing the details. |
---|
512 | ''' |
---|
513 | s = '' |
---|
514 | if len(fixedVarList) > 0: |
---|
515 | s += 'Fixed Variables:\n' |
---|
516 | for v in fixedVarList: |
---|
517 | s += ' ' + v + '\n' |
---|
518 | s += 'Variable mapping relations:\n' |
---|
519 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict |
---|
520 | for varlist,mapvars,multarr,invmultarr in zip( |
---|
521 | dependentParmList,indParmList,arrayList,invarrayList): |
---|
522 | for i,mv in enumerate(mapvars): |
---|
523 | if multarr is None: |
---|
524 | s += ' ' + str(mv) + ' defines parameter(s): ' |
---|
525 | j = 0 |
---|
526 | for v,m in zip(varlist,invmultarr): |
---|
527 | print v,m[0] |
---|
528 | if j > 0: s += ' & ' |
---|
529 | j += 1 |
---|
530 | s += str(v) |
---|
531 | if m != 1: |
---|
532 | s += " / " + str(m[0]) |
---|
533 | s += '\n' |
---|
534 | continue |
---|
535 | s += ' %s = ' % mv |
---|
536 | j = 0 |
---|
537 | for m,v in zip(multarr[i,:],varlist): |
---|
538 | if m == 0: continue |
---|
539 | if j > 0: s += ' + ' |
---|
540 | j += 1 |
---|
541 | s += '(%s * %s)' % (m,v) |
---|
542 | if mv in varyList: s += ' VARY' |
---|
543 | s += '\n' |
---|
544 | s += 'Inverse variable mapping relations:\n' |
---|
545 | for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList): |
---|
546 | for i,mv in enumerate(varlist): |
---|
547 | s += ' %s = ' % mv |
---|
548 | j = 0 |
---|
549 | for m,v in zip(invmultarr[i,:],mapvars): |
---|
550 | if m == 0: continue |
---|
551 | if j > 0: s += ' + ' |
---|
552 | j += 1 |
---|
553 | s += '(%s * %s)' % (m,v) |
---|
554 | s += '\n' |
---|
555 | return s |
---|
556 | |
---|
557 | def Dict2Deriv(varyList,derivDict,dMdv): |
---|
558 | '''Compute derivatives for Independent Parameters from the |
---|
559 | derivatives for the original parameters |
---|
560 | ''' |
---|
561 | global dependentParmList,arrayList,invarrayList,indParmList,invarrayList |
---|
562 | for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList): |
---|
563 | for i,name in enumerate(mapvars): |
---|
564 | # grouped variables: need to add in the derv. w/r |
---|
565 | # dependent variables to the independent ones |
---|
566 | if name not in varyList: continue # if independent var not varied |
---|
567 | for m,v in zip(invmultarr[:,i],varlist): |
---|
568 | # print 'add derv',v,'*',m,'to derv',name |
---|
569 | if m == 0: continue |
---|
570 | dMdv[varyList.index(name)] += m * derivDict[v] |
---|
571 | |
---|
572 | def Map2Dict(parmDict,varyList): |
---|
573 | '''Create (or update) the Independent Parameters from the original |
---|
574 | set of Parameters |
---|
575 | |
---|
576 | Removes dependent variables from the varyList |
---|
577 | |
---|
578 | This should be done once, before any variable refinement is done |
---|
579 | to complete the parameter dictionary with the Independent Parameters |
---|
580 | ''' |
---|
581 | # process the Independent vars: remove dependent ones from varylist |
---|
582 | # and then compute values for the Independent ones from their dependents |
---|
583 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict |
---|
584 | for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList): |
---|
585 | for item in varlist: |
---|
586 | try: |
---|
587 | varyList.remove(item) |
---|
588 | except ValueError: |
---|
589 | pass |
---|
590 | if multarr is None: continue |
---|
591 | valuelist = [parmDict[var] for var in varlist] |
---|
592 | parmDict.update(zip(mapvars, |
---|
593 | np.dot(multarr,np.array(valuelist))) |
---|
594 | ) |
---|
595 | # now remove fixed variables from the varyList |
---|
596 | global fixedVarList |
---|
597 | for item in fixedVarList: |
---|
598 | try: |
---|
599 | varyList.remove(item) |
---|
600 | except ValueError: |
---|
601 | pass |
---|
602 | # Set constrained parameters to their fixed values |
---|
603 | parmDict.update(fixedDict) |
---|
604 | |
---|
605 | def Dict2Map(parmDict,varyList): |
---|
606 | #I think this needs fixing to update parmDict with new values |
---|
607 | # from the constraints based on what is in varyList - RVD |
---|
608 | '''Convert the remapped values back to the original parameters |
---|
609 | |
---|
610 | This should be done to apply constraints to parameter values (use |
---|
611 | Map2Dict first). It should be done as part of the Model function |
---|
612 | evaluation, before any computation is done |
---|
613 | ''' |
---|
614 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict |
---|
615 | # reset fixed values (should not be needed, but very quick) |
---|
616 | # - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended |
---|
617 | parmDict.update(fixedDict) |
---|
618 | for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList): |
---|
619 | if invmultarr is None: continue |
---|
620 | valuelist = [parmDict[var] for var in mapvars] |
---|
621 | parmDict.update(zip(varlist, |
---|
622 | np.dot(invmultarr,np.array(valuelist))) |
---|
623 | ) |
---|
624 | |
---|
625 | #====================================================================== |
---|
626 | # internal routines follow (these routines are unlikely to be called |
---|
627 | # from outside the module) |
---|
628 | def FillInMissingRelations(arrayin,nvars): |
---|
629 | '''Fill in unspecified rows in array with non-colinear unit vectors |
---|
630 | arrayin is a square array, where only the first nvars rows are defined. |
---|
631 | |
---|
632 | Returns a new array where all rows are defined |
---|
633 | |
---|
634 | Throws an exception if there is no non-signular result |
---|
635 | (likely because two or more input rows are co-linear) |
---|
636 | ''' |
---|
637 | a = arrayin.copy() |
---|
638 | n = nvars |
---|
639 | # fill in the matrix with basis vectors not colinear with any of the starting set |
---|
640 | for i in range(n,len(a)): |
---|
641 | try: |
---|
642 | a[n] = VectorProduct(a[:n]) |
---|
643 | except: |
---|
644 | raise Exception,"VectorProduct failed, singular input?" |
---|
645 | n += 1 |
---|
646 | # use the G-S algorithm to compute a complete set of unit vectors orthogonal |
---|
647 | # to the first in array |
---|
648 | gs = GramSchmidtOrtho(a) |
---|
649 | n = nvars |
---|
650 | # now replace the created vectors with the ones least correlated to the |
---|
651 | # initial set |
---|
652 | vlist = [v for v in gs[1:]] # drop the first row |
---|
653 | for j in range(n,len(a)): |
---|
654 | minavg = None |
---|
655 | droplist = [] |
---|
656 | for k in range(len(vlist)): |
---|
657 | v = vlist[k] |
---|
658 | avgcor = 0 |
---|
659 | for i in range(n): |
---|
660 | cor = np.dot(a[i],vlist[k])/(np.linalg.norm(a[i]) * np.linalg.norm(vlist[k]) ) |
---|
661 | if np.allclose(cor, 1.0): |
---|
662 | droplist.append(k) # found a vector co-linear to one we already have |
---|
663 | # don't need it again |
---|
664 | #vlist.pop(k) |
---|
665 | break |
---|
666 | avgcor += cor |
---|
667 | else: |
---|
668 | if minavg == None: |
---|
669 | minavg = abs(avgcor/n) |
---|
670 | minnum = k |
---|
671 | elif abs(avgcor/n) < minavg: |
---|
672 | minavg = abs(avgcor/n) |
---|
673 | minnum = k |
---|
674 | if minavg == None: |
---|
675 | raise Exception,"Failed to find a non-colinear G-S vector for row %d. Should not happen!" % n |
---|
676 | a[j] = vlist[minnum] |
---|
677 | droplist.append(minnum) # don't need to use this vector again |
---|
678 | for i in sorted(droplist,reverse=True): |
---|
679 | vlist.pop(i) |
---|
680 | n += 1 |
---|
681 | return a |
---|
682 | |
---|
683 | |
---|
684 | def MakeArray(constrDict, group, varlist): |
---|
685 | """Convert the multipliers in a constraint group to an array of |
---|
686 | coefficients and place in a square numpy array, adding empty rows as |
---|
687 | needed. |
---|
688 | |
---|
689 | Throws an exception if some sanity checks fail. |
---|
690 | """ |
---|
691 | # do some error checks |
---|
692 | if len(varlist) < len(group): # too many relationships -- no can do |
---|
693 | raise Exception,'The number of relationships (%d) exceeds the number of parameters (%d):\n\t%s\n\t%s'% ( |
---|
694 | len(group),len(varlist),group,varlist) |
---|
695 | # put the coefficients into an array |
---|
696 | multarr = np.zeros(2*[len(varlist),]) |
---|
697 | i = 0 |
---|
698 | for cnum in group: |
---|
699 | cdict = constrDict[cnum] |
---|
700 | j = 0 |
---|
701 | for var in varlist: |
---|
702 | m = cdict.get(var) |
---|
703 | if m != None: |
---|
704 | multarr[i,j] = m |
---|
705 | j += 1 |
---|
706 | i += 1 |
---|
707 | return multarr |
---|
708 | |
---|
709 | def GramSchmidtOrtho(arrayin): |
---|
710 | '''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to |
---|
711 | find orthonormal unit vectors relative to first row. |
---|
712 | input: |
---|
713 | arrayin: a 2-D non-signular square array |
---|
714 | returns: |
---|
715 | a orthonormal set of unit vectors as a square array |
---|
716 | ''' |
---|
717 | def proj(a,b): |
---|
718 | 'Projection operator' |
---|
719 | return a*(np.dot(a,b)/np.dot(a,a)) |
---|
720 | a = arrayin.copy() |
---|
721 | for j in range (len(a)): |
---|
722 | for i in range(j): |
---|
723 | a[j] = a[j] - proj(a[i],a[j]) |
---|
724 | if np.allclose(np.linalg.norm(a[j]),0.0): |
---|
725 | raise Exception,"Singular input to GramSchmidtOrtho" |
---|
726 | a[j] = a[j]/np.linalg.norm(a[j]) |
---|
727 | return a |
---|
728 | |
---|
729 | def VectorProduct(vl): |
---|
730 | '''Evaluate the n-dimensional "cross" product of the list of vectors in vl |
---|
731 | vl can be any length. The length of each vector is arbitrary, but |
---|
732 | all must be the same length. See http://en.wikipedia.org/wiki/Levi-Civita_symbol |
---|
733 | |
---|
734 | This appears to return a vector perpendicular to the supplied set. |
---|
735 | |
---|
736 | Throws an exception if a vector can not be obtained because the input set of |
---|
737 | vectors is already complete or contains any redundancies. |
---|
738 | |
---|
739 | Uses LeviCitvitaMult recursively to obtain all permutations of vector elements |
---|
740 | ''' |
---|
741 | i = 0 |
---|
742 | newvec = [] |
---|
743 | for e in vl[0]: |
---|
744 | i += 1 |
---|
745 | tl = [([i,],1),] |
---|
746 | seq = LeviCitvitaMult(tl ,vl) |
---|
747 | tsum = 0 |
---|
748 | for terms,prod in seq: |
---|
749 | tsum += EvalLCterm(terms) * prod |
---|
750 | newvec.append(tsum) |
---|
751 | if np.allclose(newvec,0.0): |
---|
752 | raise Exception,"VectorProduct failed, singular or already complete input" |
---|
753 | return newvec |
---|
754 | |
---|
755 | def LeviCitvitaMult(tin,vl): |
---|
756 | '''Recursion formula to compute all permutations of vector element products |
---|
757 | The first term in each list is the indices of the Levi-Civita symbol, note |
---|
758 | that if an index is repeated, the value is zero, so the element is not included. |
---|
759 | The second term is the product of the vector elements |
---|
760 | ''' |
---|
761 | v = vl[0] |
---|
762 | vl1 = vl[1:] |
---|
763 | |
---|
764 | j = 0 |
---|
765 | tl = [] |
---|
766 | for e in vl[0]: |
---|
767 | j += 1 |
---|
768 | for ind,vals in tin: |
---|
769 | if j in ind: continue |
---|
770 | tl.append((ind+[j,],vals*e)) |
---|
771 | if len(vl1): |
---|
772 | return LeviCitvitaMult(tl,vl1) |
---|
773 | else: |
---|
774 | return tl |
---|
775 | |
---|
776 | def EvalLCterm(term): |
---|
777 | '''Evaluate the Levi-Civita symbol Epsilon(i,j,k,...)''' |
---|
778 | p = 1 |
---|
779 | for i in range(len(term)): |
---|
780 | for j in range(i+1,len(term)): |
---|
781 | p *= (term[j] - term[i]) |
---|
782 | if p == 0: return 0 |
---|
783 | return p/abs(p) |
---|
784 | |
---|
785 | |
---|
786 | if __name__ == "__main__": |
---|
787 | import sys |
---|
788 | #remap = MapParameters() # create an object (perhaps better as a module) |
---|
789 | #remap.debug = 1 |
---|
790 | #debug = 1 |
---|
791 | parmdict = {} |
---|
792 | StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) |
---|
793 | varylist = ['2::atomx:3',] |
---|
794 | #print VarRemapShow(varylist) |
---|
795 | #msg = SetVaryFlags(varylist) |
---|
796 | #if msg != "": print msg |
---|
797 | varylist = ['2::atomx:3', '2::atomy:3', '2::atomz:3',] |
---|
798 | print VarRemapShow(varylist) |
---|
799 | msg = SetVaryFlags(varylist) |
---|
800 | if msg != "": print msg |
---|
801 | #parmdict = { |
---|
802 | # '2::atomx:3':.1 , |
---|
803 | # '2::atomy:3':.2 , # conflicts with constraint |
---|
804 | # '2::atomz:3':.3 , |
---|
805 | #} |
---|
806 | |
---|
807 | varylist = [] |
---|
808 | mapList1 = [ |
---|
809 | "1 * 2::atomz:3 = 0", |
---|
810 | "1*p1 + 2e0*p2 - 1.0*p3", |
---|
811 | "1*p4 + 1*p7", |
---|
812 | "1*p1+2*p2-3.0*p2 VARY", |
---|
813 | ] |
---|
814 | mapList = [ |
---|
815 | "1 * 2::atomz:3 = 0", |
---|
816 | "1*p1 + 2e0*p2 - 1.0*p3", |
---|
817 | "1*p4 + 1*p7", |
---|
818 | "1*p1+2*p2-3.0*p2 VARY", |
---|
819 | "1*p21 + 0*p22 - 0*p23 + 0*p24 varyfree", |
---|
820 | "1*p5 + 2*p6 = 1 varyfree", |
---|
821 | "-1 * p6 + 1*p5", |
---|
822 | "-10e-1 * p1 - -2*p2 + 3.0*p4", |
---|
823 | ] |
---|
824 | constrDict,constrFlag,fixedList = InputParse(mapList1) |
---|
825 | #print constrDict |
---|
826 | #print constrFlag |
---|
827 | #print fixedList |
---|
828 | groups,parmlist = GroupConstraints(constrDict) |
---|
829 | GenerateConstraints(groups,parmlist,varylist,constrDict,constrFlag,fixedList) |
---|
830 | print VarRemapShow(varylist) |
---|
831 | sys.exit() |
---|
832 | parmdict.update( {'p1':1,'p2':2,'p3':3,'p4':4, |
---|
833 | 'p6':6,'p5':5, # conflicts with constraint |
---|
834 | 'p7':7, |
---|
835 | "p21":.1,"p22":.2,"p23":.3,"p24":.4, |
---|
836 | }) |
---|
837 | #print 'parmdict start',parmdict |
---|
838 | before = parmdict.copy() |
---|
839 | Map2Dict(parmdict,[]) |
---|
840 | print 'parmdict before and after Map2Dict' |
---|
841 | print ' key / before / after' |
---|
842 | for key in sorted(parmdict.keys()): |
---|
843 | print ' '+key,'\t',before.get(key),'\t',parmdict[key] |
---|
844 | |
---|
845 | before = parmdict.copy() |
---|
846 | Dict2Map(parmdict,[]) |
---|
847 | print 'after Dict2Map' |
---|
848 | print ' key / before / after' |
---|
849 | for key in sorted(parmdict.keys()): |
---|
850 | print ' '+key,'\t',before.get(key),'\t',parmdict[key] |
---|