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 | varyList: a list of parameters that have been flagged to be varied. |
---|
93 | |
---|
94 | arrayList: a list by group of relationship matrices to relate |
---|
95 | parameters in dependentParmList to those in indParmList. Unlikely |
---|
96 | to be used externally. |
---|
97 | invarrayList: a list by group of relationship matrices to relate |
---|
98 | parameters in indParmList to those in dependentParmList. Unlikely |
---|
99 | to be used externally. |
---|
100 | fixedDict: a dictionary containing the fixed values corresponding |
---|
101 | to parameter equations. The dict key is an ascii string, but the |
---|
102 | dict value is a float. Unlikely to be used externally. |
---|
103 | """ |
---|
104 | |
---|
105 | import numpy as np |
---|
106 | import re |
---|
107 | # data used for constraints; |
---|
108 | debug = False # turns on printing as constraint input is processed |
---|
109 | # note that constraints are stored listed by contraint groups, where each constraint |
---|
110 | # group contains those parameters that must be handled together |
---|
111 | dependentParmList = [] # contains a list of parameters in each group |
---|
112 | # note that parameters listed in dependentParmList should not be refined |
---|
113 | arrayList = [] # a list of of relationship matrices |
---|
114 | invarrayList = [] # a list of inverse relationship matrices |
---|
115 | indParmList = [] # a list of names for the new parameters |
---|
116 | fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations |
---|
117 | # key is original ascii string, value is float |
---|
118 | varyList = [] # a list of varied constraints |
---|
119 | |
---|
120 | # compile regular expressions used for parsing input |
---|
121 | rex_mult = re.compile('[+-]?[0-9.]+[eE]?[+-]?[0-9]*') |
---|
122 | rex_star = re.compile('\s*\*\s*') |
---|
123 | rex_var = re.compile('\w+') |
---|
124 | rex_plusminus = re.compile('\s*[+-=]\s*') |
---|
125 | rex_vary = re.compile('\s*Vary\s*', re.IGNORECASE) |
---|
126 | rex_varyfree = re.compile('(.*)\s*VaryFree\s*', re.IGNORECASE) |
---|
127 | |
---|
128 | # prefix for parameter names |
---|
129 | paramPrefix = "::constr:" |
---|
130 | consNum = 0 # number of the next constraint to be created |
---|
131 | |
---|
132 | def InitVars(): |
---|
133 | '''Initializes all constraint information''' |
---|
134 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,varyList,consNum |
---|
135 | dependentParmList = [] # contains a list of parameters in each group |
---|
136 | arrayList = [] # a list of of relationship matrices |
---|
137 | invarrayList = [] # a list of inverse relationship matrices |
---|
138 | indParmList = [] # a list of names for the new parameters |
---|
139 | fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations |
---|
140 | varyList = [] # a list of varied constraints |
---|
141 | consNum = 0 # number of the next constraint to be created |
---|
142 | |
---|
143 | def InputParse(mapList): |
---|
144 | '''Converts a set relationships used to remap parameters into composite |
---|
145 | parameters or to one or more constants. |
---|
146 | |
---|
147 | Input: |
---|
148 | mapList: a list of character strings where each string defines |
---|
149 | a relationship in form: |
---|
150 | ('<const11>*<prm1> [+-] <const12>*<prm2> [+-] ... = <value> [VaryFree]', |
---|
151 | '<const21>*<prm1> [+-] <const22>*<prm2> [+-] ... [Vary/VaryFree]', |
---|
152 | ...) |
---|
153 | these define either a constraint or a new independent parameter, |
---|
154 | where <constXX> is a constant containing an optional sign, digits, |
---|
155 | optionally a decimal point and/or an exponent, prefixed by e or E, |
---|
156 | and <prmN> is the name of a parameter defined in the Model. |
---|
157 | Parameters can be included in any order. Multiplying a parameter |
---|
158 | by 0 causes that parameter to be included in a group of |
---|
159 | relationships. Note that if the same parameter is included twice |
---|
160 | in a relationship/equation, the coefficients are added. |
---|
161 | |
---|
162 | When an equality is defined, a constant, <value>, is |
---|
163 | included. This then describes a constraint equation. |
---|
164 | |
---|
165 | Vary is an optional keyword that indicates the indicated remapped |
---|
166 | parameter should be varied (note, case insensitive). Should be the |
---|
167 | last item on the line. Makes no sense with an equation. |
---|
168 | |
---|
169 | VaryFree is an optional keyword that indicates all possible |
---|
170 | remapped parameters should be varied (note, case |
---|
171 | insensitive). Makes most sense with a constraint equation. Should |
---|
172 | be the last item on the line. |
---|
173 | |
---|
174 | returns |
---|
175 | constrDict: a list with a dict for each item in mapList that |
---|
176 | defines the relationship. The keys are parameter names and the |
---|
177 | values are the multiplier for the parameter name |
---|
178 | constrFlag: a list for each item in mapList with a list contains |
---|
179 | 'Vary' and/or 'VaryFree' |
---|
180 | fixedList: a list for each item in mapList. Contains the value |
---|
181 | (as a string) for each contraint equation or None for a |
---|
182 | constraint relationship. |
---|
183 | ''' |
---|
184 | i = 0 |
---|
185 | constrDict = [] |
---|
186 | fixedList = [] |
---|
187 | constrFlag = [] |
---|
188 | if debug: print 50*'-','\n(Input)' |
---|
189 | for line in mapList: |
---|
190 | line = line.strip() |
---|
191 | if len(line) == 0: continue # empty lines are ignored |
---|
192 | constrDict.append({}) |
---|
193 | constrFlag.append([]) |
---|
194 | fixedList.append(None) |
---|
195 | i += 1 |
---|
196 | fixedval = None |
---|
197 | if debug: print '\t',line |
---|
198 | try: |
---|
199 | j = 0 |
---|
200 | sign = 1.0 |
---|
201 | while len(line): |
---|
202 | j += 1 |
---|
203 | m = sign * float(rex_mult.match(line).group()) |
---|
204 | line = line[rex_mult.match(line).end():] |
---|
205 | j += 1 |
---|
206 | line = line[rex_star.match(line).end():] |
---|
207 | j += 1 |
---|
208 | v = rex_var.match(line).group() |
---|
209 | line = line[rex_var.match(line).end():] |
---|
210 | #print m,'times',v |
---|
211 | #if v not in varlist: varlist.append(v) |
---|
212 | if constrDict[i-1].get(v) == None: |
---|
213 | constrDict[i-1][v] = m |
---|
214 | else: |
---|
215 | constrDict[i-1][v] += m |
---|
216 | if len(line.strip()) == 0: break |
---|
217 | j += 1 |
---|
218 | # advance to next separator (+-=) |
---|
219 | pm = rex_plusminus.match(line) |
---|
220 | if pm != None: |
---|
221 | line = line[pm.end():] |
---|
222 | pm = pm.group() |
---|
223 | else: |
---|
224 | pm = '' |
---|
225 | if pm.strip() == '+': |
---|
226 | sign = 1.0 |
---|
227 | elif pm.strip() == '-': |
---|
228 | sign = -1.0 |
---|
229 | elif pm.strip() == '=': |
---|
230 | # found a fixed value, also check for a VaryFree flag |
---|
231 | if rex_varyfree.match(line): |
---|
232 | constrFlag[-1].append('VaryFree') |
---|
233 | #fixedval = float(rex_varyfree.split(line)[1]) |
---|
234 | fixedval = rex_varyfree.split(line)[1].strip() |
---|
235 | else: |
---|
236 | #fixedval = float(line.strip()) |
---|
237 | fixedval = line.strip() |
---|
238 | fixedList[-1] = fixedval |
---|
239 | line = "" |
---|
240 | elif rex_varyfree.match(line): |
---|
241 | constrFlag[-1].append('VaryFree') |
---|
242 | line = line[rex_varyfree.match(line).end():] |
---|
243 | elif rex_vary.match(line): |
---|
244 | constrFlag[-1].append('Vary') |
---|
245 | line = line[rex_vary.match(line).end():] |
---|
246 | else: |
---|
247 | # something else is on the line, but not a keyword |
---|
248 | raise Exception |
---|
249 | except SyntaxError: |
---|
250 | if debug: print 'Error in line',i,'token',j,'@','"'+line+'"' |
---|
251 | raise Exception,'Error in line %d token %d, beginning with %s'% (i,j, line) |
---|
252 | |
---|
253 | if debug: # on debug, show what is parsed in a semi-readable |
---|
254 | print 50*'-','\n(parsed relationship/equation & flag)' |
---|
255 | for i in range(len(constrDict)): |
---|
256 | flags = '' |
---|
257 | for f in constrFlag[i]: |
---|
258 | if flags != '': |
---|
259 | flags = flags + ' + ' + f |
---|
260 | else: |
---|
261 | flags = f |
---|
262 | if fixedList[i] is None: |
---|
263 | print '#',i+1,constrDict[i],'\t',constrFlag[i] |
---|
264 | else: |
---|
265 | print '#',i+1,constrDict[i],'=',fixedList[i],'\t',constrFlag[i] |
---|
266 | |
---|
267 | return constrDict,constrFlag,fixedList |
---|
268 | |
---|
269 | def GroupConstraints(constrDict): |
---|
270 | """divide the constraints into groups that share no parameters. |
---|
271 | Input |
---|
272 | constrDict: a list of dicts defining relationships/constraints |
---|
273 | (see InputParse) |
---|
274 | Returns two lists of lists: |
---|
275 | a list of relationship list indices for each group and |
---|
276 | a list containing the parameters used in each group |
---|
277 | """ |
---|
278 | assignedlist = [] # relationships that have been used |
---|
279 | groups = [] # contains a list of grouplists |
---|
280 | ParmList = [] |
---|
281 | for i,consi in enumerate(constrDict): |
---|
282 | if i in assignedlist: continue # already in a group, skip |
---|
283 | # starting a new group |
---|
284 | grouplist = [i,] |
---|
285 | assignedlist.append(i) |
---|
286 | groupset = set(consi.keys()) |
---|
287 | changes = True # always loop at least once |
---|
288 | while(changes): # loop until we can't find anything to add to the current group |
---|
289 | changes = False # but don't loop again unless we find something |
---|
290 | for j,consj in enumerate(constrDict): |
---|
291 | if j in assignedlist: continue # already in a group, skip |
---|
292 | if len(set(consj.keys()) & groupset) > 0: # true if this needs to be added |
---|
293 | changes = True |
---|
294 | grouplist.append(j) |
---|
295 | assignedlist.append(j) |
---|
296 | groupset = groupset | set(consj.keys()) |
---|
297 | group = sorted(grouplist) |
---|
298 | varlist = sorted(list(groupset)) |
---|
299 | groups.append(group) |
---|
300 | ParmList.append(varlist) |
---|
301 | return groups,ParmList |
---|
302 | |
---|
303 | def GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList): |
---|
304 | '''Takes a list of relationship entries comprising a group of constraints and |
---|
305 | builds the relationship lists and their inverse and stores them in global variables |
---|
306 | ''' |
---|
307 | global dependentParmList,arrayList,invarrayList,indParmList,varyList,consNum |
---|
308 | for group,varlist in zip(groups,parmlist): |
---|
309 | VaryFree = False |
---|
310 | for row in group: |
---|
311 | if 'VaryFree' in constrFlag[row]: VaryFree = True |
---|
312 | arr = MakeArray(constrDict, group, varlist) |
---|
313 | constrArr = FillInMissingRelations(arr,len(group)) |
---|
314 | mapvar = [] |
---|
315 | group = group[:] |
---|
316 | for i in range(len(varlist)): |
---|
317 | vary = VaryFree |
---|
318 | if len(group) > 0: |
---|
319 | rel = group.pop(0) |
---|
320 | fixedval = fixedList[rel] |
---|
321 | if 'Vary' in constrFlag[rel]: vary = True |
---|
322 | else: |
---|
323 | fixedval = None |
---|
324 | if fixedval is None: |
---|
325 | varname = paramPrefix + str(consNum) |
---|
326 | mapvar.append(varname) |
---|
327 | consNum += 1 |
---|
328 | if vary: varyList.append(varname) |
---|
329 | else: |
---|
330 | mapvar.append(fixedval) |
---|
331 | dependentParmList.append(varlist) |
---|
332 | arrayList.append(constrArr) |
---|
333 | invarrayList.append(np.linalg.inv(constrArr)) |
---|
334 | indParmList.append(mapvar) |
---|
335 | # setup dictionary containing the fixed values |
---|
336 | global fixedDict |
---|
337 | # key is original ascii string, value is float |
---|
338 | for fixedval in fixedList: |
---|
339 | if fixedval is not None: |
---|
340 | fixedDict[fixedval] = float(fixedval) |
---|
341 | |
---|
342 | if debug: # on debug, show what is parsed & generated, semi-readable |
---|
343 | print 50*'-' |
---|
344 | for group,varlist,multarr,inv,mapvar in zip(groups,parmlist,arrayList,invarrayList,indParmList): |
---|
345 | print '\n*** relation(s) in group:',group,'\tvars in group:',varlist |
---|
346 | print 'new parameters:', mapvar |
---|
347 | print 'Input relationship matrix' |
---|
348 | print multarr[:len(group)] |
---|
349 | added = len(group) - len(varlist) |
---|
350 | if added < 0: |
---|
351 | print 'added relationship rows' |
---|
352 | print multarr[added:] |
---|
353 | print 'Inverse relationship matrix' |
---|
354 | print inv |
---|
355 | |
---|
356 | def StoreEquivalence(independentVar,dependentList): |
---|
357 | '''Takes a list of dependent parameter(s) and stores their |
---|
358 | relationship to a single independent parameter (independentVar) |
---|
359 | input: |
---|
360 | independentVar -- name of parameter that will set others (may |
---|
361 | be varied) |
---|
362 | dependentList -- list of parameter that will set from |
---|
363 | independentVar (may not be varied). Each item can be a parameter |
---|
364 | name or a tuple containing a name and multiplier: |
---|
365 | ('parm1',('parm2',.5),) |
---|
366 | ''' |
---|
367 | |
---|
368 | global dependentParmList,arrayList,invarrayList,indParmList |
---|
369 | mapList = [] |
---|
370 | multlist = [] |
---|
371 | for var in dependentList: |
---|
372 | if isinstance(var, basestring): |
---|
373 | mult = 1.0 |
---|
374 | elif len(var) == 2: |
---|
375 | var,mult = var |
---|
376 | else: |
---|
377 | raise Exception, "Cannot parse "+repr(var) + " as var or (var,multiplier)" |
---|
378 | mapList.append(var) |
---|
379 | multlist.append(tuple((mult,))) |
---|
380 | # added relationships to stored values |
---|
381 | arrayList.append(None) |
---|
382 | invarrayList.append(np.array(multlist)) |
---|
383 | indParmList.append(tuple((independentVar,))) |
---|
384 | dependentParmList.append(mapList) |
---|
385 | return |
---|
386 | |
---|
387 | def VarRemapShow(): |
---|
388 | '''List out the saved relationships. |
---|
389 | Returns a string containing the details. |
---|
390 | ''' |
---|
391 | s = 'Mapping relations:\n' |
---|
392 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,varyList |
---|
393 | for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList): |
---|
394 | i = 0 |
---|
395 | for mv in mapvars: |
---|
396 | if multarr is None: |
---|
397 | s += ' ' + str(mv) + ' defines parameter(s): ' |
---|
398 | j = 0 |
---|
399 | for v in varlist: |
---|
400 | if j > 0: s += ' & ' |
---|
401 | j += 1 |
---|
402 | s += str(v) |
---|
403 | s += '\n' |
---|
404 | continue |
---|
405 | s += ' %s = ' % mv |
---|
406 | j = 0 |
---|
407 | for m,v in zip(multarr[i,:],varlist): |
---|
408 | if m == 0: continue |
---|
409 | if j > 0: s += ' + ' |
---|
410 | j += 1 |
---|
411 | s += '(%s * %s)' % (m,v) |
---|
412 | if mv in varyList: s += ' VARY' |
---|
413 | s += '\n' |
---|
414 | i += 1 |
---|
415 | s += 'Inverse mapping relations:\n' |
---|
416 | for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList): |
---|
417 | i = 0 |
---|
418 | for mv in varlist: |
---|
419 | s += ' %s = ' % mv |
---|
420 | j = 0 |
---|
421 | for m,v in zip(invmultarr[i,:],mapvars): |
---|
422 | if m == 0: continue |
---|
423 | if j > 0: s += ' + ' |
---|
424 | j += 1 |
---|
425 | s += '(%s * %s)' % (m,v) |
---|
426 | s += '\n' |
---|
427 | i += 1 |
---|
428 | return s |
---|
429 | |
---|
430 | def Map2Dict(parmDict,varyList): |
---|
431 | '''Create (or update) the Independent Parameters from the original |
---|
432 | set of Parameters |
---|
433 | |
---|
434 | This should be done once, before any variable refinement is done |
---|
435 | to complete the parameter dictionary with the Independent Parameters |
---|
436 | ''' |
---|
437 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict |
---|
438 | for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList): |
---|
439 | for item in varlist: |
---|
440 | try: |
---|
441 | varyList.remove(item) |
---|
442 | except ValueError: |
---|
443 | pass |
---|
444 | if multarr is None: continue |
---|
445 | valuelist = [parmdict[var] for var in varlist] |
---|
446 | parmDict.update(zip(mapvars, |
---|
447 | np.dot(multarr,np.array(valuelist))) |
---|
448 | ) |
---|
449 | # overwrite dict with constraints - why not parmDict.update(fixDict)? |
---|
450 | parmDict.update(fixedDict) |
---|
451 | |
---|
452 | def Dict2Map(parmDict): |
---|
453 | '''Convert the remapped values back to the original parameters |
---|
454 | |
---|
455 | This should be done to apply constraints to parameter values (use |
---|
456 | Map2Dict first). It should be done as part of the Model function |
---|
457 | evaluation, before any computation is done |
---|
458 | ''' |
---|
459 | global dependentParmList,arrayList,invarrayList,indParmList,fixedDict |
---|
460 | # reset fixed values (should not be needed, but very quick) |
---|
461 | parmDict.update(fixedDict) |
---|
462 | for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList): |
---|
463 | if invmultarr is None: continue |
---|
464 | valuelist = [parmDict[var] for var in mapvars] |
---|
465 | parmDict.update(zip(varlist, |
---|
466 | np.dot(invmultarr,np.array(valuelist))) |
---|
467 | ) |
---|
468 | |
---|
469 | #====================================================================== |
---|
470 | # internal routines follow (these routines are unlikely to be called |
---|
471 | # from outside the module) |
---|
472 | def FillInMissingRelations(arrayin,nvars): |
---|
473 | '''Fill in unspecified rows in array with non-colinear unit vectors |
---|
474 | arrayin is a square array, where only the first nvars rows are defined. |
---|
475 | |
---|
476 | Returns a new array where all rows are defined |
---|
477 | |
---|
478 | Throws an exception if there is no non-signular result |
---|
479 | (likely because two or more input rows are co-linear) |
---|
480 | ''' |
---|
481 | a = arrayin.copy() |
---|
482 | n = nvars |
---|
483 | # fill in the matrix with basis vectors not colinear with any of the starting set |
---|
484 | for i in range(n,len(a)): |
---|
485 | try: |
---|
486 | a[n] = VectorProduct(a[:n]) |
---|
487 | except: |
---|
488 | raise Exception,"VectorProduct failed, singular input?" |
---|
489 | n += 1 |
---|
490 | # use the G-S algorithm to compute a complete set of unit vectors orthogonal |
---|
491 | # to the first in array |
---|
492 | gs = GramSchmidtOrtho(a) |
---|
493 | n = nvars |
---|
494 | # now replace the created vectors with the ones least correlated to the |
---|
495 | # initial set |
---|
496 | vlist = [v for v in gs[1:]] # drop the first row |
---|
497 | for j in range(n,len(a)): |
---|
498 | minavg = None |
---|
499 | droplist = [] |
---|
500 | for k in range(len(vlist)): |
---|
501 | v = vlist[k] |
---|
502 | avgcor = 0 |
---|
503 | for i in range(n): |
---|
504 | cor = np.dot(a[i],vlist[k])/(np.linalg.norm(a[i]) * np.linalg.norm(vlist[k]) ) |
---|
505 | if np.allclose(cor, 1.0): |
---|
506 | droplist.append(k) # found a vector co-linear to one we already have |
---|
507 | # don't need it again |
---|
508 | #vlist.pop(k) |
---|
509 | break |
---|
510 | avgcor += cor |
---|
511 | else: |
---|
512 | if minavg == None: |
---|
513 | minavg = abs(avgcor/n) |
---|
514 | minnum = k |
---|
515 | elif abs(avgcor/n) < minavg: |
---|
516 | minavg = abs(avgcor/n) |
---|
517 | minnum = k |
---|
518 | if minavg == None: |
---|
519 | raise Exception,"Failed to find a non-colinear G-S vector for row %d. Should not happen!" % n |
---|
520 | a[j] = vlist[minnum] |
---|
521 | droplist.append(minnum) # don't need to use this vector again |
---|
522 | for i in sorted(droplist,reverse=True): |
---|
523 | vlist.pop(i) |
---|
524 | n += 1 |
---|
525 | return a |
---|
526 | |
---|
527 | |
---|
528 | def MakeArray(constrDict, group, varlist): |
---|
529 | """Convert the multipliers in a constraint group to an array of |
---|
530 | coefficients and place in a square numpy array, adding empty rows as |
---|
531 | needed. |
---|
532 | |
---|
533 | Throws an exception if some sanity checks fail. |
---|
534 | """ |
---|
535 | # do some error checks |
---|
536 | if len(varlist) < len(group): # too many relationships -- no can do |
---|
537 | raise Exception,'The number of relationships (%d) exceeds the number of parameters (%d):\n\t%s\n\t%s'% ( |
---|
538 | len(group),len(varlist),group,varlist) |
---|
539 | if len(varlist) == 1: # one to one mapping -- something is probably wrong |
---|
540 | raise Exception,'Why remap a single parameter? (%s)'% (varlist[0]) |
---|
541 | # put the coefficients into an array |
---|
542 | multarr = np.zeros(2*[len(varlist),]) |
---|
543 | i = 0 |
---|
544 | for cnum in group: |
---|
545 | cdict = constrDict[cnum] |
---|
546 | j = 0 |
---|
547 | for var in varlist: |
---|
548 | m = cdict.get(var) |
---|
549 | if m != None: |
---|
550 | multarr[i,j] = m |
---|
551 | j += 1 |
---|
552 | i += 1 |
---|
553 | return multarr |
---|
554 | |
---|
555 | def GramSchmidtOrtho(arrayin): |
---|
556 | '''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to |
---|
557 | find orthonormal unit vectors relative to first row. |
---|
558 | input: |
---|
559 | arrayin: a 2-D non-signular square array |
---|
560 | returns: |
---|
561 | a orthonormal set of unit vectors as a square array |
---|
562 | ''' |
---|
563 | def proj(a,b): |
---|
564 | 'Projection operator' |
---|
565 | return a*(np.dot(a,b)/np.dot(a,a)) |
---|
566 | a = arrayin.copy() |
---|
567 | for j in range (len(a)): |
---|
568 | for i in range(j): |
---|
569 | a[j] = a[j] - proj(a[i],a[j]) |
---|
570 | if np.allclose(np.linalg.norm(a[j]),0.0): |
---|
571 | raise Exception,"Singular input to GramSchmidtOrtho" |
---|
572 | a[j] = a[j]/np.linalg.norm(a[j]) |
---|
573 | return a |
---|
574 | |
---|
575 | def VectorProduct(vl): |
---|
576 | '''Evaluate the n-dimensional "cross" product of the list of vectors in vl |
---|
577 | vl can be any length. The length of each vector is arbitrary, but |
---|
578 | all must be the same length. See http://en.wikipedia.org/wiki/Levi-Civita_symbol |
---|
579 | |
---|
580 | This appears to return a vector perpendicular to the supplied set. |
---|
581 | |
---|
582 | Throws an exception if a vector can not be obtained because the input set of |
---|
583 | vectors is already complete or contains any redundancies. |
---|
584 | |
---|
585 | Uses LeviCitvitaMult recursively to obtain all permutations of vector elements |
---|
586 | ''' |
---|
587 | i = 0 |
---|
588 | newvec = [] |
---|
589 | for e in vl[0]: |
---|
590 | i += 1 |
---|
591 | tl = [([i,],1),] |
---|
592 | seq = LeviCitvitaMult(tl ,vl) |
---|
593 | tsum = 0 |
---|
594 | for terms,prod in seq: |
---|
595 | tsum += EvalLCterm(terms) * prod |
---|
596 | newvec.append(tsum) |
---|
597 | if np.allclose(newvec,0.0): |
---|
598 | raise Exception,"VectorProduct failed, singular or already complete input" |
---|
599 | return newvec |
---|
600 | |
---|
601 | def LeviCitvitaMult(tin,vl): |
---|
602 | '''Recursion formula to compute all permutations of vector element products |
---|
603 | The first term in each list is the indices of the Levi-Civita symbol, note |
---|
604 | that if an index is repeated, the value is zero, so the element is not included. |
---|
605 | The second term is the product of the vector elements |
---|
606 | ''' |
---|
607 | v = vl[0] |
---|
608 | vl1 = vl[1:] |
---|
609 | |
---|
610 | j = 0 |
---|
611 | tl = [] |
---|
612 | for e in vl[0]: |
---|
613 | j += 1 |
---|
614 | for ind,vals in tin: |
---|
615 | if j in ind: continue |
---|
616 | tl.append((ind+[j,],vals*e)) |
---|
617 | if len(vl1): |
---|
618 | return LeviCitvitaMult(tl,vl1) |
---|
619 | else: |
---|
620 | return tl |
---|
621 | |
---|
622 | def EvalLCterm(term): |
---|
623 | '''Evaluate the Levi-Civita symbol Epsilon(i,j,k,...)''' |
---|
624 | p = 1 |
---|
625 | for i in range(len(term)): |
---|
626 | for j in range(i+1,len(term)): |
---|
627 | p *= (term[j] - term[i]) |
---|
628 | if p == 0: return 0 |
---|
629 | return p/abs(p) |
---|
630 | |
---|
631 | |
---|
632 | if __name__ == "__main__": |
---|
633 | import sys |
---|
634 | #remap = MapParameters() # create an object (perhaps better as a module) |
---|
635 | #remap.debug = 1 |
---|
636 | #debug = 1 |
---|
637 | parmdict = {} |
---|
638 | StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) |
---|
639 | print VarRemapShow() |
---|
640 | |
---|
641 | parmdict = { |
---|
642 | '2::atomx:3':.1 , |
---|
643 | '2::atomy:3':.2 , # conflicts with constraint |
---|
644 | '2::atomz:3':.3 , |
---|
645 | } |
---|
646 | |
---|
647 | mapList = [ |
---|
648 | "1*p1 + 2e0*p2 - 1.0*p3", |
---|
649 | "1*p4 + 1*p7", |
---|
650 | "1*p1+2*p2-3.0*p2 VARY", |
---|
651 | "1*p21 + 0*p22 - 0*p23 + 0*p24 varyfree", |
---|
652 | "1*p5 + 2*p6 = 1 varyfree", |
---|
653 | "-1 * p6 + 1*p5", |
---|
654 | "-10e-1 * p1 - -2*p2 + 3.0*p4", |
---|
655 | ] |
---|
656 | constrDict,constrFlag,fixedList = InputParse(mapList) |
---|
657 | print constrDict |
---|
658 | print constrFlag |
---|
659 | print fixedList |
---|
660 | groups,parmlist = GroupConstraints(constrDict) |
---|
661 | GenerateConstraints(groups,parmlist,constrDict,constrFlag,fixedList) |
---|
662 | print VarRemapShow() |
---|
663 | parmdict.update( {'p1':1,'p2':2,'p3':3,'p4':4, |
---|
664 | 'p6':6,'p5':5, # conflicts with constraint |
---|
665 | 'p7':7, |
---|
666 | "p21":.1,"p22":.2,"p23":.3,"p24":.4, |
---|
667 | }) |
---|
668 | #print 'parmdict start',parmdict |
---|
669 | before = parmdict.copy() |
---|
670 | Map2Dict(parmdict,[]) |
---|
671 | print 'parmdict before and after Map2Dict' |
---|
672 | print ' key / before / after' |
---|
673 | for key in sorted(parmdict.keys()): |
---|
674 | print ' '+key,'\t',before.get(key),'\t',parmdict[key] |
---|
675 | |
---|
676 | before = parmdict.copy() |
---|
677 | Dict2Map(parmdict) |
---|
678 | print 'after Dict2Map' |
---|
679 | print ' key / before / after' |
---|
680 | for key in sorted(parmdict.keys()): |
---|
681 | print ' '+key,'\t',before.get(key),'\t',parmdict[key] |
---|