source: branch/MPbranch/GSASIImapvars.py

Last change on this file was 484, checked in by vondreele, 10 years ago

change authorship
split GSASIIelem into a GUI & non-GUI parts
replace MsgDialogs? in GSASIIElem.py with simple prints
cleanup & refactor distance/angle/torsion calcs.

File size: 34.5 KB
Line 
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
9and parameter simplification contraints.
10
11Parameter redefinition is done by creating one or more relationships
12between a set of parameters
13   Mx1 * Px + My1 * Py +...
14   Mx2 * Px + Mz2 * Pz + ...
15where Pj is a parameter name and Mjk is a constant.
16
17Relations are typically supplied as input to InputParse, where either
18of 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
26Relations can also be supplied in the form of an equation:
27  nx1 * Px + ny1 * Py +... = C1
28where Cn is a constant. These equations define an algebraic
29constrant. The keyword VARY makes no sense when used with a constraint
30equation, but VARYFREE can be used.  #RVD??? is VARYFREE required or default??
31
32Unconstrained relations describe a new, independent, parameter, which
33is defined in terms of dependent parameters that are defined in the
34Model, while constrained relations effectively reduce the complexity
35of the Model by removing a degree of freedom.
36
37Relationships are grouped so that a set of dependent parameters appear
38in only one group (done in routine GroupConstraints.) Note that if a
39group contains relationships/equations that involve N dependent
40parameters, there must exist N-C new parameters, where C is the number
41of contraint equations in the group. Routine GenerateConstraints takes
42the output from InputParse and GroupConstraints generates the
43"missing" relationships and saves that information in the module's
44global variables. Each generated parameter is named
45sequentially using paramPrefix.
46
47Parameter redefinition is done by equating one (independent) parameter to several
48(now dependent) parameters, in algebraic form:
49  P1 = n2 * P2 = n3 * P3 ...
50Each equality in the relationship reduces the complexity of the model
51by one potentially variable parameter. Input is provided to routine
52StoreEquivalence in the form of an independent parameter and a list of
53dependent parameters, optionally with a multiplier.
54
55Note that none of the dependent parameters in any constraint or
56reformulation can be refined (see dependentParmList, below).
57
58External 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
85Global 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
110import numpy as np
111import re
112# data used for constraints;
113debug = 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
116dependentParmList = [] # contains a list of parameters in each group
117# note that parameters listed in dependentParmList should not be refined
118arrayList = [] # a list of of relationship matrices
119invarrayList = [] # a list of inverse relationship matrices
120indParmList = [] # a list of names for the new parameters
121fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations
122               # key is original ascii string, value is float
123fixedVarList = [] # List of variables that should not be refined
124
125# compile regular expressions used for parsing input
126rex_mult = re.compile('[+-]?[0-9.]+[eE]?[+-]?[0-9]*')
127rex_star = re.compile('\s*\*\s*')
128rex_var = re.compile('[a-zA-Z0-9_:]+')
129rex_plusminus = re.compile('\s*[+-=]\s*')
130rex_vary = re.compile('\s*Vary\s*', re.IGNORECASE)
131rex_varyfree = re.compile('(.*)\s*VaryFree\s*', re.IGNORECASE)
132
133# prefix for parameter names
134paramPrefix = "::constr:"
135consNum = 0 # number of the next constraint to be created
136
137def 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
148def 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
275def 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
309def 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
366def 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
397def 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
433def 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
442def 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
453def 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
487def 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       
509def 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
557def 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
572def 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
605def 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)
628def 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
684def 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
709def 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
729def 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
755def 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
776def 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
786if __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]
Note: See TracBrowser for help on using the repository browser.