source: trunk/GSASIImapvars.py @ 420

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

small fixes including
move msg = in SetVaryFlags? outside loop
arrange so plot shows Yo-Yc (was other way around)
put Yo-Yc at bottom in multiplot plots
more on seq refinement plots

File size: 33.5 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 varlist,mapvars,multarr,invmultarr in zip(
458    #        dependentParmList,indParmList,arrayList,invarrayList):
459    for i,name in enumerate(mapvars):
460        if name in fixedDict: continue
461        if PrintAll or name in varyList:
462            sig = sigDict.get(name)
463            printlist.append([name,parmDict[name],sig])
464    if len(printlist) == 0: return
465    s1 = ''
466    print 130*'-'
467    print "Variables generated by constraints"
468    printlist.append(3*[None])
469    for name,val,esd in printlist:
470        if len(s1) > 40 or name is None:
471            print
472            print s1
473            print s2
474            print s3
475            s1 = ''
476            if name is None: break
477        if s1 == "":
478            s1 = ' name  :'
479            s2 = ' value :'
480            s3 = ' sig   :'
481        s1 += '%12s' % (name)
482        s2 += '%12.6f' % (val)
483        if esd is None:
484            s3 += '%12s' % ('n/a')
485        else:   
486            s3 += '%12.6f' % (esd)
487           
488def VarRemapShow(varyList):
489    '''List out the saved relationships.
490    Returns a string containing the details.
491    '''
492    s = ''
493    if len(fixedVarList) > 0:
494        s += 'Fixed Variables:\n'
495        for v in fixedVarList:
496            s += '    ' + v + '\n'
497    s += 'Variable mapping relations:\n'
498    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
499    for varlist,mapvars,multarr,invmultarr in zip(
500        dependentParmList,indParmList,arrayList,invarrayList):
501        for i,mv in enumerate(mapvars):
502            if multarr is None:
503                s += '  ' + str(mv) + ' defines parameter(s): '
504                j = 0
505                for v,m in zip(varlist,invmultarr):
506                    print v,m[0]
507                    if j > 0: s += '  & '
508                    j += 1
509                    s += str(v)
510                    if m != 1:
511                        s += " / " + str(m[0])                       
512                s += '\n'
513                continue
514            s += %s = ' % mv
515            j = 0
516            for m,v in zip(multarr[i,:],varlist):
517                if m == 0: continue
518                if j > 0: s += ' + '
519                j += 1
520                s += '(%s * %s)' % (m,v)
521            if mv in varyList: s += ' VARY'
522            s += '\n'
523    s += 'Inverse variable mapping relations:\n'
524    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
525        for i,mv in enumerate(varlist):
526            s += %s = ' % mv
527            j = 0
528            for m,v in zip(invmultarr[i,:],mapvars):
529                if m == 0: continue
530                if j > 0: s += ' + '
531                j += 1
532                s += '(%s * %s)' % (m,v)
533            s += '\n'
534    return s
535
536def Dict2Deriv(varyList,derivDict,dMdv):
537    '''Compute derivatives for Independent Parameters from the
538    derivatives for the original parameters
539    '''
540    global dependentParmList,arrayList,invarrayList,indParmList,invarrayList
541    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
542        for i,name in enumerate(mapvars):
543            # grouped variables: need to add in the derv. w/r
544            # dependent variables to the independent ones
545            if name not in varyList: continue # if independent var not varied
546            for m,v in zip(invmultarr[:,i],varlist):
547                print 'add derv',v,'*',m,'to derv',name
548                if m == 0: continue
549                dMdv[varyList.index(name)] += m * derivDict[v]
550
551def Map2Dict(parmDict,varyList):
552    '''Create (or update) the Independent Parameters from the original
553    set of Parameters
554
555    Removes dependent variables from the varyList
556
557    This should be done once, before any variable refinement is done
558    to complete the parameter dictionary with the Independent Parameters
559    '''
560    # process the Independent vars: remove dependent ones from varylist
561    # and then compute values for the Independent ones from their dependents
562    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
563    for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
564        for item in varlist:
565            try:
566                varyList.remove(item)
567            except ValueError:
568                pass
569        if multarr is None: continue
570        valuelist = [parmDict[var] for var in varlist]
571        parmDict.update(zip(mapvars,
572                            np.dot(multarr,np.array(valuelist)))
573                        )
574    # now remove fixed variables from the varyList
575    global fixedVarList
576    for item in fixedVarList:
577        try:
578            varyList.remove(item)
579        except ValueError:
580            pass
581    # Set constrained parameters to their fixed values
582    parmDict.update(fixedDict)
583
584def Dict2Map(parmDict,varyList):
585    #I think this needs fixing to update parmDict with new values
586    #   from the constraints based on what is in varyList - RVD
587    '''Convert the remapped values back to the original parameters
588   
589    This should be done to apply constraints to parameter values (use
590    Map2Dict first). It should be done as part of the Model function
591    evaluation, before any computation is done
592    '''
593    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
594    # reset fixed values (should not be needed, but very quick)
595    # - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended
596    parmDict.update(fixedDict)
597    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
598        if invmultarr is None: continue
599        valuelist = [parmDict[var] for var in mapvars]
600        parmDict.update(zip(varlist,
601                            np.dot(invmultarr,np.array(valuelist)))
602                        )
603
604#======================================================================
605# internal routines follow (these routines are unlikely to be called
606# from outside the module)
607def FillInMissingRelations(arrayin,nvars):
608    '''Fill in unspecified rows in array with non-colinear unit vectors
609    arrayin is a square array, where only the first nvars rows are defined.
610   
611    Returns a new array where all rows are defined
612
613    Throws an exception if there is no non-signular result
614    (likely because two or more input rows are co-linear)
615    '''
616    a = arrayin.copy()
617    n = nvars
618    # fill in the matrix with basis vectors not colinear with any of the starting set
619    for i in range(n,len(a)):
620        try:
621            a[n] = VectorProduct(a[:n])
622        except:
623            raise Exception,"VectorProduct failed, singular input?"
624        n += 1
625    # use the G-S algorithm to compute a complete set of unit vectors orthogonal
626    # to the first in array
627    gs = GramSchmidtOrtho(a) 
628    n = nvars
629    # now replace the created vectors with the ones least correlated to the
630    # initial set
631    vlist = [v for v in gs[1:]] # drop the first row
632    for j in range(n,len(a)):
633        minavg = None
634        droplist = []
635        for k in range(len(vlist)):
636            v = vlist[k]
637            avgcor = 0
638            for i in range(n):
639                cor = np.dot(a[i],vlist[k])/(np.linalg.norm(a[i]) * np.linalg.norm(vlist[k]) )
640                if np.allclose(cor, 1.0):
641                    droplist.append(k) # found a vector co-linear to one we already have
642                                       # don't need it again
643                    #vlist.pop(k)
644                    break 
645                avgcor += cor
646            else:
647                if minavg == None:
648                    minavg = abs(avgcor/n)
649                    minnum = k
650                elif abs(avgcor/n) < minavg:
651                    minavg = abs(avgcor/n)
652                    minnum = k
653        if minavg == None:
654            raise Exception,"Failed to find a non-colinear G-S vector for row %d. Should not happen!" % n
655        a[j] = vlist[minnum]
656        droplist.append(minnum) # don't need to use this vector again
657        for i in sorted(droplist,reverse=True):
658            vlist.pop(i) 
659        n += 1
660    return a
661
662
663def MakeArray(constrDict, group, varlist):
664    """Convert the multipliers in a constraint group to an array of
665    coefficients and place in a square numpy array, adding empty rows as
666    needed.
667
668    Throws an exception if some sanity checks fail.
669    """
670    # do some error checks
671    if len(varlist) < len(group): # too many relationships -- no can do
672        raise Exception,'The number of relationships (%d) exceeds the number of parameters (%d):\n\t%s\n\t%s'% (
673            len(group),len(varlist),group,varlist)
674    # put the coefficients into an array
675    multarr = np.zeros(2*[len(varlist),])
676    i = 0
677    for cnum in group:
678        cdict = constrDict[cnum]
679        j = 0
680        for var in varlist:
681            m = cdict.get(var)
682            if m != None:
683                multarr[i,j] = m
684            j += 1
685        i += 1
686    return multarr
687
688def GramSchmidtOrtho(arrayin):
689    '''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to
690    find orthonormal unit vectors relative to first row.
691    input:
692       arrayin: a 2-D non-signular square array
693    returns:
694       a orthonormal set of unit vectors as a square array
695    '''
696    def proj(a,b):
697        'Projection operator'
698        return a*(np.dot(a,b)/np.dot(a,a))
699    a = arrayin.copy()
700    for j in range (len(a)):
701        for i in range(j):
702            a[j] = a[j] - proj(a[i],a[j])
703        if np.allclose(np.linalg.norm(a[j]),0.0):
704            raise Exception,"Singular input to GramSchmidtOrtho"
705        a[j] = a[j]/np.linalg.norm(a[j])
706    return a
707
708def VectorProduct(vl):
709    '''Evaluate the n-dimensional "cross" product of the list of vectors in vl
710    vl can be any length. The length of each vector is arbitrary, but
711    all must be the same length. See http://en.wikipedia.org/wiki/Levi-Civita_symbol
712
713    This appears to return a vector perpendicular to the supplied set.
714
715    Throws an exception if a vector can not be obtained because the input set of
716    vectors is already complete or contains any redundancies.
717   
718    Uses LeviCitvitaMult recursively to obtain all permutations of vector elements
719    '''
720    i = 0
721    newvec = []
722    for e in vl[0]:
723        i += 1
724        tl = [([i,],1),]
725        seq = LeviCitvitaMult(tl ,vl)
726        tsum = 0
727        for terms,prod in seq:
728            tsum += EvalLCterm(terms) * prod
729        newvec.append(tsum)
730    if np.allclose(newvec,0.0):
731        raise Exception,"VectorProduct failed, singular or already complete input"
732    return newvec
733
734def LeviCitvitaMult(tin,vl):
735    '''Recursion formula to compute all permutations of vector element products
736    The first term in each list is the indices of the Levi-Civita symbol, note
737    that if an index is repeated, the value is zero, so the element is not included.
738    The second term is the product of the vector elements
739    '''
740    v = vl[0]
741    vl1 = vl[1:]       
742
743    j = 0
744    tl = []
745    for e in vl[0]:
746        j += 1
747        for ind,vals in tin:
748            if j in ind: continue
749            tl.append((ind+[j,],vals*e))
750    if len(vl1):
751        return LeviCitvitaMult(tl,vl1)
752    else:
753        return tl
754
755def EvalLCterm(term):
756    '''Evaluate the Levi-Civita symbol Epsilon(i,j,k,...)'''
757    p = 1
758    for i in range(len(term)):
759        for j in range(i+1,len(term)):
760            p *= (term[j] - term[i])
761            if p == 0: return 0
762    return p/abs(p)
763
764
765if __name__ == "__main__":
766    import sys
767    #remap = MapParameters() # create an object (perhaps better as a module)
768    #remap.debug = 1
769    #debug = 1
770    parmdict = {}
771    StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), ))
772    varylist = ['2::atomx:3',]
773    #print VarRemapShow(varylist)
774    #msg = SetVaryFlags(varylist)
775    #if msg != "": print msg
776    varylist = ['2::atomx:3', '2::atomy:3', '2::atomz:3',]
777    print VarRemapShow(varylist)
778    msg = SetVaryFlags(varylist)
779    if msg != "": print msg
780    #parmdict = {
781    #    '2::atomx:3':.1 ,
782    #    '2::atomy:3':.2 , # conflicts with constraint
783    #    '2::atomz:3':.3 ,
784    #}
785
786    varylist = []
787    mapList1 = [
788        "1 * 2::atomz:3 = 0",
789        "1*p1 + 2e0*p2 - 1.0*p3",
790        "1*p4 + 1*p7",
791        "1*p1+2*p2-3.0*p2 VARY",
792        ]
793    mapList = [
794        "1 * 2::atomz:3 = 0",
795        "1*p1 + 2e0*p2 - 1.0*p3",
796        "1*p4 + 1*p7",
797        "1*p1+2*p2-3.0*p2 VARY",
798        "1*p21 + 0*p22 - 0*p23 + 0*p24 varyfree",
799        "1*p5 + 2*p6 = 1 varyfree",
800        "-1 * p6 + 1*p5",
801        "-10e-1 * p1 - -2*p2 + 3.0*p4",
802        ]
803    constrDict,constrFlag,fixedList = InputParse(mapList1)
804    #print constrDict
805    #print constrFlag
806    #print fixedList
807    groups,parmlist = GroupConstraints(constrDict)
808    GenerateConstraints(groups,parmlist,varylist,constrDict,constrFlag,fixedList)
809    print VarRemapShow(varylist)
810    sys.exit()
811    parmdict.update( {'p1':1,'p2':2,'p3':3,'p4':4,
812                      'p6':6,'p5':5,  # conflicts with constraint
813                      'p7':7,
814                      "p21":.1,"p22":.2,"p23":.3,"p24":.4,
815                      })
816    #print 'parmdict start',parmdict
817    before = parmdict.copy()
818    Map2Dict(parmdict,[])
819    print 'parmdict before and after Map2Dict'
820    print '  key / before / after'
821    for key in sorted(parmdict.keys()):
822        print '  '+key,'\t',before.get(key),'\t',parmdict[key]
823
824    before = parmdict.copy()
825    Dict2Map(parmdict,[])
826    print 'after Dict2Map'
827    print '  key / before / after'
828    for key in sorted(parmdict.keys()):
829        print '  '+key,'\t',before.get(key),'\t',parmdict[key]
Note: See TracBrowser for help on using the repository browser.