source: trunk/GSASIImapvars.py @ 417

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

compute derivs for constrained vars

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