source: trunk/GSASIImapvars.py @ 415

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

implement refinement of constrained parameters; fix minor bugs; add fast interpolate

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