source: trunk/GSASIImapvars.py @ 577

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

fix error in reporting constraint errors

File size: 46.6 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 sorted(inboth):
276            if s != "": s+= ", "
277            s += str(var)
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 sorted(set(multdepVarList)):
283            if s != "": s+= ", "
284            s += str(var)           
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 sorted(inboth):
293            if s != "": s+= ", "
294            s += str(var)
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                #s = ''
439                varied = 0
440                notvaried = ''
441                if mv in varyList:
442                    varied += 1
443                else:
444                    if notvaried: notvaried += ', '
445                    notvaried += mv
446                if mv not in indepVarList: indepVarList.append(mv)
447                for v,m in zip(varlist,invmultarr):
448                    #if len(s): s += '  & '
449                    #s += str(v)
450                    #if m != 1:
451                    #    s += " / " + str(m[0])                       
452                    if m == 0: zeromult = True
453                    if v in varyList:
454                        varied += 1
455                    else:
456                        if notvaried: notvaried += ', '
457                        notvaried += v
458                    if v in depVarList:
459                        multdepVarList.append(v)
460                    else:
461                        depVarList.append(v)
462                #print str(mv) + ' is equivalent to parameter(s): '+s
463            if varied > 0 and varied != len(varlist)+1:
464                msg += "\nNot all variables refined in equivalence:\n\t"
465                s = ""
466                for v in varlist:
467                    if s != "": s+= " & "
468                    s += str(v)           
469                msg += str(mv) + " => " + s
470                msg += '\nNot refined: ' + notvaried + '\n'
471            if zeromult:
472                msg += "\nZero multiplier is invalid in equivalence:\n\t"
473                s = ""
474                for v in varlist:
475                    if s != "": s+= " & "
476                    s += str(v)           
477                msg += str(mv) + " => " + s + '\n'
478
479    #print 'indepVarList',indepVarList
480    #print 'depVarList',depVarList
481    # check for errors:
482    inboth = set(indepVarList).intersection(set(depVarList))
483    if len(inboth) > 0:
484        msg += "\nThe following parameters(s) are used as both dependent and independent variables in Equivalence relations:\n"
485        s = ''
486        for var in sorted(inboth):
487            if s != "": s+= ", "
488            s += str(var)
489        msg += '\t'+ s + '\n'
490    if len(multdepVarList) > 0:
491        msg += "\nThe following parameters(s) are used in multiple Equivalence relations as dependent variables:\n"
492        s = ''
493        for var in sorted(set(multdepVarList)):
494            if s != "": s+= ", "
495            s += str(var)           
496        msg += '\t'+ s + '\n'
497    equivVarList = list(set(indepVarList).union(set(depVarList)))
498    #print 'equivVarList',equivVarList
499    inboth = set(fixedVarList).intersection(set(equivVarList))
500    if len(inboth) > 0:
501        msg += "\nError! The following variables are used in both Equivalence and Fixed constraints:\n"
502        s = ''
503        for var in sorted(inboth):
504            if s != "": s+= ", "
505            s += str(var)
506        msg += '\t'+ s + '\n'
507
508    # scan through parameters in each relationship. Are all varied? If only some are
509    # varied, create an error message.
510    # If all are varied and this is a constraint equation, then set VaryFree flag
511    # so that newly created relationships will be varied
512    for group,varlist in zip(groups,parmlist):
513        if len(varlist) == 1: continue
514        VaryFree = False
515        for rel in group:
516            varied = 0
517            notvaried = ''
518            for var in constrDict[rel]:
519                if var in varyList:
520                    varied += 1
521                else:
522                    if notvaried: notvaried += ', '
523                    notvaried += var
524                if var in fixedVarList:
525                    msg += '\nError: parameter '+var+" is Fixed and used in a constraint:\n\t"
526                    msg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
527                if var in equivVarList:
528                    msg += '\nError: parameter '+var+" is Equivalenced and used in a constraint:\n\t"
529                    msg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
530            if varied > 0 and varied != len(constrDict[rel]):
531                msg += "\nNot all variables refined in constraint:\n\t"
532                msg += FormatConstraint(constrDict[rel],fixedList[rel])
533                msg += '\nNot refined: ' + notvaried + '\n'
534            if fixedList[rel] is not None and varied > 0:
535                VaryFree = True
536    # if there were errors found, go no farther
537    if msg:
538        print ' *** ERROR in constraint definitions! ***'
539        print msg
540        raise Exception
541               
542    # now process each group and create the relations that are needed to form
543    # non-singular square matrix
544    for group,varlist in zip(groups,parmlist):
545        if len(varlist) == 1: continue
546        arr = MakeArray(constrDict, group, varlist)
547        constrArr = FillInMissingRelations(arr,len(group))
548        mapvar = []
549        group = group[:]
550        # scan through all generated and input variables, add to the varied list
551        # all the new parameters where VaryFree has been set or where all the
552        # dependent parameters are varied. Check again for inconsistent variable use
553        # for new variables -- where varied and unvaried parameters get grouped
554        # together. I don't think this can happen when not flagged before, but
555        # it does not hurt to check again.
556        for i in range(len(varlist)):
557            varied = 0
558            notvaried = ''
559            if len(group) > 0:
560                rel = group.pop(0)
561                fixedval = fixedList[rel]
562                for var in constrDict[rel]:
563                    if var in varyList:
564                        varied += 1
565                    else:
566                        if notvaried: notvaried += ', '
567                        notvaried += var
568            else:
569                fixedval = None
570            if fixedval is None:
571                varname = paramPrefix + str(consNum)
572                mapvar.append(varname)
573                consNum += 1
574                if VaryFree or varied > 0:
575                    varyList.append(varname)
576            else:
577                mapvar.append(fixedval)
578            if varied > 0 and notvaried != '':
579                msg += "\nNot all variables refined in generated constraint\n\t"
580                msg += '\nNot refined: ' + notvaried + '\n'
581        dependentParmList.append(varlist)
582        arrayList.append(constrArr)
583        invarrayList.append(np.linalg.inv(constrArr))
584        indParmList.append(mapvar)
585    if msg:
586        print ' *** ERROR in constraint definitions! ***'
587        print msg
588        print VarRemapShow(varyList)
589        raise Exception
590    # setup dictionary containing the fixed values
591    global fixedDict
592    # key is original ascii string, value is float
593    for fixedval in fixedList:
594        if fixedval:
595            fixedDict[fixedval] = float(fixedval)
596
597    if debug: # on debug, show what is parsed & generated, semi-readable
598        print 50*'-'
599        for group,varlist,multarr,inv,mapvar in zip(groups,parmlist,arrayList,invarrayList,indParmList):
600            print '\n*** relation(s) in group:',group,'\tvars in group:',varlist
601            print 'new parameters:', mapvar
602            print 'Input relationship matrix'
603            print multarr[:len(group)]
604            added = len(group) - len(varlist)
605            if added < 0:
606                print 'added relationship rows'
607                print multarr[added:]
608            print 'Inverse relationship matrix'
609            print inv
610
611def StoreEquivalence(independentVar,dependentList):
612    '''Takes a list of dependent parameter(s) and stores their
613    relationship to a single independent parameter (independentVar)
614    input:
615       independentVar -- name of parameter that will set others (may
616         be varied)
617       dependentList -- list of parameter that will set from
618         independentVar (may not be varied). Each item can be a parameter
619         name or a tuple containing a name and multiplier:
620         ('parm1',('parm2',.5),)
621    '''
622   
623    global dependentParmList,arrayList,invarrayList,indParmList
624    mapList = []
625    multlist = []
626    for var in dependentList:
627        if isinstance(var, basestring):
628            mult = 1.0
629        elif len(var) == 2:
630            var,mult = var
631        else:
632            raise Exception, "Cannot parse "+repr(var) + " as var or (var,multiplier)"
633        mapList.append(var)
634        multlist.append(tuple((mult,)))
635    # added relationships to stored values
636    arrayList.append(None)
637    invarrayList.append(np.array(multlist))
638    indParmList.append(tuple((independentVar,)))
639    dependentParmList.append(mapList)
640    return
641
642def GetDependentVars():
643    '''Return a list of dependent variables: e.g. variables that are
644    constrained in terms of other variables'''
645    dependentVars = []
646    global dependentParmList
647    for lst in dependentParmList:
648        for itm in lst: dependentVars.append(itm)
649    return dependentVars
650
651def GetIndependentVars():
652    '''Return a list of independent variables: e.g. variables that are
653    created by constraints of other variables'''
654    independentVars = []
655    global indParmList,fixedDict
656    for lst in indParmList:
657        for name in lst:
658            if name in fixedDict: continue
659            independentVars.append(name)
660    return independentVars
661
662def PrintIndependentVars(parmDict,varyList,sigDict,PrintAll=False):
663    '''Print the values and uncertainties on the independent variables'''
664    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
665    printlist = []
666    mapvars = GetIndependentVars()
667    for i,name in enumerate(mapvars):
668        if name in fixedDict: continue
669        if PrintAll or name in varyList:
670            sig = sigDict.get(name)
671            printlist.append([name,parmDict[name],sig])
672    if len(printlist) == 0: return
673    s1 = ''
674    print 130*'-'
675    print "Variables generated by constraints"
676    printlist.append(3*[None])
677    for name,val,esd in printlist:
678        if len(s1) > 110 or name is None:
679            print
680            print s1
681            print s2
682            print s3
683            s1 = ''
684            if name is None: break
685        if s1 == "":
686            s1 = ' name  :'
687            s2 = ' value :'
688            s3 = ' sig   :'
689        s1 += '%12s' % (name)
690        s2 += '%12.6f' % (val)
691        if esd is None:
692            s3 += '%12s' % ('n/a')
693        else:   
694            s3 += '%12.6f' % (esd)
695
696def ComputeDepESD(covMatrix,varyList,parmDict):
697    '''Compute uncertainties for dependent parameters from independent ones
698    returns a dictionary containing the esd values for dependent parameters
699    '''
700    sigmaDict = {}
701    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
702        #if invmultarr is None: continue # probably not needed
703        valuelist = [parmDict[var] for var in mapvars]
704        # get the v-covar matrix for independent parameters
705        vcov = np.zeros((len(mapvars),len(mapvars)))
706        for i1,name1 in enumerate(mapvars):
707            if name1 not in varyList: continue
708            iv1 = varyList.index(name1)
709            for i2,name2 in enumerate(mapvars):
710                if name2 not in varyList: continue
711                iv2 = varyList.index(name2)
712                vcov[i1][i2] = covMatrix[iv1][iv2]
713        # vec is the vector that multiplies each of the independent values
714        for v,vec in zip(varlist,invmultarr):
715            sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec)))
716    return sigmaDict
717
718def FormatConstraint(RelDict,RelVal):
719    '''Formats a Constraint or Function for use in a convenient way'''
720    linelen = 45
721    s = [""]
722    for var,val in RelDict.items():
723        if len(s[-1]) > linelen: s.append(' ')
724        m = val
725        if s[-1] != "" and m >= 0:
726            s[-1] += ' + '
727        elif s[-1] != "":
728            s[-1] += ' - '
729            m = abs(m)
730        s[-1] += '%.3f*%s '%(m,var)
731    if len(s[-1]) > linelen: s.append(' ')
732    if RelVal is None:
733        s[-1] += ' = New variable'
734    else:
735        s[-1] += ' = ' + RelVal
736    s1 = ''
737    for s2 in s:
738        if s1 != '': s1 += '\n\t'
739        s1 += s2
740    return s1
741
742def VarRemapShow(varyList):
743    '''List out the saved relationships.
744    Returns a string containing the details.
745    '''
746    s = ''
747    if len(fixedVarList) > 0:
748        s += 'Fixed Variables:\n'
749        for v in fixedVarList:
750            s += '    ' + v + '\n'
751    s += 'Variable mapping relations:\n'
752    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
753    for varlist,mapvars,multarr,invmultarr in zip(
754        dependentParmList,indParmList,arrayList,invarrayList):
755        for i,mv in enumerate(mapvars):
756            if multarr is None:
757                s += '  ' + str(mv) + ' is equivalent to parameter(s): '
758                j = 0
759                for v,m in zip(varlist,invmultarr):
760                    #print v,m[0]
761                    if j > 0: s += '  & '
762                    j += 1
763                    s += str(v)
764                    if m != 1:
765                        s += " / " + str(m[0])                       
766                s += '\n'
767                continue
768            s += %s = ' % mv
769            j = 0
770            for m,v in zip(multarr[i,:],varlist):
771                if m == 0: continue
772                if j > 0: s += ' + '
773                j += 1
774                s += '(%s * %s)' % (m,v)
775            if mv in varyList: s += ' VARY'
776            s += '\n'
777    s += 'Inverse variable mapping relations:\n'
778    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
779        for i,mv in enumerate(varlist):
780            s += %s = ' % mv
781            j = 0
782            for m,v in zip(invmultarr[i,:],mapvars):
783                if m == 0: continue
784                if j > 0: s += ' + '
785                j += 1
786                s += '(%s * %s)' % (m,v)
787            s += '\n'
788    return s
789
790def Dict2Deriv(varyList,derivDict,dMdv):
791    '''Compute derivatives for Independent Parameters from the
792    derivatives for the original parameters
793    '''
794    global dependentParmList,arrayList,invarrayList,indParmList,invarrayList
795    for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
796        for i,name in enumerate(mapvars):
797            # grouped variables: need to add in the derv. w/r
798            # dependent variables to the independent ones
799            if name not in varyList: continue # skip if independent var not varied
800            if multarr is None:
801                for v,m in zip(varlist,invmultarr):
802                    #print 'add derv',v,'/',m[0],'to derv',name
803                    if m == 0: continue
804                    dMdv[varyList.index(name)] += derivDict[v] / m[0]
805            else:
806                for v,m in zip(varlist,multarr[i,:]):
807                    #print 'add derv',v,'*',m,'to derv',name
808                    if m == 0: continue
809                    dMdv[varyList.index(name)] += m * derivDict[v]
810
811def Map2Dict(parmDict,varyList):
812    '''Create (or update) the Independent Parameters from the original
813    set of Parameters
814
815    Removes dependent variables from the varyList
816
817    This should be done once, before any variable refinement is done
818    to complete the parameter dictionary with the Independent Parameters
819    '''
820    # process the Independent vars: remove dependent ones from varylist
821    # and then compute values for the Independent ones from their dependents
822    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
823    for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
824        for item in varlist:
825            try:
826                varyList.remove(item)
827            except ValueError:
828                pass
829        if multarr is None: continue
830        valuelist = [parmDict[var] for var in varlist]
831        parmDict.update(zip(mapvars,
832                            np.dot(multarr,np.array(valuelist)))
833                        )
834    # now remove fixed variables from the varyList
835    global fixedVarList
836    for item in fixedVarList:
837        try:
838            varyList.remove(item)
839        except ValueError:
840            pass
841    # Set constrained parameters to their fixed values
842    parmDict.update(fixedDict)
843
844def Dict2Map(parmDict,varyList):
845    '''Convert the remapped values back to the original parameters
846   
847    This should be done to apply constraints to parameter values (use
848    Map2Dict once first). It should be done as part of the Model function
849    evaluation, before any computation is done
850    '''
851    #I think this needs fixing to update parmDict with new values
852    #   from the constraints based on what is in varyList - RVD. Don't think so -- BHT
853    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
854    # reset fixed values (should not be needed, but very quick)
855    # - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended
856    # not needed, but also not a problem - BHT
857    parmDict.update(fixedDict)
858    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
859        #if invmultarr is None: continue
860        valuelist = [parmDict[var] for var in mapvars]
861        parmDict.update(zip(varlist,
862                            np.dot(invmultarr,np.array(valuelist)))
863                        )
864
865#======================================================================
866# internal routines follow (these routines are unlikely to be called
867# from outside the module)
868def FillInMissingRelations(arrayin,nvars):
869    '''Fill in unspecified rows in array with non-colinear unit vectors
870    arrayin is a square array, where only the first nvars rows are defined.
871   
872    Returns a new array where all rows are defined
873
874    Throws an exception if there is no non-signular result
875    (likely because two or more input rows are co-linear)
876    '''
877    a = arrayin.copy()
878    n = nvars
879    # fill in the matrix with basis vectors not colinear with any of the starting set
880    for i in range(n,len(a)):
881        try:
882            a[n] = VectorProduct(a[:n])
883        except:
884            raise Exception,"VectorProduct failed, singular input?"
885        n += 1
886    # use the G-S algorithm to compute a complete set of unit vectors orthogonal
887    # to the first in array
888    gs = GramSchmidtOrtho(a) 
889    n = nvars
890    # now replace the created vectors with the ones least correlated to the
891    # initial set
892    vlist = [v for v in gs[1:]] # drop the first row
893    for j in range(n,len(a)):
894        minavg = None
895        droplist = []
896        for k in range(len(vlist)):
897            v = vlist[k]
898            avgcor = 0
899            for i in range(n):
900                cor = np.dot(a[i],vlist[k])/(np.linalg.norm(a[i]) * np.linalg.norm(vlist[k]) )
901                if np.allclose(cor, 1.0):
902                    droplist.append(k) # found a vector co-linear to one we already have
903                                       # don't need it again
904                    #vlist.pop(k)
905                    break 
906                avgcor += cor
907            else:
908                if minavg == None:
909                    minavg = abs(avgcor/n)
910                    minnum = k
911                elif abs(avgcor/n) < minavg:
912                    minavg = abs(avgcor/n)
913                    minnum = k
914        if minavg == None:
915            raise Exception,"Failed to find a non-colinear G-S vector for row %d. Should not happen!" % n
916        a[j] = vlist[minnum]
917        droplist.append(minnum) # don't need to use this vector again
918        for i in sorted(droplist,reverse=True):
919            vlist.pop(i) 
920        n += 1
921    return a
922
923
924def MakeArray(constrDict, group, varlist):
925    """Convert the multipliers in a constraint group to an array of
926    coefficients and place in a square numpy array, adding empty rows as
927    needed.
928
929    Throws an exception if some sanity checks fail.
930    """
931    # do some error checks
932    if len(varlist) < len(group): # too many relationships -- no can do
933        raise Exception,'The number of relationships (%d) exceeds the number of parameters (%d):\n\t%s\n\t%s'% (
934            len(group),len(varlist),group,varlist)
935    # put the coefficients into an array
936    multarr = np.zeros(2*[len(varlist),])
937    i = 0
938    for cnum in group:
939        cdict = constrDict[cnum]
940        j = 0
941        for var in varlist:
942            m = cdict.get(var)
943            if m != None:
944                multarr[i,j] = m
945            j += 1
946        i += 1
947    return multarr
948
949def GramSchmidtOrtho(arrayin):
950    '''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to
951    find orthonormal unit vectors relative to first row.
952    input:
953       arrayin: a 2-D non-signular square array
954    returns:
955       a orthonormal set of unit vectors as a square array
956    '''
957    def proj(a,b):
958        'Projection operator'
959        return a*(np.dot(a,b)/np.dot(a,a))
960    a = arrayin.copy()
961    for j in range (len(a)):
962        for i in range(j):
963            a[j] = a[j] - proj(a[i],a[j])
964        if np.allclose(np.linalg.norm(a[j]),0.0):
965            raise Exception,"Singular input to GramSchmidtOrtho"
966        a[j] = a[j]/np.linalg.norm(a[j])
967    return a
968
969def VectorProduct(vl):
970    '''Evaluate the n-dimensional "cross" product of the list of vectors in vl
971    vl can be any length. The length of each vector is arbitrary, but
972    all must be the same length. See http://en.wikipedia.org/wiki/Levi-Civita_symbol
973
974    This appears to return a vector perpendicular to the supplied set.
975
976    Throws an exception if a vector can not be obtained because the input set of
977    vectors is already complete or contains any redundancies.
978   
979    Uses LeviCitvitaMult recursively to obtain all permutations of vector elements
980    '''
981    i = 0
982    newvec = []
983    for e in vl[0]:
984        i += 1
985        tl = [([i,],1),]
986        seq = LeviCitvitaMult(tl ,vl)
987        tsum = 0
988        for terms,prod in seq:
989            tsum += EvalLCterm(terms) * prod
990        newvec.append(tsum)
991    if np.allclose(newvec,0.0):
992        raise Exception,"VectorProduct failed, singular or already complete input"
993    return newvec
994
995def LeviCitvitaMult(tin,vl):
996    '''Recursion formula to compute all permutations of vector element products
997    The first term in each list is the indices of the Levi-Civita symbol, note
998    that if an index is repeated, the value is zero, so the element is not included.
999    The second term is the product of the vector elements
1000    '''
1001    v = vl[0]
1002    vl1 = vl[1:]       
1003
1004    j = 0
1005    tl = []
1006    for e in vl[0]:
1007        j += 1
1008        for ind,vals in tin:
1009            if j in ind: continue
1010            tl.append((ind+[j,],vals*e))
1011    if len(vl1):
1012        return LeviCitvitaMult(tl,vl1)
1013    else:
1014        return tl
1015
1016def EvalLCterm(term):
1017    '''Evaluate the Levi-Civita symbol Epsilon(i,j,k,...)'''
1018    p = 1
1019    for i in range(len(term)):
1020        for j in range(i+1,len(term)):
1021            p *= (term[j] - term[i])
1022            if p == 0: return 0
1023    return p/abs(p)
1024
1025
1026if __name__ == "__main__":
1027    import sys
1028    parmdict = {}
1029    constrDict = [
1030        {'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},
1031        {'0:0:eA': 0.0},
1032        {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
1033        {'1::C(10,0,1)': 1.0, '2::C(10,0,1)': 1.0},
1034        {'1::AUiso:0': 1.0, '0::AUiso:0': 1.0},
1035        {'0::A0': 0.0}
1036        ]
1037    fixedList = ['5.0', '0', None, None, '1.0', '0']
1038    StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), ))
1039    #StoreEquivalence('1::atomx:3',('2::atomx:3', ('2::atomz:3',2,), )) # error: dependent & independent vars mixed
1040    #StoreEquivalence('1::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) # error: dependent vars repeated
1041    #StoreEquivalence('0:1:eA',('0:0:eA',)) # error: equiv & fixed
1042    #StoreEquivalence('0:99:Scale',('0:12:Scale',)) # error: equiv & constrained
1043    #StoreEquivalence('0:12:Scale',('0:99:Scale',)) # error: equiv & constrained
1044    varylist = ['2::atomx:3',
1045                '2::C(10,6,1)', '1::C(10,6,1)',
1046                '2::atomy:3', '2::atomz:3',
1047                '0:12:Scale', '0:11:Scale', '0:14:Scale', '0:13:Scale', '0:0:Scale']
1048    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']
1049    constrDict = [
1050        {'0::Afrac:4': 24.0, '0::Afrac:1': 16.0, '0::Afrac:3': 24.0, '0::Afrac:2': 16.0},
1051        {'0::Afrac:1': 1.0, '0::Afrac:2': 1.0},
1052        {'0::Afrac:4': 1.0, '0::Afrac:3': 1.0}]
1053    fixedList = ['40.0', '1.0', '1.0']
1054
1055    errmsg, warnmsg = CheckConstraints(varylist,constrDict,fixedList)
1056    if errmsg:
1057        print "*** Error ********************"
1058        print errmsg
1059    if warnmsg:
1060        print "*** Warning ********************"
1061        print warnmsg
1062    if errmsg or warnmsg:
1063        sys.exit()
1064    groups,parmlist = GroupConstraints(constrDict)
1065    GenerateConstraints(groups,parmlist,varylist,constrDict,fixedList)
1066    print VarRemapShow(varylist)
1067    parmdict.update( {
1068        '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,
1069        '0:0:eA': 0.0,
1070        '2::C(10,6,1)': 0.2, '1::C(10,6,1)': 0.3,
1071        '1::C(10,0,1)': 0.2, '2::C(10,0,1)': 0.3,
1072        '1::AUiso:0': 0.02, '0::AUiso:0': 0.03,
1073        '0::A0': 0.0,
1074        '2::atomx:3':0.23,'2::atomy:3':-.23, '2::atomz:3':-0.11,
1075        })
1076    print 'parmdict start',parmdict
1077    print 'varylist start',varylist
1078    before = parmdict.copy()
1079    Map2Dict(parmdict,varylist)
1080    print 'parmdict before and after Map2Dict'
1081    print '  key / before / after'
1082    for key in sorted(parmdict.keys()):
1083        print '  '+key,'\t',before.get(key),'\t',parmdict[key]
1084    print 'varylist after',varylist
1085    before = parmdict.copy()
1086    Dict2Map(parmdict,varylist)
1087    print 'after Dict2Map'
1088    print '  key / before / after'
1089    for key in sorted(parmdict.keys()):
1090        print '  '+key,'\t',before.get(key),'\t',parmdict[key]
1091    dMdv = len(varylist)*[0]
1092    deriv = {}
1093    for i,v in enumerate(parmdict.keys()): deriv[v]=i
1094    Dict2Deriv(varylist,deriv,dMdv)
Note: See TracBrowser for help on using the repository browser.