source: trunk/GSASIImapvars.py @ 567

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

more work on constraints

File size: 37.3 KB
Line 
1########### SVN repository information ###################
2# $Date: 2011-01-07 13:27:24 -0600 (Fri, 07 Jan 2011) $
3# $Author: vondreele & toby $
4# $Revision: 230 $
5# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIImapvars.py $
6# $Id: GSASIImapvars.py 230 2011-01-07 19:27:24Z vondreele $
7########### SVN repository information ###################
8"""Module that implements algebraic contraints, parameter redefinition
9and parameter simplification contraints.
10
11Parameter redefinition (new vars) is done by creating one or more relationships
12between a set of parameters
13   Mx1 * Px + My1 * Py +...
14   Mx2 * Px + Mz2 * Pz + ...
15where Pj is a parameter name and Mjk is a constant.
16
17Constant constraint Relations can also be supplied in the form of an equation:
18  nx1 * Px + ny1 * Py +... = C1
19where Cn is a constant. These equations define an algebraic
20constrant.
21
22Parameters can also be "fixed" (held), which prevents them from being refined.
23
24All of the above three cases is supplied as input using routines
25GroupConstraints and GenerateConstraints. The input consists of a list of
26relationship dictionaries:
27    constrDict = [
28         {'0:12:Scale': 2.0, '0:14:Scale': 4.0, '0:13:Scale': 3.0, '0:0:Scale': 0.5},
29         {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
30         {'0::A0': 0.0}]
31
32    fixedList = ['5.0', None, '0']
33
34Where the dictionary defines the first part of an expression and the corresponding fixedList
35item is either None (for parameter redefinition) or the constant value, for a constant
36constraint equation. A dictionary that contains a single term defines a variable that
37will be fixed (held). The multiplier and the fixedList value in this case are ignored.
38
39Parameters can also be equivalenced or "slaved" to another parameter, such that one
40(independent) parameter is equated to several (now dependent) parameters. In
41algebraic form this is:
42   P0 = M1 * P1 = M2 * P2 = ...
43Thus parameters P0, P1 and P2,... are linearly equivalent. Routine StoreEquivalence is
44used to specify these equivalences.
45
46Parameter redefinition (new vars) describes a new, independent, parameter, which
47is defined in terms of dependent parameters that are defined in the
48Model, while fixed constrained relations effectively reduce the complexity
49of the Model by removing a degree of freedom. It is possible for a parameter to appear
50in both a parameter redefinition expression and a fixed constraint equation, but a
51parameter cannot be used a parameter equivalance cannot be used elsewhere (not fixed,
52constrained or redefined). Likewise a fixed parameter cannot be used elsewhere (not
53equivalanced, constrained or redefined).
54
55Relationships are grouped so that a set of dependent parameters appear
56in only one group (done in routine GroupConstraints.) Note that if a
57group contains relationships/equations that involve N dependent
58parameters, there must exist N-C new parameters, where C is the number
59of contraint equations in the group. Routine GenerateConstraints takes
60the output from InputParse and GroupConstraints generates the
61"missing" relationships and saves that information in the module's
62global variables. Each generated parameter is named sequentially using paramPrefix.
63
64A list of parameters that will be varied is specified as input to GenerateConstraints
65(varyList). A fixed parameter will simply be removed from this list preventing that
66parameter from being varied. Note that all parameters in a relationship must specified as
67varied (appear in varyList) or none can be varied. This is checked in GenerateConstraints
68(as well as for generated relationships in SetVaryFlags).
69* If all parameters in a parameter redefinition (new var) relationship are varied, the
70  parameter assigned to this expression (::constr:n, see paramPrefix) newly generated
71  parameter is varied. Note that any generated "missing" relations are not varied. Only
72  the input relations are varied.
73* If all parameters in a fixed constraint equation are varied, the generated "missing"
74  relations in the group are all varied. This provides the N-C degrees of freedom.
75
76External Routines:
77   To define a set of constrained and unconstrained relations, one
78     calls InputParse, GroupConstraints and GenerateConstraints. It
79     is best to supply all relations/equations in one call to these
80     routines so that grouping is done correctly.
81
82   To implement parameter redefinition, one calls
83     StoreEquivalence. This should be called for every set of
84     equivalence relationships. There is no harm in using
85     StoreEquivalence with the same independent variable:
86       StoreEquivalence('x',('y',))
87       StoreEquivalence('x',('z',))
88     (though StoreEquivalence('x',('y','z')) will run more
89     efficiently) but mixing independent and dependent variables is
90     problematic. This may not work properly:
91        StoreEquivalence('x',('y',))
92        StoreEquivalence('y',('z',))
93
94   To show a summary of the parameter remapping, one calls
95      VarRemapShow
96
97   To determine values for the parameters created in this module, one
98      calls Map2Dict. This will not apply contraints.
99
100   To take values from the new independent parameters and constraints,
101      one calls Dict2Map. This will apply contraints.
102
103Global Variables:
104   dependentParmList: contains a list by group of lists of
105     parameters used in the group. Note that parameters listed in
106     dependentParmList should not be refined as they will not affect
107     the model
108   indParmList: a list of groups of Independent parameters defined in
109     each group. This contains both parameters used in parameter
110     redefinitions as well as names of generated new parameters.
111
112   fixedVarList: a list of variables that have been 'fixed'
113     by defining them as equal to a constant (::var: = 0). Note that
114     the constant value is ignored at present. These variables are
115     later removed from varyList which prevents them from being refined.
116     Unlikely to be used externally.
117   arrayList: a list by group of relationship matrices to relate
118     parameters in dependentParmList to those in indParmList. Unlikely
119     to be used externally.
120   invarrayList: a list by group of relationship matrices to relate
121     parameters in indParmList to those in dependentParmList. Unlikely
122     to be used externally.
123   fixedDict: a dictionary containing the fixed values corresponding
124     to parameter equations.  The dict key is an ascii string, but the
125     dict value is a float.  Unlikely to be used externally.
126"""
127
128import numpy as np
129# data used for constraints;
130debug = False # turns on printing as constraint input is processed
131# note that constraints are stored listed by contraint groups, where each constraint
132# group contains those parameters that must be handled together
133dependentParmList = [] # contains a list of parameters in each group
134# note that parameters listed in dependentParmList should not be refined
135arrayList = [] # a list of of relationship matrices
136invarrayList = [] # a list of inverse relationship matrices
137indParmList = [] # a list of names for the new parameters
138fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations
139               # key is original ascii string, value is float
140fixedVarList = [] # List of variables that should not be refined
141
142# prefix for parameter names
143paramPrefix = "::constr:"
144consNum = 0 # number of the next constraint to be created
145
146def InitVars():
147    '''Initializes all constraint information'''
148    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,consNum
149    dependentParmList = [] # contains a list of parameters in each group
150    arrayList = [] # a list of of relationship matrices
151    invarrayList = [] # a list of inverse relationship matrices
152    indParmList = [] # a list of names for the new parameters
153    fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations
154    consNum = 0 # number of the next constraint to be created
155    fixedVarList = []
156
157def GroupConstraints(constrDict):
158    """divide the constraints into groups that share no parameters.
159    Input
160       constrDict: a list of dicts defining relationships/constraints
161       (see InputParse)
162    Returns two lists of lists:
163      a list of relationship list indices for each group and
164      a list containing the parameters used in each group
165      """
166    assignedlist = [] # relationships that have been used
167    groups = [] # contains a list of grouplists
168    ParmList = []
169    for i,consi in enumerate(constrDict):
170        if i in assignedlist: continue # already in a group, skip
171        # starting a new group
172        grouplist = [i,]
173        assignedlist.append(i)
174        groupset = set(consi.keys())
175        changes = True # always loop at least once
176        while(changes): # loop until we can't find anything to add to the current group
177            changes = False # but don't loop again unless we find something
178            for j,consj in enumerate(constrDict):
179                if j in assignedlist: continue # already in a group, skip
180                if len(set(consj.keys()) & groupset) > 0: # true if this needs to be added
181                    changes = True
182                    grouplist.append(j)
183                    assignedlist.append(j)
184                    groupset = groupset | set(consj.keys())
185        group = sorted(grouplist)
186        varlist = sorted(list(groupset))
187        groups.append(group)
188        ParmList.append(varlist)
189    return groups,ParmList
190
191def GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList):
192    '''Takes a list of relationship entries comprising a group of
193    constraints and builds the relationship lists and their inverse
194    and stores them in global variables Also checks for internal
195    conflicts or inconsistencies in parameter/variable definitions.
196    '''
197    global dependentParmList,arrayList,invarrayList,indParmList,consNum
198    msg = ''
199
200    # process fixed (held) variables
201    for cdict in constrDict:
202        if len(cdict) == 1:
203            fixedVarList.append(cdict.keys()[0])
204    #print 'fixedVarList',fixedVarList
205   
206    # process equivalences: make a list of dependent and independent vars
207    #    and check for repeated uses (repetition of a parameter as an
208    #    independent var is OK)
209    indepVarList = []
210    depVarList = []
211    multdepVarList = []
212    for varlist,mapvars,multarr,invmultarr in zip(
213        dependentParmList,indParmList,arrayList,invarrayList):
214        if multarr is None:
215            zeromult = False
216            for mv in mapvars:
217                varied = 0
218                notvaried = ''
219                if mv in varyList:
220                    varied += 1
221                else:
222                    if notvaried: notvaried += ', '
223                    notvaried += mv
224                if mv not in indepVarList: indepVarList.append(mv)
225                for v,m in zip(varlist,invmultarr):
226                    if m == 0: zeromult = True
227                    if v in varyList:
228                        varied += 1
229                    else:
230                        if notvaried: notvaried += ', '
231                        notvaried += v
232                    if v in depVarList:
233                        multdepVarList.append(v)
234                    else:
235                        depVarList.append(v)
236            if varied > 0 and varied != len(varlist)+1:
237                msg += "\nNot all variables refined in equivalence:\n\t"
238                s = ""
239                for v in varlist:
240                    if s != "": s+= " & "
241                    s += str(v)           
242                msg += str(mv) + " => " + s
243                msg += '\nNot refined: ' + notvaried + '\n'
244            if zeromult:
245                msg += "\nZero multiplier is invalid in equivalence:\n\t"
246                s = ""
247                for v in varlist:
248                    if s != "": s+= " & "
249                    s += str(v)           
250                msg += str(mv) + " => " + s + '\n'
251
252    #print 'indepVarList',indepVarList
253    #print 'depVarList',depVarList
254    # check for errors:
255    inboth = set(indepVarList).intersection(set(depVarList))
256    if len(inboth) > 0:
257        msg += "\nThe following parameters(s) are used as both dependent and independent variables in Equivalence relations:\n"
258        s = ''
259        for var in inboth:
260            if s != "": s+= ", "
261            s += str(v)
262        msg += '\t'+ s + '\n'
263    if len(multdepVarList) > 0:
264        msg += "\nThe following parameters(s) are used in multiple Equivalence relations as dependent variables:\n"
265        s = ''
266        for var in multdepVarList:
267            if s != "": s+= ", "
268            s += str(v)           
269        msg += '\t'+ s + '\n'
270    equivVarList = list(set(indepVarList).union(set(depVarList)))
271    #print 'equivVarList',equivVarList
272    inboth = set(fixedVarList).intersection(set(equivVarList))
273    if len(inboth) > 0:
274        msg += "\nError! The following variables are used in both Equivalence and Fixed constraints:\n"
275        s = ''
276        for var in inboth:
277            if s != "": s+= ", "
278            s += str(v)
279        msg += '\t'+ s + '\n'
280
281    # scan through parameters in each relationship. Are all varied? If only some are
282    # varied, create an error message.
283    # If all are varied and this is a constraint equation, then set VaryFree flag
284    # so that newly created relationships will be varied
285    for group,varlist in zip(groups,parmlist):
286        if len(varlist) == 1: continue
287        VaryFree = False
288        for rel in group:
289            varied = 0
290            notvaried = ''
291            for var in constrDict[rel]:
292                if var in varyList:
293                    varied += 1
294                else:
295                    if notvaried: notvaried += ', '
296                    notvaried += var
297                if var in fixedVarList:
298                    msg += '\nError: parameter '+var+" is Fixed and used in a constraint:\n\t"
299                    msg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
300                if var in equivVarList:
301                    msg += '\nError: parameter '+var+" is Equivalenced and used in a constraint:\n\t"
302                    msg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
303            if varied > 0 and varied != len(constrDict[rel]):
304                msg += "\nNot all variables refined in constraint:\n\t"
305                msg += FormatConstraint(constrDict[rel],fixedList[rel])
306                msg += '\nNot refined: ' + notvaried + '\n'
307            if fixedList[rel] is not None and varied > 0:
308                VaryFree = True
309    # if there were errors found, go no farther
310    if msg:
311        print ' *** ERROR in constraint definitions! ***'
312        print msg
313        raise Exception
314               
315    # now process each group and create the relations that are needed to form
316    # non-singular square matrix
317    for group,varlist in zip(groups,parmlist):
318        if len(varlist) == 1: continue
319        arr = MakeArray(constrDict, group, varlist)
320        constrArr = FillInMissingRelations(arr,len(group))
321        mapvar = []
322        group = group[:]
323        # scan through all generated and input variables, add to the varied list
324        # all the new parameters where VaryFree has been set or where all the
325        # dependent parameters are varied. Check again for inconsistent variable use
326        # for new variables -- where varied and unvaried parameters get grouped
327        # together. I don't think this can happen when not flagged before, but
328        # it does not hurt to check again.
329        for i in range(len(varlist)):
330            varied = 0
331            notvaried = ''
332            if len(group) > 0:
333                rel = group.pop(0)
334                fixedval = fixedList[rel]
335                for var in constrDict[rel]:
336                    if var in varyList:
337                        varied += 1
338                    else:
339                        if notvaried: notvaried += ', '
340                        notvaried += var
341            else:
342                fixedval = None
343            if fixedval is None:
344                varname = paramPrefix + str(consNum)
345                mapvar.append(varname)
346                consNum += 1
347                if VaryFree or varied > 0:
348                    varyList.append(varname)
349            else:
350                mapvar.append(fixedval)
351            if varied > 0 and notvaried != '':
352                msg += "\nNot all variables refined in generated constraint\n\t"
353                msg += '\nNot refined: ' + notvaried + '\n'
354        dependentParmList.append(varlist)
355        arrayList.append(constrArr)
356        invarrayList.append(np.linalg.inv(constrArr))
357        indParmList.append(mapvar)
358    if msg:
359        print ' *** ERROR in constraint definitions! ***'
360        print msg
361        print VarRemapShow(varyList)
362        raise Exception
363    # setup dictionary containing the fixed values
364    global fixedDict
365    # key is original ascii string, value is float
366    for fixedval in fixedList:
367        if fixedval:
368            fixedDict[fixedval] = float(fixedval)
369
370    if debug: # on debug, show what is parsed & generated, semi-readable
371        print 50*'-'
372        for group,varlist,multarr,inv,mapvar in zip(groups,parmlist,arrayList,invarrayList,indParmList):
373            print '\n*** relation(s) in group:',group,'\tvars in group:',varlist
374            print 'new parameters:', mapvar
375            print 'Input relationship matrix'
376            print multarr[:len(group)]
377            added = len(group) - len(varlist)
378            if added < 0:
379                print 'added relationship rows'
380                print multarr[added:]
381            print 'Inverse relationship matrix'
382            print inv
383
384def StoreEquivalence(independentVar,dependentList):
385    '''Takes a list of dependent parameter(s) and stores their
386    relationship to a single independent parameter (independentVar)
387    input:
388       independentVar -- name of parameter that will set others (may
389         be varied)
390       dependentList -- list of parameter that will set from
391         independentVar (may not be varied). Each item can be a parameter
392         name or a tuple containing a name and multiplier:
393         ('parm1',('parm2',.5),)
394    '''
395   
396    global dependentParmList,arrayList,invarrayList,indParmList
397    mapList = []
398    multlist = []
399    for var in dependentList:
400        if isinstance(var, basestring):
401            mult = 1.0
402        elif len(var) == 2:
403            var,mult = var
404        else:
405            raise Exception, "Cannot parse "+repr(var) + " as var or (var,multiplier)"
406        mapList.append(var)
407        multlist.append(tuple((mult,)))
408    # added relationships to stored values
409    arrayList.append(None)
410    invarrayList.append(np.array(multlist))
411    indParmList.append(tuple((independentVar,)))
412    dependentParmList.append(mapList)
413    return
414
415# def SetVaryFlags(varyList,fixedList):
416#     '''Adds independent variables to the varyList provided that all
417#     dependent variables are being varied.
418#     Ignores independent variables where no dependent variables are
419#     being varied.
420#     Returns a non-empty text message where some but not all dependent
421#     variables are being varied.
422#     '''
423#     global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
424#     msg = ""
425#     for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
426#         for mv in mapvars:
427#             varied = []
428#             notvaried = []
429#             if mv in fixedDict: continue
430#             if multarr is None:
431#                 if mv in varyList:
432#                     varied.append(mv)
433#                 else:
434#                     notvaried.append(mv)
435#             for v in varlist:
436#                 if v in varyList:
437#                     varied.append(v)
438#                 else:
439#                     notvaried.append(v)
440#             if len(varied) > 0 and len(notvaried) > 0:
441#                 if msg != "": msg += "\n"
442#                 msg += "Error: inconsistent use of parameter " + mv
443#                 msg += "\n  varied: "
444#                 for v in varied: msg += v
445#                 msg += "\n  not varied: "
446#                 for v in notvaried: msg += v   
447#             #elif len(varied) > 0 and multarr is not None:
448#             #    varyList.append(mv)
449#     return msg
450
451def GetDependentVars():
452    '''Return a list of dependent variables: e.g. variables that are
453    constrained in terms of other variables'''
454    dependentVars = []
455    global dependentParmList
456    for lst in dependentParmList:
457        for itm in lst: dependentVars.append(itm)
458    return dependentVars
459
460def GetIndependentVars():
461    '''Return a list of independent variables: e.g. variables that are
462    created by constrains of other variables'''
463    independentVars = []
464    global indParmList,fixedDict
465    for lst in indParmList:
466        for name in lst:
467            if name in fixedDict: continue
468            independentVars.append(name)
469    return independentVars
470
471def PrintIndependentVars(parmDict,varyList,sigDict,PrintAll=False):
472    '''Print the values and uncertainties on the independent variables'''
473    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
474    printlist = []
475    mapvars = GetIndependentVars()
476    for i,name in enumerate(mapvars):
477        if name in fixedDict: continue
478        if PrintAll or name in varyList:
479            sig = sigDict.get(name)
480            printlist.append([name,parmDict[name],sig])
481    if len(printlist) == 0: return
482    s1 = ''
483    print 130*'-'
484    print "Variables generated by constraints"
485    printlist.append(3*[None])
486    for name,val,esd in printlist:
487        if len(s1) > 110 or name is None:
488            print
489            print s1
490            print s2
491            print s3
492            s1 = ''
493            if name is None: break
494        if s1 == "":
495            s1 = ' name  :'
496            s2 = ' value :'
497            s3 = ' sig   :'
498        s1 += '%12s' % (name)
499        s2 += '%12.6f' % (val)
500        if esd is None:
501            s3 += '%12s' % ('n/a')
502        else:   
503            s3 += '%12.6f' % (esd)
504
505def ComputeDepESD(covMatrix,varyList,parmDict):
506    '''Compute uncertainties for dependent parameters from independent ones
507    returns a dictionary containing the esd values for dependent parameters
508    '''
509    # I think this fails for equivalencies: the ESD for the master (independent variable)
510    # needs to be adjusted.
511    sigmaDict = {}
512    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
513        #if invmultarr is None: continue # probably not needed
514        valuelist = [parmDict[var] for var in mapvars]
515        # get the v-covar matrix for independent parameters
516        vcov = np.zeros((len(mapvars),len(mapvars)))
517        for i1,name1 in enumerate(mapvars):
518            if name1 not in varyList: continue
519            iv1 = varyList.index(name1)
520            for i2,name2 in enumerate(mapvars):
521                if name2 not in varyList: continue
522                iv2 = varyList.index(name2)
523                vcov[i1][i2] = covMatrix[iv1][iv2]
524        # vec is the vector that multiplies each of the independent values
525        for v,vec in zip(varlist,invmultarr):
526            sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec)))
527    return sigmaDict
528
529def FormatConstraint(RelDict,RelVal):
530    '''Formats a Constraint or Function for use in a convenient way'''
531    s = [""]
532    for var,val in RelDict.items():
533        if len(s[-1]) > 60: s.append(' ')
534        m = val
535        if s[-1] != "" and m >= 0:
536            s[-1] += ' + '
537        elif s[-1] != "":
538            s[-1] += ' - '
539            m = abs(m)
540        s[-1] += '%.3f*%s '%(m,var)
541    if len(s[-1]) > 60: s.append(' ')
542    if RelVal is None:
543        s[-1] += ' = New variable'
544    else:
545        s[-1] += ' = ' + RelVal
546    s1 = ''
547    for s2 in s:
548        if s1 != '': s1 += '\n\t'
549        s1 += s2
550    return s1
551
552def VarRemapShow(varyList):
553    '''List out the saved relationships.
554    Returns a string containing the details.
555    '''
556    s = ''
557    if len(fixedVarList) > 0:
558        s += 'Fixed Variables:\n'
559        for v in fixedVarList:
560            s += '    ' + v + '\n'
561    s += 'Variable mapping relations:\n'
562    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
563    for varlist,mapvars,multarr,invmultarr in zip(
564        dependentParmList,indParmList,arrayList,invarrayList):
565        for i,mv in enumerate(mapvars):
566            if multarr is None:
567                s += '  ' + str(mv) + ' is equivalent to parameter(s): '
568                j = 0
569                for v,m in zip(varlist,invmultarr):
570                    #print v,m[0]
571                    if j > 0: s += '  & '
572                    j += 1
573                    s += str(v)
574                    if m != 1:
575                        s += " / " + str(m[0])                       
576                s += '\n'
577                continue
578            s += %s = ' % mv
579            j = 0
580            for m,v in zip(multarr[i,:],varlist):
581                if m == 0: continue
582                if j > 0: s += ' + '
583                j += 1
584                s += '(%s * %s)' % (m,v)
585            if mv in varyList: s += ' VARY'
586            s += '\n'
587    s += 'Inverse variable mapping relations:\n'
588    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
589        for i,mv in enumerate(varlist):
590            s += %s = ' % mv
591            j = 0
592            for m,v in zip(invmultarr[i,:],mapvars):
593                if m == 0: continue
594                if j > 0: s += ' + '
595                j += 1
596                s += '(%s * %s)' % (m,v)
597            s += '\n'
598    return s
599
600def Dict2Deriv(varyList,derivDict,dMdv):
601    '''Compute derivatives for Independent Parameters from the
602    derivatives for the original parameters
603    '''
604    global dependentParmList,arrayList,invarrayList,indParmList,invarrayList
605    for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
606        for i,name in enumerate(mapvars):
607            # grouped variables: need to add in the derv. w/r
608            # dependent variables to the independent ones
609            if name not in varyList: continue # skip if independent var not varied
610            if multarr is None:
611                for v,m in zip(varlist,invmultarr):
612                    #print 'add derv',v,'/',m[0],'to derv',name
613                    if m == 0: continue
614                    dMdv[varyList.index(name)] += derivDict[v] / m[0]
615            else:
616                for v,m in zip(varlist,multarr[i,:]):
617                    #print 'add derv',v,'*',m,'to derv',name
618                    if m == 0: continue
619                    dMdv[varyList.index(name)] += m * derivDict[v]
620
621def Map2Dict(parmDict,varyList):
622    '''Create (or update) the Independent Parameters from the original
623    set of Parameters
624
625    Removes dependent variables from the varyList
626
627    This should be done once, before any variable refinement is done
628    to complete the parameter dictionary with the Independent Parameters
629    '''
630    # process the Independent vars: remove dependent ones from varylist
631    # and then compute values for the Independent ones from their dependents
632    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
633    for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
634        for item in varlist:
635            try:
636                varyList.remove(item)
637            except ValueError:
638                pass
639        if multarr is None: continue
640        valuelist = [parmDict[var] for var in varlist]
641        parmDict.update(zip(mapvars,
642                            np.dot(multarr,np.array(valuelist)))
643                        )
644    # now remove fixed variables from the varyList
645    global fixedVarList
646    for item in fixedVarList:
647        try:
648            varyList.remove(item)
649        except ValueError:
650            pass
651    # Set constrained parameters to their fixed values
652    parmDict.update(fixedDict)
653
654def Dict2Map(parmDict,varyList):
655    '''Convert the remapped values back to the original parameters
656   
657    This should be done to apply constraints to parameter values (use
658    Map2Dict once first). It should be done as part of the Model function
659    evaluation, before any computation is done
660    '''
661    #I think this needs fixing to update parmDict with new values
662    #   from the constraints based on what is in varyList - RVD. Don't think so -- BHT
663    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
664    # reset fixed values (should not be needed, but very quick)
665    # - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended
666    # not needed, but also not a problem - BHT
667    parmDict.update(fixedDict)
668    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
669        #if invmultarr is None: continue
670        valuelist = [parmDict[var] for var in mapvars]
671        parmDict.update(zip(varlist,
672                            np.dot(invmultarr,np.array(valuelist)))
673                        )
674
675#======================================================================
676# internal routines follow (these routines are unlikely to be called
677# from outside the module)
678def FillInMissingRelations(arrayin,nvars):
679    '''Fill in unspecified rows in array with non-colinear unit vectors
680    arrayin is a square array, where only the first nvars rows are defined.
681   
682    Returns a new array where all rows are defined
683
684    Throws an exception if there is no non-signular result
685    (likely because two or more input rows are co-linear)
686    '''
687    a = arrayin.copy()
688    n = nvars
689    # fill in the matrix with basis vectors not colinear with any of the starting set
690    for i in range(n,len(a)):
691        try:
692            a[n] = VectorProduct(a[:n])
693        except:
694            raise Exception,"VectorProduct failed, singular input?"
695        n += 1
696    # use the G-S algorithm to compute a complete set of unit vectors orthogonal
697    # to the first in array
698    gs = GramSchmidtOrtho(a) 
699    n = nvars
700    # now replace the created vectors with the ones least correlated to the
701    # initial set
702    vlist = [v for v in gs[1:]] # drop the first row
703    for j in range(n,len(a)):
704        minavg = None
705        droplist = []
706        for k in range(len(vlist)):
707            v = vlist[k]
708            avgcor = 0
709            for i in range(n):
710                cor = np.dot(a[i],vlist[k])/(np.linalg.norm(a[i]) * np.linalg.norm(vlist[k]) )
711                if np.allclose(cor, 1.0):
712                    droplist.append(k) # found a vector co-linear to one we already have
713                                       # don't need it again
714                    #vlist.pop(k)
715                    break 
716                avgcor += cor
717            else:
718                if minavg == None:
719                    minavg = abs(avgcor/n)
720                    minnum = k
721                elif abs(avgcor/n) < minavg:
722                    minavg = abs(avgcor/n)
723                    minnum = k
724        if minavg == None:
725            raise Exception,"Failed to find a non-colinear G-S vector for row %d. Should not happen!" % n
726        a[j] = vlist[minnum]
727        droplist.append(minnum) # don't need to use this vector again
728        for i in sorted(droplist,reverse=True):
729            vlist.pop(i) 
730        n += 1
731    return a
732
733
734def MakeArray(constrDict, group, varlist):
735    """Convert the multipliers in a constraint group to an array of
736    coefficients and place in a square numpy array, adding empty rows as
737    needed.
738
739    Throws an exception if some sanity checks fail.
740    """
741    # do some error checks
742    if len(varlist) < len(group): # too many relationships -- no can do
743        raise Exception,'The number of relationships (%d) exceeds the number of parameters (%d):\n\t%s\n\t%s'% (
744            len(group),len(varlist),group,varlist)
745    # put the coefficients into an array
746    multarr = np.zeros(2*[len(varlist),])
747    i = 0
748    for cnum in group:
749        cdict = constrDict[cnum]
750        j = 0
751        for var in varlist:
752            m = cdict.get(var)
753            if m != None:
754                multarr[i,j] = m
755            j += 1
756        i += 1
757    return multarr
758
759def GramSchmidtOrtho(arrayin):
760    '''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to
761    find orthonormal unit vectors relative to first row.
762    input:
763       arrayin: a 2-D non-signular square array
764    returns:
765       a orthonormal set of unit vectors as a square array
766    '''
767    def proj(a,b):
768        'Projection operator'
769        return a*(np.dot(a,b)/np.dot(a,a))
770    a = arrayin.copy()
771    for j in range (len(a)):
772        for i in range(j):
773            a[j] = a[j] - proj(a[i],a[j])
774        if np.allclose(np.linalg.norm(a[j]),0.0):
775            raise Exception,"Singular input to GramSchmidtOrtho"
776        a[j] = a[j]/np.linalg.norm(a[j])
777    return a
778
779def VectorProduct(vl):
780    '''Evaluate the n-dimensional "cross" product of the list of vectors in vl
781    vl can be any length. The length of each vector is arbitrary, but
782    all must be the same length. See http://en.wikipedia.org/wiki/Levi-Civita_symbol
783
784    This appears to return a vector perpendicular to the supplied set.
785
786    Throws an exception if a vector can not be obtained because the input set of
787    vectors is already complete or contains any redundancies.
788   
789    Uses LeviCitvitaMult recursively to obtain all permutations of vector elements
790    '''
791    i = 0
792    newvec = []
793    for e in vl[0]:
794        i += 1
795        tl = [([i,],1),]
796        seq = LeviCitvitaMult(tl ,vl)
797        tsum = 0
798        for terms,prod in seq:
799            tsum += EvalLCterm(terms) * prod
800        newvec.append(tsum)
801    if np.allclose(newvec,0.0):
802        raise Exception,"VectorProduct failed, singular or already complete input"
803    return newvec
804
805def LeviCitvitaMult(tin,vl):
806    '''Recursion formula to compute all permutations of vector element products
807    The first term in each list is the indices of the Levi-Civita symbol, note
808    that if an index is repeated, the value is zero, so the element is not included.
809    The second term is the product of the vector elements
810    '''
811    v = vl[0]
812    vl1 = vl[1:]       
813
814    j = 0
815    tl = []
816    for e in vl[0]:
817        j += 1
818        for ind,vals in tin:
819            if j in ind: continue
820            tl.append((ind+[j,],vals*e))
821    if len(vl1):
822        return LeviCitvitaMult(tl,vl1)
823    else:
824        return tl
825
826def EvalLCterm(term):
827    '''Evaluate the Levi-Civita symbol Epsilon(i,j,k,...)'''
828    p = 1
829    for i in range(len(term)):
830        for j in range(i+1,len(term)):
831            p *= (term[j] - term[i])
832            if p == 0: return 0
833    return p/abs(p)
834
835
836if __name__ == "__main__":
837    import sys
838    parmdict = {}
839    constrDict = [
840        {'0:12:Scale': 2.0, '0:11:Scale': 1.0, '0:14:Scale': 4.0, '0:13:Scale': 3.0, '0:0:Scale': 0.5},
841        {'0:0:eA': 0.0},
842        {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
843        {'1::C(10,0,1)': 1.0, '2::C(10,0,1)': 1.0},
844        {'1::AUiso:0': 1.0, '0::AUiso:0': 1.0},
845        {'0::A0': 0.0}
846        ]
847    fixedList = ['5.0', '0', None, None, '1.0', '0']
848    StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), ))
849    #StoreEquivalence('1::atomx:3',('2::atomx:3', ('2::atomz:3',2,), )) # error: dependent & independent vars mixed
850    #StoreEquivalence('1::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) # error: dependent vars repeated
851    #StoreEquivalence('0:1:eA',('0:0:eA',)) # error: equiv & fixed
852    #StoreEquivalence('0:99:Scale',('0:12:Scale',)) # error: equiv & constrained
853    #StoreEquivalence('0:12:Scale',('0:99:Scale',)) # error: equiv & constrained
854    varylist = ['2::atomx:3',
855                '2::C(10,6,1)', '1::C(10,6,1)',
856                '2::atomy:3', '2::atomz:3',
857                '0:12:Scale', '0:11:Scale', '0:14:Scale', '0:13:Scale', '0:0:Scale']
858    groups,parmlist = GroupConstraints(constrDict)
859    GenerateConstraints(groups,parmlist,varylist,constrDict,fixedList)
860    print VarRemapShow(varylist)
861    parmdict.update( {
862        '0:12:Scale': 1.0, '0:11:Scale': 1.0, '0:14:Scale': 1.0, '0:13:Scale': 1.0, '0:0:Scale': 2.0,
863        '0:0:eA': 0.0,
864        '2::C(10,6,1)': 0.2, '1::C(10,6,1)': 0.3,
865        '1::C(10,0,1)': 0.2, '2::C(10,0,1)': 0.3,
866        '1::AUiso:0': 0.02, '0::AUiso:0': 0.03,
867        '0::A0': 0.0,
868        '2::atomx:3':0.23,'2::atomy:3':-.23, '2::atomz:3':-0.11,
869        })
870    print 'parmdict start',parmdict
871    print 'varylist start',varylist
872    before = parmdict.copy()
873    Map2Dict(parmdict,varylist)
874    print 'parmdict before and after Map2Dict'
875    print '  key / before / after'
876    for key in sorted(parmdict.keys()):
877        print '  '+key,'\t',before.get(key),'\t',parmdict[key]
878    print 'varylist after',varylist
879    before = parmdict.copy()
880    Dict2Map(parmdict,varylist)
881    print 'after Dict2Map'
882    print '  key / before / after'
883    for key in sorted(parmdict.keys()):
884        print '  '+key,'\t',before.get(key),'\t',parmdict[key]
885    dMdv = len(varylist)*[0]
886    deriv = {}
887    for i,v in enumerate(parmdict.keys()): deriv[v]=i
888    Dict2Deriv(varylist,deriv,dMdv)
Note: See TracBrowser for help on using the repository browser.