source: trunk/GSASIImapvars.py @ 421

Last change on this file since 421 was 421, checked in by toby, 10 years ago

determine uncertainties on constrained variables

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