source: trunk/GSASIImapvars.py @ 418

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

fix compute derivs for constrained vars

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