source: trunk/GSASIImapvars.py @ 570

Last change on this file since 570 was 570, checked in by toby, 9 years ago

Add checks on constraints: if problems display error window when displayed, added/edited & on Refine; handle GPX w/o doPawley flag
File size: 46.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 GroupConstraints and 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 InitVars, 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 is not allowed:
91        StoreEquivalence('x',('y',))
92        StoreEquivalence('y',('z',))
93   Use StoreEquivalence before calling GenerateConstraints or
94      CheckConstraints
95
96   To check that input in internally consistent, use CheckConstraints
97
98   To show a summary of the parameter remapping, one calls
99      VarRemapShow
100
101   To determine values for the parameters created in this module, one
102      calls Map2Dict. This will not apply contraints.
103
104   To take values from the new independent parameters and constraints,
105      one calls Dict2Map. This will apply contraints.
106
107   Use Dict2Deriv to determine derivatives on independent parameters
108      from those on dependent ones
109     
110   Use ComputeDepESD to compute uncertainties on dependent variables
111
112Global Variables:
113   dependentParmList: contains a list by group of lists of
114     parameters used in the group. Note that parameters listed in
115     dependentParmList should not be refined as they will not affect
116     the model
117   indParmList: a list of groups of Independent parameters defined in
118     each group. This contains both parameters used in parameter
119     redefinitions as well as names of generated new parameters.
120
121   fixedVarList: a list of variables that have been 'fixed'
122     by defining them as equal to a constant (::var: = 0). Note that
123     the constant value is ignored at present. These variables are
124     later removed from varyList which prevents them from being refined.
125     Unlikely to be used externally.
126   arrayList: a list by group of relationship matrices to relate
127     parameters in dependentParmList to those in indParmList. Unlikely
128     to be used externally.
129   invarrayList: a list by group of relationship matrices to relate
130     parameters in indParmList to those in dependentParmList. Unlikely
131     to be used externally.
132   fixedDict: a dictionary containing the fixed values corresponding
133     to parameter equations.  The dict key is an ascii string, but the
134     dict value is a float.  Unlikely to be used externally.
135"""
136
137import numpy as np
138# data used for constraints;
139debug = False # turns on printing as constraint input is processed
140# note that constraints are stored listed by contraint groups, where each constraint
141# group contains those parameters that must be handled together
142dependentParmList = [] # contains a list of parameters in each group
143# note that parameters listed in dependentParmList should not be refined
144arrayList = [] # a list of of relationship matrices
145invarrayList = [] # a list of inverse relationship matrices
146indParmList = [] # a list of names for the new parameters
147fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations
148               # key is original ascii string, value is float
149fixedVarList = [] # List of variables that should not be refined
150
151# prefix for parameter names
152paramPrefix = "::constr:"
153consNum = 0 # number of the next constraint to be created
154
155def InitVars():
156    '''Initializes all constraint information'''
157    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,consNum
158    dependentParmList = [] # contains a list of parameters in each group
159    arrayList = [] # a list of of relationship matrices
160    invarrayList = [] # a list of inverse relationship matrices
161    indParmList = [] # a list of names for the new parameters
162    fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations
163    consNum = 0 # number of the next constraint to be created
164    fixedVarList = []
165
166def GroupConstraints(constrDict):
167    """divide the constraints into groups that share no parameters.
168    Input
169       constrDict: a list of dicts defining relationships/constraints
170       constrDict = [{<constr1>}, {<constr2>}, ...]
171       where {<constr1>} is {'param1': mult1, 'param2': mult2,...}
172    Returns two lists of lists:
173      a list of relationship list indices for each group pointing to
174        elements in constrDict and
175      a list containing the parameters used in each group
176      """
177    assignedlist = [] # relationships that have been used
178    groups = [] # contains a list of grouplists
179    ParmList = []
180    for i,consi in enumerate(constrDict):
181        if i in assignedlist: continue # already in a group, skip
182        # starting a new group
183        grouplist = [i,]
184        assignedlist.append(i)
185        groupset = set(consi.keys())
186        changes = True # always loop at least once
187        while(changes): # loop until we can't find anything to add to the current group
188            changes = False # but don't loop again unless we find something
189            for j,consj in enumerate(constrDict):
190                if j in assignedlist: continue # already in a group, skip
191                if len(set(consj.keys()) & groupset) > 0: # true if this needs to be added
192                    changes = True
193                    grouplist.append(j)
194                    assignedlist.append(j)
195                    groupset = groupset | set(consj.keys())
196        group = sorted(grouplist)
197        varlist = sorted(list(groupset))
198        groups.append(group)
199        ParmList.append(varlist)
200    return groups,ParmList
201
202def CheckConstraints(varyList,constrDict,fixedList):
203    '''Takes a list of relationship entries comprising a group of
204    constraints and checks for inconsistencies such as conflicts in
205    parameter/variable definitions and or inconsistently varied parameters.
206    Input: see GenerateConstraints
207    Output: returns two strings
208      the first lists conflicts internal to the specified constraints
209      the second lists conflicts where the varyList specifies some
210        parameters in a constraint, but not all
211      If there are no errors, both strings will be empty
212    '''
213    global dependentParmList,arrayList,invarrayList,indParmList,consNum
214    errmsg = ''
215    warnmsg = ''
216    fixVlist = []
217    # process fixed (held) variables
218    for cdict in constrDict:
219        if len(cdict) == 1:
220            fixVlist.append(cdict.keys()[0])
221   
222    # process equivalences: make a list of dependent and independent vars
223    #    and check for repeated uses (repetition of a parameter as an
224    #    independent var is OK)
225    indepVarList = []
226    depVarList = []
227    multdepVarList = []
228    for varlist,mapvars,multarr,invmultarr in zip(
229        dependentParmList,indParmList,arrayList,invarrayList):
230        if multarr is None: # an equivalence
231            zeromult = False
232            for mv in mapvars:
233                varied = 0
234                notvaried = ''
235                if mv in varyList:
236                    varied += 1
237                else:
238                    if notvaried: notvaried += ', '
239                    notvaried += mv
240                if mv not in indepVarList: indepVarList.append(mv)
241                for v,m in zip(varlist,invmultarr):
242                    if m == 0: zeromult = True
243                    if v in varyList:
244                        varied += 1
245                    else:
246                        if notvaried: notvaried += ', '
247                        notvaried += v
248                    if v in depVarList:
249                        multdepVarList.append(v)
250                    else:
251                        depVarList.append(v)
252            if varied > 0 and varied != len(varlist)+1:
253                warnmsg += "\nNot all variables refined in equivalence:\n\t"
254                s = ""
255                for v in varlist:
256                    if s != "": s+= " & "
257                    s += str(v)           
258                warnmsg += str(mv) + " => " + s
259                warnmsg += '\nNot refined: ' + notvaried + '\n'
260            if zeromult:
261                errmsg += "\nZero multiplier is invalid in equivalence:\n\t"
262                s = ""
263                for v in varlist:
264                    if s != "": s+= " & "
265                    s += str(v)           
266                errmsg += str(mv) + " => " + s + '\n'
267
268    #print 'indepVarList',indepVarList
269    #print 'depVarList',depVarList
270    # check for errors:
271    inboth = set(indepVarList).intersection(set(depVarList))
272    if len(inboth) > 0:
273        errmsg += "\nThe following parameters(s) are used as both dependent and independent variables in Equivalence relations:\n"
274        s = ''
275        for var in inboth:
276            if s != "": s+= ", "
277            s += str(v)
278        errmsg += '\t'+ s + '\n'
279    if len(multdepVarList) > 0:
280        errmsg += "\nThe following parameters(s) are used in multiple Equivalence relations as dependent variables:\n"
281        s = ''
282        for var in multdepVarList:
283            if s != "": s+= ", "
284            s += str(v)           
285        errmsg += '\t'+ s + '\n'
286    equivVarList = list(set(indepVarList).union(set(depVarList)))
287    #print 'equivVarList',equivVarList
288    inboth = set(fixVlist).intersection(set(equivVarList))
289    if len(inboth) > 0:
290        errmsg += "\nThe following parameter(s) are used in both Equivalence and Fixed constraints:\n"
291        s = ''
292        for var in inboth:
293            if s != "": s+= ", "
294            s += str(v)
295        errmsg += '\t'+ s + '\n'
296
297    groups,parmlist = GroupConstraints(constrDict)
298    # scan through parameters in each relationship. Are all varied? If only some are
299    # varied, create a warning message.
300    for group,varlist in zip(groups,parmlist):
301        if len(varlist) == 1: continue
302        VaryFree = False
303        for rel in group:
304            varied = 0
305            notvaried = ''
306            for var in constrDict[rel]:
307                if var in varyList:
308                    varied += 1
309                else:
310                    if notvaried: notvaried += ', '
311                    notvaried += var
312                if var in fixVlist:
313                    errmsg += '\nParameter '+var+" is Fixed and used in a constraint:\n\t"
314                    errmsg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
315                if var in equivVarList:
316                    errmsg += '\nParameter '+var+" is Equivalenced and used in a constraint:\n\t"
317                    errmsg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
318            if varied > 0 and varied != len(constrDict[rel]):
319                warnmsg += "\nNot all variables refined in constraint:\n\t"
320                warnmsg += FormatConstraint(constrDict[rel],fixedList[rel])
321                warnmsg += '\nNot refined: ' + notvaried + '\n'
322    if errmsg or warnmsg: return errmsg,warnmsg
323
324    # now look for process each group and create the relations that are needed to form
325    # non-singular square matrix
326    for group,varlist in zip(groups,parmlist):
327        if len(varlist) == 1: continue
328        try: 
329            arr = MakeArray(constrDict, group, varlist)
330        except:
331            errmsg += "\nOver-constrained input. "
332            errmsg += "There are more constraints" + str(len(group))
333            errmsg += "\n\tthan variables" + str(len(varlist)) + "\n"
334            for rel in group:
335                errmsg += FormatConstraint(constrDict[rel],fixedList[rel])
336                errmsg += "\n"
337        try:
338            constrArr = FillInMissingRelations(arr,len(group))
339        except Exception,msg:
340            if msg.message.startswith('VectorProduct'):
341                errmsg += "\nSingular input. "
342                errmsg += "This set of constraints is not linearly independent:\n\n"
343            else:
344                errmsg += "\nInconsistent input. "
345                errmsg += "Cound not generate a set of non-singular constraints"
346                errmsg += "\n\tfor this set of input:\n\n"
347            for rel in group:
348                errmsg += FormatConstraint(constrDict[rel],fixedList[rel])
349                errmsg += "\n"
350            #import traceback
351            #print traceback.format_exc()
352            continue
353       
354        mapvar = []
355        group = group[:]
356        # scan through all generated and input variables
357        # Check again for inconsistent variable use
358        # for new variables -- where varied and unvaried parameters get grouped
359        # together. I don't think this can happen when not flagged before, but
360        # it does not hurt to check again.
361        for i in range(len(varlist)):
362            varied = 0
363            notvaried = ''
364            if len(group) > 0:
365                rel = group.pop(0)
366                fixedval = fixedList[rel]
367                for var in constrDict[rel]:
368                    if var in varyList:
369                        varied += 1
370                    else:
371                        if notvaried: notvaried += ', '
372                        notvaried += var
373            else:
374                fixedval = None
375            if fixedval is None:
376                varname = paramPrefix + str(consNum)
377                mapvar.append(varname)
378                consNum += 1
379                if VaryFree or varied > 0:
380                    varyList.append(varname)
381            else:
382                mapvar.append(fixedval)
383            if varied > 0 and notvaried != '':
384                warnmsg += "\nNot all variables refined in generated constraint"
385                warnmsg += '\nPlease report this unexpected error\n'
386                for rel in group:
387                    warnmsg += FormatConstraint(constrDict[rel],fixedList[rel])
388                    warnmsg += "\n"
389                warnmsg += '\n\tNot refined: ' + notvaried + '\n'
390        try:
391            np.linalg.inv(constrArr)
392        except:
393            errmsg += "\nSingular input. "
394            errmsg += "The following constraints are not "
395            errmsg += "linearly independent\n\tor do not "
396            errmsg += "allow for generation of a non-singular set\n"
397            errmsg += 'Please report this unexpected error\n'
398            for rel in group:
399                errmsg += FormatConstraint(constrDict[rel],fixedList[rel])
400                errmsg += "\n"
401    return errmsg,warnmsg
402
403def GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList):
404    '''Takes a list of relationship entries comprising a group of
405    constraints and builds the relationship lists and their inverse
406    and stores them in global variables Also checks for internal
407    conflicts or inconsistencies in parameter/variable definitions.
408    Input:
409       groups,parmlist: see GroupConstraints
410       varyList: a list of parameters that will be varied
411       constrDict: a list of dicts defining relationships/constraints
412         constrDict = [{<constr1>}, {<constr2>}, ...]
413         where {<constr1>} is {'param1': mult1, 'param2': mult2,...}
414       fixedList: a list of values for each constraint/variable in
415          constrDict, value is either a float (contraint) or None (for
416          a new variable).
417    '''
418    global dependentParmList,arrayList,invarrayList,indParmList,consNum
419    msg = ''
420
421    # process fixed (held) variables
422    for cdict in constrDict:
423        if len(cdict) == 1:
424            fixedVarList.append(cdict.keys()[0])
425    #print 'fixedVarList',fixedVarList
426   
427    # process equivalences: make a list of dependent and independent vars
428    #    and check for repeated uses (repetition of a parameter as an
429    #    independent var is OK)
430    indepVarList = []
431    depVarList = []
432    multdepVarList = []
433    for varlist,mapvars,multarr,invmultarr in zip(
434        dependentParmList,indParmList,arrayList,invarrayList):
435        if multarr is None:
436            zeromult = False
437            for mv in mapvars:
438                varied = 0
439                notvaried = ''
440                if mv in varyList:
441                    varied += 1
442                else:
443                    if notvaried: notvaried += ', '
444                    notvaried += mv
445                if mv not in indepVarList: indepVarList.append(mv)
446                for v,m in zip(varlist,invmultarr):
447                    if m == 0: zeromult = True
448                    if v in varyList:
449                        varied += 1
450                    else:
451                        if notvaried: notvaried += ', '
452                        notvaried += v
453                    if v in depVarList:
454                        multdepVarList.append(v)
455                    else:
456                        depVarList.append(v)
457            if varied > 0 and varied != len(varlist)+1:
458                msg += "\nNot all variables refined in equivalence:\n\t"
459                s = ""
460                for v in varlist:
461                    if s != "": s+= " & "
462                    s += str(v)           
463                msg += str(mv) + " => " + s
464                msg += '\nNot refined: ' + notvaried + '\n'
465            if zeromult:
466                msg += "\nZero multiplier is invalid in equivalence:\n\t"
467                s = ""
468                for v in varlist:
469                    if s != "": s+= " & "
470                    s += str(v)           
471                msg += str(mv) + " => " + s + '\n'
472
473    #print 'indepVarList',indepVarList
474    #print 'depVarList',depVarList
475    # check for errors:
476    inboth = set(indepVarList).intersection(set(depVarList))
477    if len(inboth) > 0:
478        msg += "\nThe following parameters(s) are used as both dependent and independent variables in Equivalence relations:\n"
479        s = ''
480        for var in inboth:
481            if s != "": s+= ", "
482            s += str(v)
483        msg += '\t'+ s + '\n'
484    if len(multdepVarList) > 0:
485        msg += "\nThe following parameters(s) are used in multiple Equivalence relations as dependent variables:\n"
486        s = ''
487        for var in multdepVarList:
488            if s != "": s+= ", "
489            s += str(v)           
490        msg += '\t'+ s + '\n'
491    equivVarList = list(set(indepVarList).union(set(depVarList)))
492    #print 'equivVarList',equivVarList
493    inboth = set(fixedVarList).intersection(set(equivVarList))
494    if len(inboth) > 0:
495        msg += "\nError! The following variables are used in both Equivalence and Fixed constraints:\n"
496        s = ''
497        for var in inboth:
498            if s != "": s+= ", "
499            s += str(v)
500        msg += '\t'+ s + '\n'
501
502    # scan through parameters in each relationship. Are all varied? If only some are
503    # varied, create an error message.
504    # If all are varied and this is a constraint equation, then set VaryFree flag
505    # so that newly created relationships will be varied
506    for group,varlist in zip(groups,parmlist):
507        if len(varlist) == 1: continue
508        VaryFree = False
509        for rel in group:
510            varied = 0
511            notvaried = ''
512            for var in constrDict[rel]:
513                if var in varyList:
514                    varied += 1
515                else:
516                    if notvaried: notvaried += ', '
517                    notvaried += var
518                if var in fixedVarList:
519                    msg += '\nError: parameter '+var+" is Fixed and used in a constraint:\n\t"
520                    msg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
521                if var in equivVarList:
522                    msg += '\nError: parameter '+var+" is Equivalenced and used in a constraint:\n\t"
523                    msg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
524            if varied > 0 and varied != len(constrDict[rel]):
525                msg += "\nNot all variables refined in constraint:\n\t"
526                msg += FormatConstraint(constrDict[rel],fixedList[rel])
527                msg += '\nNot refined: ' + notvaried + '\n'
528            if fixedList[rel] is not None and varied > 0:
529                VaryFree = True
530    # if there were errors found, go no farther
531    if msg:
532        print ' *** ERROR in constraint definitions! ***'
533        print msg
534        raise Exception
535               
536    # now process each group and create the relations that are needed to form
537    # non-singular square matrix
538    for group,varlist in zip(groups,parmlist):
539        if len(varlist) == 1: continue
540        arr = MakeArray(constrDict, group, varlist)
541        constrArr = FillInMissingRelations(arr,len(group))
542        mapvar = []
543        group = group[:]
544        # scan through all generated and input variables, add to the varied list
545        # all the new parameters where VaryFree has been set or where all the
546        # dependent parameters are varied. Check again for inconsistent variable use
547        # for new variables -- where varied and unvaried parameters get grouped
548        # together. I don't think this can happen when not flagged before, but
549        # it does not hurt to check again.
550        for i in range(len(varlist)):
551            varied = 0
552            notvaried = ''
553            if len(group) > 0:
554                rel = group.pop(0)
555                fixedval = fixedList[rel]
556                for var in constrDict[rel]:
557                    if var in varyList:
558                        varied += 1
559                    else:
560                        if notvaried: notvaried += ', '
561                        notvaried += var
562            else:
563                fixedval = None
564            if fixedval is None:
565                varname = paramPrefix + str(consNum)
566                mapvar.append(varname)
567                consNum += 1
568                if VaryFree or varied > 0:
569                    varyList.append(varname)
570            else:
571                mapvar.append(fixedval)
572            if varied > 0 and notvaried != '':
573                msg += "\nNot all variables refined in generated constraint\n\t"
574                msg += '\nNot refined: ' + notvaried + '\n'
575        dependentParmList.append(varlist)
576        arrayList.append(constrArr)
577        invarrayList.append(np.linalg.inv(constrArr))
578        indParmList.append(mapvar)
579    if msg:
580        print ' *** ERROR in constraint definitions! ***'
581        print msg
582        print VarRemapShow(varyList)
583        raise Exception
584    # setup dictionary containing the fixed values
585    global fixedDict
586    # key is original ascii string, value is float
587    for fixedval in fixedList:
588        if fixedval:
589            fixedDict[fixedval] = float(fixedval)
590
591    if debug: # on debug, show what is parsed & generated, semi-readable
592        print 50*'-'
593        for group,varlist,multarr,inv,mapvar in zip(groups,parmlist,arrayList,invarrayList,indParmList):
594            print '\n*** relation(s) in group:',group,'\tvars in group:',varlist
595            print 'new parameters:', mapvar
596            print 'Input relationship matrix'
597            print multarr[:len(group)]
598            added = len(group) - len(varlist)
599            if added < 0:
600                print 'added relationship rows'
601                print multarr[added:]
602            print 'Inverse relationship matrix'
603            print inv
604
605def StoreEquivalence(independentVar,dependentList):
606    '''Takes a list of dependent parameter(s) and stores their
607    relationship to a single independent parameter (independentVar)
608    input:
609       independentVar -- name of parameter that will set others (may
610         be varied)
611       dependentList -- list of parameter that will set from
612         independentVar (may not be varied). Each item can be a parameter
613         name or a tuple containing a name and multiplier:
614         ('parm1',('parm2',.5),)
615    '''
616   
617    global dependentParmList,arrayList,invarrayList,indParmList
618    mapList = []
619    multlist = []
620    for var in dependentList:
621        if isinstance(var, basestring):
622            mult = 1.0
623        elif len(var) == 2:
624            var,mult = var
625        else:
626            raise Exception, "Cannot parse "+repr(var) + " as var or (var,multiplier)"
627        mapList.append(var)
628        multlist.append(tuple((mult,)))
629    # added relationships to stored values
630    arrayList.append(None)
631    invarrayList.append(np.array(multlist))
632    indParmList.append(tuple((independentVar,)))
633    dependentParmList.append(mapList)
634    return
635
636def GetDependentVars():
637    '''Return a list of dependent variables: e.g. variables that are
638    constrained in terms of other variables'''
639    dependentVars = []
640    global dependentParmList
641    for lst in dependentParmList:
642        for itm in lst: dependentVars.append(itm)
643    return dependentVars
644
645def GetIndependentVars():
646    '''Return a list of independent variables: e.g. variables that are
647    created by constraints of other variables'''
648    independentVars = []
649    global indParmList,fixedDict
650    for lst in indParmList:
651        for name in lst:
652            if name in fixedDict: continue
653            independentVars.append(name)
654    return independentVars
655
656def PrintIndependentVars(parmDict,varyList,sigDict,PrintAll=False):
657    '''Print the values and uncertainties on the independent variables'''
658    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
659    printlist = []
660    mapvars = GetIndependentVars()
661    for i,name in enumerate(mapvars):
662        if name in fixedDict: continue
663        if PrintAll or name in varyList:
664            sig = sigDict.get(name)
665            printlist.append([name,parmDict[name],sig])
666    if len(printlist) == 0: return
667    s1 = ''
668    print 130*'-'
669    print "Variables generated by constraints"
670    printlist.append(3*[None])
671    for name,val,esd in printlist:
672        if len(s1) > 110 or name is None:
673            print
674            print s1
675            print s2
676            print s3
677            s1 = ''
678            if name is None: break
679        if s1 == "":
680            s1 = ' name  :'
681            s2 = ' value :'
682            s3 = ' sig   :'
683        s1 += '%12s' % (name)
684        s2 += '%12.6f' % (val)
685        if esd is None:
686            s3 += '%12s' % ('n/a')
687        else:   
688            s3 += '%12.6f' % (esd)
689
690def ComputeDepESD(covMatrix,varyList,parmDict):
691    '''Compute uncertainties for dependent parameters from independent ones
692    returns a dictionary containing the esd values for dependent parameters
693    '''
694    sigmaDict = {}
695    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
696        #if invmultarr is None: continue # probably not needed
697        valuelist = [parmDict[var] for var in mapvars]
698        # get the v-covar matrix for independent parameters
699        vcov = np.zeros((len(mapvars),len(mapvars)))
700        for i1,name1 in enumerate(mapvars):
701            if name1 not in varyList: continue
702            iv1 = varyList.index(name1)
703            for i2,name2 in enumerate(mapvars):
704                if name2 not in varyList: continue
705                iv2 = varyList.index(name2)
706                vcov[i1][i2] = covMatrix[iv1][iv2]
707        # vec is the vector that multiplies each of the independent values
708        for v,vec in zip(varlist,invmultarr):
709            sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec)))
710    return sigmaDict
711
712def FormatConstraint(RelDict,RelVal):
713    '''Formats a Constraint or Function for use in a convenient way'''
714    linelen = 45
715    s = [""]
716    for var,val in RelDict.items():
717        if len(s[-1]) > linelen: s.append(' ')
718        m = val
719        if s[-1] != "" and m >= 0:
720            s[-1] += ' + '
721        elif s[-1] != "":
722            s[-1] += ' - '
723            m = abs(m)
724        s[-1] += '%.3f*%s '%(m,var)
725    if len(s[-1]) > linelen: s.append(' ')
726    if RelVal is None:
727        s[-1] += ' = New variable'
728    else:
729        s[-1] += ' = ' + RelVal
730    s1 = ''
731    for s2 in s:
732        if s1 != '': s1 += '\n\t'
733        s1 += s2
734    return s1
735
736def VarRemapShow(varyList):
737    '''List out the saved relationships.
738    Returns a string containing the details.
739    '''
740    s = ''
741    if len(fixedVarList) > 0:
742        s += 'Fixed Variables:\n'
743        for v in fixedVarList:
744            s += '    ' + v + '\n'
745    s += 'Variable mapping relations:\n'
746    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
747    for varlist,mapvars,multarr,invmultarr in zip(
748        dependentParmList,indParmList,arrayList,invarrayList):
749        for i,mv in enumerate(mapvars):
750            if multarr is None:
751                s += '  ' + str(mv) + ' is equivalent to parameter(s): '
752                j = 0
753                for v,m in zip(varlist,invmultarr):
754                    #print v,m[0]
755                    if j > 0: s += '  & '
756                    j += 1
757                    s += str(v)
758                    if m != 1:
759                        s += " / " + str(m[0])                       
760                s += '\n'
761                continue
762            s += '  %s = ' % mv
763            j = 0
764            for m,v in zip(multarr[i,:],varlist):
765                if m == 0: continue
766                if j > 0: s += ' + '
767                j += 1
768                s += '(%s * %s)' % (m,v)
769            if mv in varyList: s += ' VARY'
770            s += '\n'
771    s += 'Inverse variable mapping relations:\n'
772    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
773        for i,mv in enumerate(varlist):
774            s += '  %s = ' % mv
775            j = 0
776            for m,v in zip(invmultarr[i,:],mapvars):
777                if m == 0: continue
778                if j > 0: s += ' + '
779                j += 1
780                s += '(%s * %s)' % (m,v)
781            s += '\n'
782    return s
783
784def Dict2Deriv(varyList,derivDict,dMdv):
785    '''Compute derivatives for Independent Parameters from the
786    derivatives for the original parameters
787    '''
788    global dependentParmList,arrayList,invarrayList,indParmList,invarrayList
789    for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
790        for i,name in enumerate(mapvars):
791            # grouped variables: need to add in the derv. w/r
792            # dependent variables to the independent ones
793            if name not in varyList: continue # skip if independent var not varied
794            if multarr is None:
795                for v,m in zip(varlist,invmultarr):
796                    #print 'add derv',v,'/',m[0],'to derv',name
797                    if m == 0: continue
798                    dMdv[varyList.index(name)] += derivDict[v] / m[0]
799            else:
800                for v,m in zip(varlist,multarr[i,:]):
801                    #print 'add derv',v,'*',m,'to derv',name
802                    if m == 0: continue
803                    dMdv[varyList.index(name)] += m * derivDict[v]
804
805def Map2Dict(parmDict,varyList):
806    '''Create (or update) the Independent Parameters from the original
807    set of Parameters
808
809    Removes dependent variables from the varyList
810
811    This should be done once, before any variable refinement is done
812    to complete the parameter dictionary with the Independent Parameters
813    '''
814    # process the Independent vars: remove dependent ones from varylist
815    # and then compute values for the Independent ones from their dependents
816    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
817    for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
818        for item in varlist:
819            try:
820                varyList.remove(item)
821            except ValueError:
822                pass
823        if multarr is None: continue
824        valuelist = [parmDict[var] for var in varlist]
825        parmDict.update(zip(mapvars,
826                            np.dot(multarr,np.array(valuelist)))
827                        )
828    # now remove fixed variables from the varyList
829    global fixedVarList
830    for item in fixedVarList:
831        try:
832            varyList.remove(item)
833        except ValueError:
834            pass
835    # Set constrained parameters to their fixed values
836    parmDict.update(fixedDict)
837
838def Dict2Map(parmDict,varyList):
839    '''Convert the remapped values back to the original parameters
840   
841    This should be done to apply constraints to parameter values (use
842    Map2Dict once first). It should be done as part of the Model function
843    evaluation, before any computation is done
844    '''
845    #I think this needs fixing to update parmDict with new values
846    #   from the constraints based on what is in varyList - RVD. Don't think so -- BHT
847    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
848    # reset fixed values (should not be needed, but very quick)
849    # - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended
850    # not needed, but also not a problem - BHT
851    parmDict.update(fixedDict)
852    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
853        #if invmultarr is None: continue
854        valuelist = [parmDict[var] for var in mapvars]
855        parmDict.update(zip(varlist,
856                            np.dot(invmultarr,np.array(valuelist)))
857                        )
858
859#======================================================================
860# internal routines follow (these routines are unlikely to be called
861# from outside the module)
862def FillInMissingRelations(arrayin,nvars):
863    '''Fill in unspecified rows in array with non-colinear unit vectors
864    arrayin is a square array, where only the first nvars rows are defined.
865   
866    Returns a new array where all rows are defined
867
868    Throws an exception if there is no non-signular result
869    (likely because two or more input rows are co-linear)
870    '''
871    a = arrayin.copy()
872    n = nvars
873    # fill in the matrix with basis vectors not colinear with any of the starting set
874    for i in range(n,len(a)):
875        try:
876            a[n] = VectorProduct(a[:n])
877        except:
878            raise Exception,"VectorProduct failed, singular input?"
879        n += 1
880    # use the G-S algorithm to compute a complete set of unit vectors orthogonal
881    # to the first in array
882    gs = GramSchmidtOrtho(a) 
883    n = nvars
884    # now replace the created vectors with the ones least correlated to the
885    # initial set
886    vlist = [v for v in gs[1:]] # drop the first row
887    for j in range(n,len(a)):
888        minavg = None
889        droplist = []
890        for k in range(len(vlist)):
891            v = vlist[k]
892            avgcor = 0
893            for i in range(n):
894                cor = np.dot(a[i],vlist[k])/(np.linalg.norm(a[i]) * np.linalg.norm(vlist[k]) )
895                if np.allclose(cor, 1.0):
896                    droplist.append(k) # found a vector co-linear to one we already have
897                                       # don't need it again
898                    #vlist.pop(k)
899                    break 
900                avgcor += cor
901            else:
902                if minavg == None:
903                    minavg = abs(avgcor/n)
904                    minnum = k
905                elif abs(avgcor/n) < minavg:
906                    minavg = abs(avgcor/n)
907                    minnum = k
908        if minavg == None:
909            raise Exception,"Failed to find a non-colinear G-S vector for row %d. Should not happen!" % n
910        a[j] = vlist[minnum]
911        droplist.append(minnum) # don't need to use this vector again
912        for i in sorted(droplist,reverse=True):
913            vlist.pop(i) 
914        n += 1
915    return a
916
917
918def MakeArray(constrDict, group, varlist):
919    """Convert the multipliers in a constraint group to an array of
920    coefficients and place in a square numpy array, adding empty rows as
921    needed.
922
923    Throws an exception if some sanity checks fail.
924    """
925    # do some error checks
926    if len(varlist) < len(group): # too many relationships -- no can do
927        raise Exception,'The number of relationships (%d) exceeds the number of parameters (%d):\n\t%s\n\t%s'% (
928            len(group),len(varlist),group,varlist)
929    # put the coefficients into an array
930    multarr = np.zeros(2*[len(varlist),])
931    i = 0
932    for cnum in group:
933        cdict = constrDict[cnum]
934        j = 0
935        for var in varlist:
936            m = cdict.get(var)
937            if m != None:
938                multarr[i,j] = m
939            j += 1
940        i += 1
941    return multarr
942
943def GramSchmidtOrtho(arrayin):
944    '''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to
945    find orthonormal unit vectors relative to first row.
946    input:
947       arrayin: a 2-D non-signular square array
948    returns:
949       a orthonormal set of unit vectors as a square array
950    '''
951    def proj(a,b):
952        'Projection operator'
953        return a*(np.dot(a,b)/np.dot(a,a))
954    a = arrayin.copy()
955    for j in range (len(a)):
956        for i in range(j):
957            a[j] = a[j] - proj(a[i],a[j])
958        if np.allclose(np.linalg.norm(a[j]),0.0):
959            raise Exception,"Singular input to GramSchmidtOrtho"
960        a[j] = a[j]/np.linalg.norm(a[j])
961    return a
962
963def VectorProduct(vl):
964    '''Evaluate the n-dimensional "cross" product of the list of vectors in vl
965    vl can be any length. The length of each vector is arbitrary, but
966    all must be the same length. See http://en.wikipedia.org/wiki/Levi-Civita_symbol
967
968    This appears to return a vector perpendicular to the supplied set.
969
970    Throws an exception if a vector can not be obtained because the input set of
971    vectors is already complete or contains any redundancies.
972   
973    Uses LeviCitvitaMult recursively to obtain all permutations of vector elements
974    '''
975    i = 0
976    newvec = []
977    for e in vl[0]:
978        i += 1
979        tl = [([i,],1),]
980        seq = LeviCitvitaMult(tl ,vl)
981        tsum = 0
982        for terms,prod in seq:
983            tsum += EvalLCterm(terms) * prod
984        newvec.append(tsum)
985    if np.allclose(newvec,0.0):
986        raise Exception,"VectorProduct failed, singular or already complete input"
987    return newvec
988
989def LeviCitvitaMult(tin,vl):
990    '''Recursion formula to compute all permutations of vector element products
991    The first term in each list is the indices of the Levi-Civita symbol, note
992    that if an index is repeated, the value is zero, so the element is not included.
993    The second term is the product of the vector elements
994    '''
995    v = vl[0]
996    vl1 = vl[1:]       
997
998    j = 0
999    tl = []
1000    for e in vl[0]:
1001        j += 1
1002        for ind,vals in tin:
1003            if j in ind: continue
1004            tl.append((ind+[j,],vals*e))
1005    if len(vl1):
1006        return LeviCitvitaMult(tl,vl1)
1007    else:
1008        return tl
1009
1010def EvalLCterm(term):
1011    '''Evaluate the Levi-Civita symbol Epsilon(i,j,k,...)'''
1012    p = 1
1013    for i in range(len(term)):
1014        for j in range(i+1,len(term)):
1015            p *= (term[j] - term[i])
1016            if p == 0: return 0
1017    return p/abs(p)
1018
1019
1020if __name__ == "__main__":
1021    import sys
1022    parmdict = {}
1023    constrDict = [
1024        {'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},
1025        {'0:0:eA': 0.0},
1026        {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
1027        {'1::C(10,0,1)': 1.0, '2::C(10,0,1)': 1.0},
1028        {'1::AUiso:0': 1.0, '0::AUiso:0': 1.0},
1029        {'0::A0': 0.0}
1030        ]
1031    fixedList = ['5.0', '0', None, None, '1.0', '0']
1032    StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), ))
1033    #StoreEquivalence('1::atomx:3',('2::atomx:3', ('2::atomz:3',2,), )) # error: dependent & independent vars mixed
1034    #StoreEquivalence('1::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) # error: dependent vars repeated
1035    #StoreEquivalence('0:1:eA',('0:0:eA',)) # error: equiv & fixed
1036    #StoreEquivalence('0:99:Scale',('0:12:Scale',)) # error: equiv & constrained
1037    #StoreEquivalence('0:12:Scale',('0:99:Scale',)) # error: equiv & constrained
1038    varylist = ['2::atomx:3',
1039                '2::C(10,6,1)', '1::C(10,6,1)',
1040                '2::atomy:3', '2::atomz:3',
1041                '0:12:Scale', '0:11:Scale', '0:14:Scale', '0:13:Scale', '0:0:Scale']
1042    varyList = ['0::A0', '0::AUiso:0', '0::Afrac:1', '0::Afrac:2', '0::Afrac:3', '0::Afrac:4', '0::dAx:5', '0::dAy:5', '0::dAz:5', '0::AUiso:5', ':0:Back:0', ':0:Back:1', ':0:Back:2', ':0:Back:3', ':0:Back:4', ':0:Back:5', ':0:Back:6', ':0:Back:7', ':0:Back:8', ':0:Back:9', ':0:Back:10', ':0:Back:11', ':0:U', ':0:V', ':0:W', ':0:X', ':0:Y', ':0:Scale', ':0:DisplaceX', ':0:DisplaceY']
1043    constrDict = [
1044        {'0::Afrac:4': 24.0, '0::Afrac:1': 16.0, '0::Afrac:3': 24.0, '0::Afrac:2': 16.0},
1045        {'0::Afrac:1': 1.0, '0::Afrac:2': 1.0},
1046        {'0::Afrac:4': 1.0, '0::Afrac:3': 1.0}]
1047    fixedList = ['40.0', '1.0', '1.0']
1048
1049    errmsg, warnmsg = CheckConstraints(varylist,constrDict,fixedList)
1050    if errmsg:
1051        print "*** Error ********************"
1052        print errmsg
1053    if warnmsg:
1054        print "*** Warning ********************"
1055        print warnmsg
1056    if errmsg or warnmsg:
1057        sys.exit()
1058    groups,parmlist = GroupConstraints(constrDict)
1059    GenerateConstraints(groups,parmlist,varylist,constrDict,fixedList)
1060    print VarRemapShow(varylist)
1061    parmdict.update( {
1062        '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,
1063        '0:0:eA': 0.0,
1064        '2::C(10,6,1)': 0.2, '1::C(10,6,1)': 0.3,
1065        '1::C(10,0,1)': 0.2, '2::C(10,0,1)': 0.3,
1066        '1::AUiso:0': 0.02, '0::AUiso:0': 0.03,
1067        '0::A0': 0.0,
1068        '2::atomx:3':0.23,'2::atomy:3':-.23, '2::atomz:3':-0.11,
1069        })
1070    print 'parmdict start',parmdict
1071    print 'varylist start',varylist
1072    before = parmdict.copy()
1073    Map2Dict(parmdict,varylist)
1074    print 'parmdict before and after Map2Dict'
1075    print '  key / before / after'
1076    for key in sorted(parmdict.keys()):
1077        print '  '+key,'\t',before.get(key),'\t',parmdict[key]
1078    print 'varylist after',varylist
1079    before = parmdict.copy()
1080    Dict2Map(parmdict,varylist)
1081    print 'after Dict2Map'
1082    print '  key / before / after'
1083    for key in sorted(parmdict.keys()):
1084        print '  '+key,'\t',before.get(key),'\t',parmdict[key]
1085    dMdv = len(varylist)*[0]
1086    deriv = {}
1087    for i,v in enumerate(parmdict.keys()): deriv[v]=i
1088    Dict2Deriv(varylist,deriv,dMdv)
Note: See TracBrowser for help on using the repository browser.