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