Source code for GSASIImapvars
# -*- coding: utf-8 -*-
########### SVN repository information ###################
# $Date: 2015-02-28 12:46:37 +0900 (Sat, 28 Feb 2015) $
# $Author: toby $
# $Revision: 1676 $
# $URL: https://subversion.xray.aps.anl.gov/pyGSAS/trunk/GSASIImapvars.py $
# $Id: GSASIImapvars.py 1676 2015-02-28 03:46:37Z toby $
########### SVN repository information ###################
"""
*GSASIImapvars: Parameter constraints*
======================================
Module to implements algebraic contraints, parameter redefinition
and parameter simplification contraints.
Parameter redefinition (new vars) is done by creating one or more relationships
between a set of parameters
::
Mx1 * Px + My1 * Py +...
Mx2 * Px + Mz2 * Pz + ...
where Pj is a parameter name and Mjk is a constant.
Constant constraint Relations can also be supplied in the form of an equation:
::
nx1 * Px + ny1 * Py +... = C1
where Cn is a constant. These equations define an algebraic
constrant.
Parameters can also be "fixed" (held), which prevents them from being refined.
All of the above three cases are input using routines
GroupConstraints and GenerateConstraints. The input consists of a list of
relationship dictionaries:
.. code-block:: python
constrDict = [
{'0:12:Scale': 2.0, '0:14:Scale': 4.0, '0:13:Scale': 3.0, '0:0:Scale': 0.5},
{'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
{'0::A0': 0.0}]
fixedList = ['5.0', None, '0']
Where the dictionary defines the first part of an expression and the corresponding fixedList
item is either None (for parameter redefinition) or the constant value, for a constant
constraint equation. A dictionary that contains a single term defines a variable that
will be fixed (held). The multiplier and the fixedList value in this case are ignored.
Parameters can also be equivalenced or "slaved" to another parameter, such that one
(independent) parameter is equated to several (now dependent) parameters. In
algebraic form this is:
::
P0 = M1 * P1 = M2 * P2 = ...
Thus parameters P0, P1 and P2,... are linearly equivalent. Routine StoreEquivalence is
used to specify these equivalences.
Parameter redefinition (new vars) describes a new, independent, parameter, which
is defined in terms of dependent parameters that are defined in the
Model, while fixed constrained relations effectively reduce the complexity
of the Model by removing a degree of freedom. It is possible for a parameter to appear
in both a parameter redefinition expression and a fixed constraint equation, but a
parameter cannot be used a parameter equivalance cannot be used elsewhere (not fixed,
constrained or redefined). Likewise a fixed parameter cannot be used elsewhere (not
equivalanced, constrained or redefined).
Relationships are grouped so that a set of dependent parameters appear
in only one group (done in routine GroupConstraints.) Note that if a
group contains relationships/equations that involve N dependent
parameters, there must exist N-C new parameters, where C is the number
of contraint equations in the group. Routine GenerateConstraints takes
the output from GroupConstraints and generates the
"missing" relationships and saves that information in the module's
global variables. Each generated parameter is named sequentially using paramPrefix.
A list of parameters that will be varied is specified as input to GenerateConstraints
(varyList). A fixed parameter will simply be removed from this list preventing that
parameter from being varied. Note that all parameters in a constraint relationship
must specified as varied (appear in varyList) or none can be varied. This is
checked in GenerateConstraints. Likewise, if all parameters in a constraint are
not referenced in a refinement, the constraint is ignored, but if some parameters
in a constraint group are not referenced in a refinement, but others are this
constitutes and error.
* When a new variable is created, the variable is assigned the name associated
in the constraint definition or it is assigned a default name of form
``::constr<n>`` (see paramPrefix). The vary setting for variables used in the
constraint are ignored.
Note that any generated "missing" relations are not varied. Only
the input relations can be are varied.
* If all parameters in a fixed constraint equation are varied, the generated "missing"
relations in the group are all varied. This provides the N-C degrees of freedom.
*External Routines*
-------------------
To define a set of constrained and unconstrained relations, one
defines a list of dictionary defining constraint parameters and their
values, a list of fixed values for each constraint and a list of
parameters to be varied. In addition, one uses
:func:`StoreEquivalence` to define parameters that are equivalent. One
can then use :func:`CheckConstraints` to check that the input is
internally consistent and finally :func:`GroupConstraints` and
:func:`GenerateConstraints` to generate the internally used
tables. Routines :func:`Map2Dict` is used to initialize the parameter
dictionary and :func:`Dict2Map`, :func:`Dict2Deriv`, and
:func:`ComputeDepESD` are used to apply constraints. Routine
:func:`VarRemapShow` is used to print out the constraint information,
as stored by :func:`GenerateConstraints`.
:func:`InitVars`
This is optionally used to clear out all defined previously defined constraint information
:func:`StoreEquivalence`
To implement parameter redefinition, one calls StoreEquivalence. This should be called for every set of
equivalence relationships. There is no harm in using StoreEquivalence with the same independent variable:
.. code-block:: python
StoreEquivalence('x',('y',))
StoreEquivalence('x',('z',))
or equivalently
.. code-block:: python
StoreEquivalence('x',('y','z'))
The latter will run more efficiently. Note that mixing independent and dependent variables is
problematic. This is not allowed:
.. code-block:: python
StoreEquivalence('x',('y',))
StoreEquivalence('y',('z',))
Use StoreEquivalence before calling GenerateConstraints or CheckConstraints
:func:`CheckConstraints`
To check that input in internally consistent, use CheckConstraints
:func:`Map2Dict`
To determine values for the parameters created in this module, one
calls Map2Dict. This will not apply contraints.
:func:`Dict2Map`
To take values from the new independent parameters and constraints,
one calls Dict2Map. This will apply contraints.
:func:`Dict2Deriv`
Use Dict2Deriv to determine derivatives on independent parameters
from those on dependent ones
:func:`ComputeDepESD`
Use ComputeDepESD to compute uncertainties on dependent variables
:func:`VarRemapShow`
To show a summary of the parameter remapping, one calls VarRemapShow
*Global Variables*
------------------
dependentParmList:
contains a list by group of lists of
parameters used in the group. Note that parameters listed in
dependentParmList should not be refined as they will not affect
the model
indParmList:
a list of groups of Independent parameters defined in
each group. This contains both parameters used in parameter
redefinitions as well as names of generated new parameters.
fixedVarList:
a list of variables that have been 'fixed'
by defining them as equal to a constant (::var: = 0). Note that
the constant value is ignored at present. These variables are
later removed from varyList which prevents them from being refined.
Unlikely to be used externally.
arrayList:
a list by group of relationship matrices to relate
parameters in dependentParmList to those in indParmList. Unlikely
to be used externally.
invarrayList:
a list by group of relationship matrices to relate
parameters in indParmList to those in dependentParmList. Unlikely
to be used externally.
fixedDict:
a dictionary containing the fixed values corresponding
to parameter equations. The dict key is an ascii string, but the
dict value is a float. Unlikely to be used externally.
*Routines*
----------
Note that parameter names in GSAS-II are strings of form ``<ph>:<hst>:<nam>``
"""
import numpy as np
import GSASIIpath
GSASIIpath.SetVersionNumber("$Revision: 1676 $")
# data used for constraints;
debug = False # turns on printing as constraint input is processed
# note that constraints are stored listed by contraint groups, where each constraint
# group contains those parameters that must be handled together
dependentParmList = [] # contains a list of parameters in each group
# note that parameters listed in dependentParmList should not be refined
arrayList = [] # a list of of relationship matrices
invarrayList = [] # a list of inverse relationship matrices
indParmList = [] # a list of names for the new parameters
fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations
# key is original ascii string, value is float
fixedVarList = [] # List of variables that should not be refined
# prefix for parameter names
paramPrefix = "::constr"
consNum = 0 # number of the next constraint to be created
[docs]def InitVars():
'''Initializes all constraint information'''
global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,consNum
dependentParmList = [] # contains a list of parameters in each group
arrayList = [] # a list of of relationship matrices
invarrayList = [] # a list of inverse relationship matrices
indParmList = [] # a list of names for the new parameters
fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations
consNum = 0 # number of the next constraint to be created
fixedVarList = []
[docs]def VarKeys(constr):
"""Finds the keys in a constraint that represent variables
e.g. eliminates any that start with '_'
:param dict constr: a single constraint entry of form::
{'var1': mult1, 'var2': mult2,... '_notVar': val,...}
(see :func:`GroupConstraints`)
:returns: a list of keys where any keys beginning with '_' are
removed.
"""
return [i for i in constr.keys() if not i.startswith('_')]
[docs]def GroupConstraints(constrDict):
"""divide the constraints into groups that share no parameters.
:param dict constrDict: a list of dicts defining relationships/constraints
::
constrDict = [{<constr1>}, {<constr2>}, ...]
where {<constr1>} is {'var1': mult1, 'var2': mult2,... }
:returns: two lists of lists:
* a list of grouped contraints where each constraint grouped containts a list
of indices for constraint constrDict entries
* a list containing lists of parameter names contained in each group
"""
assignedlist = [] # relationships that have been used
groups = [] # contains a list of grouplists
ParmList = []
for i,consi in enumerate(constrDict):
if i in assignedlist: continue # already in a group, skip
# starting a new group
grouplist = [i,]
assignedlist.append(i)
groupset = set(VarKeys(consi))
changes = True # always loop at least once
while(changes): # loop until we can't find anything to add to the current group
changes = False # but don't loop again unless we find something
for j,consj in enumerate(constrDict):
if j in assignedlist: continue # already in a group, skip
if len(set(VarKeys(consj)) & groupset) > 0: # true if this needs to be added
changes = True
grouplist.append(j)
assignedlist.append(j)
groupset = groupset | set(VarKeys(consj))
group = sorted(grouplist)
varlist = sorted(list(groupset))
groups.append(group)
ParmList.append(varlist)
return groups,ParmList
[docs]def CheckConstraints(varyList,constrDict,fixedList):
'''Takes a list of relationship entries comprising a group of
constraints and checks for inconsistencies such as conflicts in
parameter/variable definitions and or inconsistently varied parameters.
:param list varyList: a list of parameters names that will be varied
:param dict constrDict: a list of dicts defining relationships/constraints
(as created in :func:`GSASIIstrIO.ProcessConstraints` and
documented in :func:`GroupConstraints`)
:param list fixedList: a list of values specifying a fixed value for each
dict in constrDict. Values are either strings that can be converted to
floats or ``None`` if the constraint defines a new parameter rather
than a constant.
:returns: two strings:
* the first lists conflicts internal to the specified constraints
* the second lists conflicts where the varyList specifies some
parameters in a constraint, but not all
If there are no errors, both strings will be empty
'''
import re
global dependentParmList,arrayList,invarrayList,indParmList,consNum
errmsg = ''
warnmsg = ''
fixVlist = []
# process fixed variables (holds)
for cdict in constrDict:
# N.B. No "_" names in holds
if len(cdict) == 1:
fixVlist.append(cdict.keys()[0])
# process equivalences: make a list of dependent and independent vars
# and check for repeated uses (repetition of a parameter as an
# independent var is OK)
indepVarList = []
depVarList = []
multdepVarList = []
for varlist,mapvars,multarr,invmultarr in zip(
dependentParmList,indParmList,arrayList,invarrayList):
if multarr is None: # an equivalence
zeromult = False
for mv in mapvars:
varied = 0
notvaried = ''
if mv in varyList:
varied += 1
else:
if notvaried: notvaried += ', '
notvaried += mv
if mv not in indepVarList: indepVarList.append(mv)
for v,m in zip(varlist,invmultarr):
if m == 0: zeromult = True
if v in varyList:
varied += 1
else:
if notvaried: notvaried += ', '
notvaried += v
if v in depVarList:
multdepVarList.append(v)
else:
depVarList.append(v)
if varied > 0 and varied != len(varlist)+1:
warnmsg += "\nNot all variables refined in equivalence:\n\t"
s = ""
for v in varlist:
if s != "": s+= " & "
s += str(v)
warnmsg += str(mv) + " => " + s
warnmsg += '\nNot refined: ' + notvaried + '\n'
if zeromult:
errmsg += "\nZero multiplier is invalid in equivalence:\n\t"
s = ""
for v in varlist:
if s != "": s+= " & "
s += str(v)
errmsg += str(mv) + " => " + s + '\n'
# check for errors:
inboth = set(indepVarList).intersection(set(depVarList))
if len(inboth) > 0:
errmsg += "\nThe following parameters(s) are used as both dependent and independent variables in Equivalence relations:\n"
s = ''
for var in sorted(inboth):
if s != "": s+= ", "
s += str(var)
errmsg += '\t'+ s + '\n'
if len(multdepVarList) > 0:
errmsg += "\nThe following parameters(s) are used in multiple Equivalence relations as dependent variables:\n"
s = ''
for var in sorted(set(multdepVarList)):
if s != "": s+= ", "
s += str(var)
errmsg += '\t'+ s + '\n'
equivVarList = list(set(indepVarList).union(set(depVarList)))
if debug: print 'equivVarList',equivVarList
inboth = set(fixVlist).intersection(set(equivVarList))
if len(inboth) > 0:
errmsg += "\nThe following parameter(s) are used in both Equivalence and Fixed constraints:\n"
s = ''
for var in sorted(inboth):
if s != "": s+= ", "
s += str(var)
errmsg += '\t'+ s + '\n'
groups,parmlist = GroupConstraints(constrDict)
# scan through parameters in each relationship. Are all varied? If only some are
# varied, create a warning message.
for group,varlist in zip(groups,parmlist):
if len(varlist) == 1: continue
for rel in group:
varied = 0
notvaried = ''
for var in constrDict[rel]:
if var.startswith('_'): continue
if not re.match('[0-9]*:[0-9\*]*:',var):
warnmsg += "\nVariable "+str(var)+" does not begin with a ':'"
if var in varyList:
varied += 1
else:
if notvaried: notvaried += ', '
notvaried += var
if var in fixVlist:
errmsg += '\nParameter '+var+" is Fixed and used in a constraint:\n\t"
errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
if varied > 0 and varied != len(VarKeys(constrDict[rel])):
warnmsg += "\nNot all variables refined in constraint:\n\t"
warnmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
warnmsg += '\nNot refined: ' + notvaried + '\n'
if errmsg or warnmsg:
return errmsg,warnmsg
# now look for process each group and create the relations that are needed to form
# non-singular square matrix
for group,varlist in zip(groups,parmlist):
if len(varlist) == 1: continue # a constraint group with a single variable can be ignored
if len(varlist) < len(group): # too many relationships -- no can do
errmsg += "\nOver-constrained input. "
errmsg += "There are more constraints " + str(len(group))
errmsg += "\n\tthan variables " + str(len(varlist)) + "\n"
for rel in group:
errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
errmsg += "\n"
continue
try:
multarr = _FillArray(group,constrDict,varlist)
_RowEchelon(len(group),multarr,varlist)
except:
errmsg += "\nSingular input. "
errmsg += "There are internal inconsistencies in these constraints\n"
for rel in group:
errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
errmsg += "\n"
continue
try:
multarr = _FillArray(group,constrDict,varlist,FillDiagonals=True)
GramSchmidtOrtho(multarr,len(group))
except:
errmsg += "\nUnexpected singularity with constraints (in Gram-Schmidt)\n"
for rel in group:
errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
errmsg += "\n"
continue
mapvar = []
group = group[:]
# scan through all generated and input variables
# Check again for inconsistent variable use
# for new variables -- where varied and unvaried parameters get grouped
# together. I don't think this can happen when not flagged before, but
# it does not hurt to check again.
for i in range(len(varlist)):
varied = 0
notvaried = ''
if len(group) > 0:
rel = group.pop(0)
fixedval = fixedList[rel]
for var in VarKeys(constrDict[rel]):
if var in varyList:
varied += 1
else:
if notvaried: notvaried += ', '
notvaried += var
else:
fixedval = None
if fixedval is None:
varname = paramPrefix + str(consNum) # assign a name to a variable
mapvar.append(varname)
consNum += 1
else:
mapvar.append(fixedval)
if varied > 0 and notvaried != '':
warnmsg += "\nNot all variables refined in generated constraint"
warnmsg += '\nPlease report this unexpected error\n'
for rel in group:
warnmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
warnmsg += "\n"
warnmsg += '\n\tNot refined: ' + notvaried + '\n'
try:
np.linalg.inv(multarr)
except:
errmsg += "\nSingular input. "
errmsg += "The following constraints are not "
errmsg += "linearly independent\n\tor do not "
errmsg += "allow for generation of a non-singular set\n"
errmsg += 'Please report this unexpected error\n'
for rel in group:
errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
errmsg += "\n"
return errmsg,warnmsg
[docs]def GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict=None,SeqHist=None):
'''Takes a list of relationship entries comprising a group of
constraints and builds the relationship lists and their inverse
and stores them in global variables Also checks for internal
conflicts or inconsistencies in parameter/variable definitions.
:param list groups: a list of grouped contraints where each constraint
grouped containts a list of indices for constraint constrDict entries,
created in :func:`GroupConstraints` (returned as 1st value)
:param list parmlist: a list containing lists of parameter names
contained in each group, created in :func:`GroupConstraints`
(returned as 2nd value)
:param list varyList: a list of parameters names (strings of form
``<ph>:<hst>:<nam>``) that will be varied. Note that this is changed here.
:param dict constrDict: a list of dicts defining relationships/constraints
(as defined in :func:`GroupConstraints`)
:param list fixedList: a list of values specifying a fixed value for each
dict in constrDict. Values are either strings that can be converted to
floats, float values or None if the constraint defines a new parameter.
:param dict parmDict: a dict containing all parameters defined in current
refinement.
:param int SeqHist: number of current histogram, when used in a sequential
refinement. None (default) otherwise. Wildcard variable names are
set to the current histogram, when found if not None.
'''
global dependentParmList,arrayList,invarrayList,indParmList,consNum
msg = ''
# process fixed (held) variables
for cdict in constrDict:
if len(cdict) == 1:
fixedVarList.append(cdict.keys()[0])
# process equivalences: make a list of dependent and independent vars
# and check for repeated uses (repetition of a parameter as an
# independent var is OK [A=B; A=C], but chaining: [A=B; B=C] is not good)
indepVarList = []
depVarList = []
multdepVarList = []
translateTable = {} # lookup table for wildcard referenced variables
for varlist,mapvars,multarr,invmultarr in zip( # process equivalences
dependentParmList,indParmList,arrayList,invarrayList):
if multarr is None: # true only if an equivalence
zeromult = False
for mv in mapvars:
#s = ''
varied = 0
notvaried = ''
if mv in varyList:
varied += 1
else:
if notvaried: notvaried += ', '
notvaried += mv
if parmDict is not None and mv not in parmDict:
print "Dropping equivalence for variable "+str(mv)+". Not defined in this refinement"
#msg += "\nCannot equivalence to variable "+str(mv)+". Not defined in this refinement"
#continue
else:
if mv not in indepVarList: indepVarList.append(mv)
for v,m in zip(varlist,invmultarr):
if parmDict is not None and v not in parmDict:
print "Dropping equivalence for variable "+str(v)+". Not defined in this refinement"
continue
if m == 0: zeromult = True
if v in varyList:
varied += 1
else:
if notvaried: notvaried += ', '
notvaried += v
if v in depVarList:
multdepVarList.append(v)
else:
depVarList.append(v)
if varied > 0 and varied != len(varlist)+1:
msg += "\nNot all variables refined in equivalence:\n\t"
s = ""
for v in varlist:
if s != "": s+= " & "
s += str(v)
msg += str(mv) + " => " + s
msg += '\nNot refined: ' + notvaried + '\n'
if zeromult:
msg += "\nZero multiplier is invalid in equivalence:\n\t"
s = ""
for v in varlist:
if s != "": s+= " & "
s += str(v)
msg += str(mv) + " => " + s + '\n'
if debug: print 'indepVarList',indepVarList
if debug: print 'depVarList',depVarList
# check for errors:
inboth = set(indepVarList).intersection(set(depVarList))
if len(inboth) > 0:
msg += "\nThe following parameters(s) are used as both dependent and independent variables in Equivalence relations:\n"
s = ''
for var in sorted(inboth):
if s != "": s+= ", "
s += str(var)
msg += '\t'+ s + '\n'
if len(multdepVarList) > 0:
msg += "\nThe following parameters(s) are used in multiple Equivalence relations as dependent variables:\n"
s = ''
for var in sorted(set(multdepVarList)):
if s != "": s+= ", "
s += str(var)
msg += '\t'+ s + '\n'
equivVarList = list(set(indepVarList).union(set(depVarList)))
# scan through parameters in each relationship. Are all varied? If only some are
# varied, create an error message.
for group,varlist in zip(groups,parmlist):
if len(varlist) == 1: continue
for rel in group:
varied = 0
notvaried = ''
unused = 0
notused = ''
for var in constrDict[rel]:
if var.startswith('_'): continue
if var.split(':')[1] == '*' and SeqHist is not None:
# convert wildcard var to reference current histogram; save translation in table
sv = var.split(':')
sv[1] = str(SeqHist)
translateTable[var] = ':'.join(sv)
var = translateTable[var]
if parmDict is not None and var not in parmDict:
unused += 1
if notvaried: notused += ', '
notused += var
if var in varyList:
varied += 1
else:
if notvaried: notvaried += ', '
notvaried += var
if var in fixedVarList:
msg += '\nError: parameter '+var+" is Fixed and used in a constraint:\n\t"
msg += _FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
#if unused > 0:# and unused != len(VarKeys(constrDict[rel])):
if unused > 0 and unused != len(VarKeys(constrDict[rel])):
msg += "\nSome (but not all) variables in constraint are not defined:\n\t"
msg += _FormatConstraint(constrDict[rel],fixedList[rel])
msg += '\nNot used: ' + notused + '\n'
if varied > 0 and varied != len(VarKeys(constrDict[rel])):
msg += "\nNot all variables refined in constraint:\n\t"
msg += _FormatConstraint(constrDict[rel],fixedList[rel])
msg += '\nNot refined: ' + notvaried + '\n'
# if there were errors found, go no farther
if msg:
print ' *** ERROR in constraint definitions! ***'
print msg
raise Exception
# now process each group and create the relations that are needed to form
# a non-singular square matrix
# If all are varied and this is a constraint equation, then set VaryFree flag
# so that the newly created relationships will be varied
for group,varlist in zip(groups,parmlist):
if len(varlist) == 1: continue
# for constraints, if all included variables are refined,
# set the VaryFree flag, and remaining degrees of freedom will be
# varied (since consistency was checked, if any one variable is
# refined, then assume that all are)
varsList = [] # make a list of all the referenced variables as well
VaryFree = False
for rel in group:
varied = 0
unused = 0
for var in VarKeys(constrDict[rel]):
var = translateTable.get(var,var) # replace wildcards
if parmDict is not None and var not in parmDict:
unused += 1
if var not in varsList: varsList.append(var)
if var in varyList: varied += 1
if fixedList[rel] is not None and varied > 0:
VaryFree = True
if len(varlist) < len(group): # too many relationships -- no can do
msg = 'too many relationships'
break
# Since we checked before, if any variables are unused, then all must be.
# If so, this set of relationships can be ignored
if unused:
if debug: print('Constraint ignored (all variables undefined)')
if debug: print (' '+_FormatConstraint(constrDict[rel],fixedList[rel]))
continue
# fill in additional degrees of freedom
try:
arr = _FillArray(group,constrDict,varlist)
_RowEchelon(len(group),arr,varlist)
constrArr = _FillArray(group,constrDict,varlist,FillDiagonals=True)
GramSchmidtOrtho(constrArr,len(group))
except:
msg = 'Singular relationships'
break
mapvar = []
group = group[:]
# scan through all generated and input relationships, we need to add to the varied list
# all the new parameters where VaryFree has been set or where a New Var is varied.
#
# If a group does not contain any fixed values (constraint equations)
# and nothing in the group is varied, drop this group, so that the
# dependent parameters can be refined individually.
unused = True
for i in range(len(varlist)):
if len(group) > 0: # get the original equation reference
rel = group.pop(0)
fixedval = fixedList[rel]
varyflag = constrDict[rel].get('_vary',False)
varname = constrDict[rel].get('_name','')
else: # this relationship has been generated
varyflag = False
varname = ''
fixedval = None
if fixedval is None: # this is a new variable, not a constraint
if not varname:
varname = paramPrefix + str(consNum) # no assigned name, create one
consNum += 1
mapvar.append(varname)
# vary the new relationship if it is a degree of freedom in
# a set of contraint equations or if a New Var is flagged to be varied.
if VaryFree or varyflag:
unused = False
varyList.append(varname)
# fix (prevent varying) of all the variables inside the constraint group
# (dependent vars)
for var in varsList:
if var in varyList: varyList.remove(var)
else:
unused = False
mapvar.append(fixedval)
if unused: # skip over constraints that don't matter (w/o fixed value or any refined variables)
if debug: print('Constraint ignored (all variables unrefined)')
if debug: print (' '+_FormatConstraint(constrDict[rel],fixedList[rel]))
continue
dependentParmList.append([translateTable.get(var,var) for var in varlist])
arrayList.append(constrArr)
invarrayList.append(np.linalg.inv(constrArr))
indParmList.append(mapvar)
if msg:
print ' *** ERROR in constraint definitions! ***'
print msg
print VarRemapShow(varyList)
raise Exception
# setup dictionary containing the fixed values
global fixedDict
# key is original ascii string, value is float
for fixedval in fixedList:
if fixedval:
fixedDict[fixedval] = float(fixedval)
if debug: # on debug, show what is parsed & generated, semi-readable
print 50*'-'
print VarRemapShow(varyList)
print 'Varied: ',varyList
print 'Not Varied: ',fixedVarList
[docs]def StoreEquivalence(independentVar,dependentList):
'''Takes a list of dependent parameter(s) and stores their
relationship to a single independent parameter (independentVar)
:param str independentVar: name of master parameter that will be used to determine the value
to set the dependent variables
:param list dependentList: a list of parameters that will set from
independentVar. Each item in the list can be a string with the parameter
name or a tuple containing a name and multiplier:
``['parm1',('parm2',.5),]``
'''
global dependentParmList,arrayList,invarrayList,indParmList
mapList = []
multlist = []
for var in dependentList:
if isinstance(var, basestring):
mult = 1.0
elif len(var) == 2:
var,mult = var
else:
raise Exception("Cannot parse "+repr(var) + " as var or (var,multiplier)")
mapList.append(var)
multlist.append(tuple((mult,)))
# added relationships to stored values
arrayList.append(None)
invarrayList.append(np.array(multlist))
indParmList.append(tuple((independentVar,)))
dependentParmList.append(mapList)
return
[docs]def GetDependentVars():
'''Return a list of dependent variables: e.g. variables that are
constrained in terms of other variables
:returns: a list of variable names
'''
dependentVars = []
global dependentParmList
for lst in dependentParmList:
for itm in lst: dependentVars.append(itm)
return dependentVars
[docs]def GetIndependentVars():
'''Return a list of independent variables: e.g. variables that are
created by constraints of other variables
:returns: a list of variable names
'''
independentVars = []
global indParmList,fixedDict
for lst in indParmList:
for name in lst:
if name in fixedDict: continue
independentVars.append(name)
return independentVars
[docs]def PrintIndependentVars(parmDict,varyList,sigDict,PrintAll=False,pFile=None):
'''Print the values and uncertainties on the independent variables'''
global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
printlist = []
mapvars = GetIndependentVars()
for i,name in enumerate(mapvars):
if name in fixedDict: continue
if PrintAll or name in varyList:
sig = sigDict.get(name)
printlist.append([name,parmDict[name],sig])
if len(printlist) == 0: return
s1 = ''
print >>pFile,130*'-'
print >>pFile,"Variables generated by constraints"
printlist.append(3*[None])
for name,val,esd in printlist:
if len(s1) > 120 or name is None:
print >>pFile,''
print >>pFile,s1
print >>pFile,s2
print >>pFile,s3
s1 = ''
if name is None: break
if s1 == "":
s1 = ' name :'
s2 = ' value :'
s3 = ' sig :'
s1 += '%15s' % (name)
s2 += '%15.4f' % (val)
if esd is None:
s3 += '%15s' % ('n/a')
else:
s3 += '%15.4f' % (esd)
[docs]def ComputeDepESD(covMatrix,varyList,parmDict):
'''Compute uncertainties for dependent parameters from independent ones
returns a dictionary containing the esd values for dependent parameters
'''
sigmaDict = {}
for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
#if invmultarr is None: continue # probably not needed
try:
valuelist = [parmDict[var] for var in mapvars]
except KeyError:
continue
# get the v-covar matrix for independent parameters
vcov = np.zeros((len(mapvars),len(mapvars)))
for i1,name1 in enumerate(mapvars):
if name1 not in varyList: continue
iv1 = varyList.index(name1)
for i2,name2 in enumerate(mapvars):
if name2 not in varyList: continue
iv2 = varyList.index(name2)
vcov[i1][i2] = covMatrix[iv1][iv2]
# vec is the vector that multiplies each of the independent values
for v,vec in zip(varlist,invmultarr):
sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec)))
return sigmaDict
def _FormatConstraint(RelDict,RelVal):
'''Formats a Constraint or Function for use in a convenient way'''
linelen = 45
s = [""]
for var,val in RelDict.items():
if var.startswith('_'): continue
if len(s[-1]) > linelen: s.append(' ')
m = val
if s[-1] != "" and m >= 0:
s[-1] += ' + '
elif s[-1] != "":
s[-1] += ' - '
m = abs(m)
s[-1] += '%.3f*%s '%(m,var)
if len(s[-1]) > linelen: s.append(' ')
if RelVal is None:
s[-1] += ' = New variable'
else:
s[-1] += ' = ' + RelVal
s1 = ''
for s2 in s:
if s1 != '': s1 += '\n\t'
s1 += s2
return s1
[docs]def VarRemapShow(varyList):
'''List out the saved relationships. This should be done after the constraints have been
defined using :func:`StoreEquivalence`, :func:`GroupConstraints` and :func:`GenerateConstraints`.
:returns: a string containing the details of the contraint relationships
'''
s = ''
if len(fixedVarList) > 0:
s += 'Fixed Variables:\n'
for v in fixedVarList:
s += ' ' + v + '\n'
s += 'Variable mapping relations:\n'
global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
for varlist,mapvars,multarr,invmultarr in zip(
dependentParmList,indParmList,arrayList,invarrayList):
for i,mv in enumerate(mapvars):
if multarr is None:
s += ' ' + str(mv) + ' is equivalent to parameter(s): '
j = 0
for v,m in zip(varlist,invmultarr):
if debug: print 'v,m[0]: ',v,m[0]
if j > 0: s += ' & '
j += 1
s += str(v)
if m != 1:
s += " / " + str(m[0])
s += '\n'
continue
s += ' %s = ' % mv
j = 0
for m,v in zip(multarr[i,:],varlist):
if m == 0: continue
if j > 0: s += ' + '
j += 1
s += '(%s * %s)' % (m,v)
if mv in varyList: s += ' VARY'
s += '\n'
s += 'Inverse variable mapping relations:\n'
for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
for i,mv in enumerate(varlist):
s += ' %s = ' % mv
j = 0
for m,v in zip(invmultarr[i,:],mapvars):
if m == 0: continue
if j > 0: s += ' + '
j += 1
s += '(%s * %s)' % (m,v)
s += '\n'
return s
[docs]def Dict2Deriv(varyList,derivDict,dMdv):
'''Compute derivatives for Independent Parameters from the
derivatives for the original parameters
:param list varyList: a list of parameters names that will be varied
:param dict derivDict: a dict containing derivatives for parameter values keyed by the
parameter names.
:param list dMdv: a Jacobian, as a list of np.array containing derivatives for dependent
parameter computed from derivDict
'''
global dependentParmList,arrayList,invarrayList,indParmList,invarrayList
for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
for i,name in enumerate(mapvars):
# grouped variables: need to add in the derv. w/r
# dependent variables to the independent ones
if name not in varyList: continue # skip if independent var not varied
if multarr is None:
for v,m in zip(varlist,invmultarr):
if debug: print 'start dMdv',dMdv[varyList.index(name)]
if debug: print 'add derv',v,'/',m[0],'to derv',name,'add=',derivDict[v] / m[0]
if m == 0: continue
dMdv[varyList.index(name)] += derivDict[v] / m[0]
else:
for v,m in zip(varlist,multarr[i,:]):
if debug: print 'start dMdv',dMdv[varyList.index(name)]
if debug: print 'add derv',v,'*',m,'to derv',name,'add=',m * derivDict[v]
if m == 0: continue
dMdv[varyList.index(name)] += m * derivDict[v]
[docs]def Map2Dict(parmDict,varyList):
'''Create (or update) the Independent Parameters from the original
set of Parameters
Removes dependent variables from the varyList
This should be done once, after the constraints have been
defined using :func:`StoreEquivalence`,
:func:`GroupConstraints` and :func:`GenerateConstraints` and
before any variable refinement is done
to complete the parameter dictionary by defining independent
parameters and satisfying the constraint equations.
:param dict parmDict: a dict containing parameter values keyed by the
parameter names.
This will contain updated values for both dependent and independent
parameters after Dict2Map is called. It will also contain some
unexpected entries of every constant value {'0':0.0} & {'1.0':1.0},
which do not cause any problems.
:param list varyList: a list of parameters names that will be varied
'''
# process the Independent vars: remove dependent ones from varylist
# and then compute values for the Independent ones from their dependents
global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
for item in varlist:
try:
varyList.remove(item)
except ValueError:
pass
if multarr is None: continue
valuelist = [parmDict[var] for var in varlist]
parmDict.update(zip(mapvars,
np.dot(multarr,np.array(valuelist)))
)
# now remove fixed variables from the varyList
global fixedVarList
for item in fixedVarList:
try:
varyList.remove(item)
except ValueError:
pass
# Set constrained parameters to their fixed values
parmDict.update(fixedDict)
[docs]def Dict2Map(parmDict,varyList):
'''Applies the constraints defined using :func:`StoreEquivalence`,
:func:`GroupConstraints` and :func:`GenerateConstraints` by changing
values in a dict containing the parameters. This should be
done before the parameters are used for any computations
:param dict parmDict: a dict containing parameter values keyed by the
parameter names.
This will contain updated values for both dependent and independent
parameters after Dict2Map is called. It will also contain some
unexpected entries of every constant value {'0':0.0} & {'1.0':1.0},
which do not cause any problems.
:param list varyList: a list of parameters names that will be varied
'''
global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
# reset fixed values (should not be needed, but very quick)
# - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended
# not needed, but also not a problem - BHT
parmDict.update(fixedDict)
for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
#if invmultarr is None: continue
try:
valuelist = [parmDict[var] for var in mapvars]
except KeyError:
continue
parmDict.update(zip(varlist,
np.dot(invmultarr,np.array(valuelist)))
)
#======================================================================
# internal routines follow (these routines are unlikely to be called
# from outside the module)
[docs]def GramSchmidtOrtho(a,nkeep=0):
'''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to
find orthonormal unit vectors relative to first row.
If nkeep is non-zero, the first nkeep rows in the array are not changed
input:
arrayin: a 2-D non-singular square array
returns:
a orthonormal set of unit vectors as a square array
'''
def proj(a,b):
'Projection operator'
return a*(np.dot(a,b)/np.dot(a,a))
for j in range(nkeep,len(a)):
for i in range(j):
a[j] -= proj(a[i],a[j])
if np.allclose(np.linalg.norm(a[j]),0.0):
raise Exception,"Singular input to GramSchmidtOrtho"
a[j] /= np.linalg.norm(a[j])
return a
def _FillArray(sel,dict,collist,FillDiagonals=False):
'''Construct a n by n matrix (n = len(collist)
filling in the rows using the relationships defined
in the dictionaries selected by sel
If FillDiagonals is True, diagonal elements in the
array are set to 1.0
'''
n = len(collist)
if FillDiagonals:
arr = np.eye(n)
else:
arr = np.zeros(2*[n,])
# fill the top rows
for i,cnum in enumerate(sel):
for j,var in enumerate(collist):
arr[i,j] = dict[cnum].get(var,0)
return arr
def _SwapColumns(i,m,v):
'''Swap columns in matrix m as well as the labels in v
so that element (i,i) is replaced by the first non-zero element in row i after that element
Throws an exception if there are no non-zero elements in that row
'''
for j in range(i+1,len(v)):
if not np.allclose(m[i,j],0):
m[:,(i,j)] = m[:,(j,i)]
v[i],v[j] = v[j],v[i]
return
else:
raise Exception,'Singular input'
def _RowEchelon(m,arr,collist):
'''Convert the first m rows in Matrix arr to row-echelon form
exchanging columns in the matrix and collist as needed.
throws an exception if the matrix is singular because
the first m rows are not linearly independent
'''
n = len(collist)
for i in range(m):
if np.allclose(arr[i,i],0):
_SwapColumns(i,arr,collist)
arr[i,:] /= arr[i,i] # normalize row
# subtract current row from subsequent rows to set values to left of diagonal to 0
for j in range(i+1,m):
arr[j,:] -= arr[i,:] * arr[j,i]
if __name__ == "__main__":
parmdict = {}
constrDict = [
{'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},
{'0:0:eA': 0.0},
{'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
{'1::C(10,0,1)': 1.0, '2::C(10,0,1)': 1.0},
{'1::AUiso:0': 1.0, '0::AUiso:0': 1.0},
{'0::A0': 0.0}
]
fixedList = ['5.0', '0', None, None, '1.0', '0']
StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), ))
#StoreEquivalence('1::atomx:3',('2::atomx:3', ('2::atomz:3',2,), )) # error: dependent & independent vars mixed
#StoreEquivalence('1::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) # error: dependent vars repeated
#StoreEquivalence('0:1:eA',('0:0:eA',)) # error: equiv & fixed
#StoreEquivalence('0:99:Scale',('0:12:Scale',)) # error: equiv & constrained
#StoreEquivalence('0:12:Scale',('0:99:Scale',)) # error: equiv & constrained
varylist = ['2::atomx:3',
'2::C(10,6,1)', '1::C(10,6,1)',
'2::atomy:3', '2::atomz:3',
'0:12:Scale', '0:11:Scale', '0:14:Scale', '0:13:Scale', '0:0:Scale']
# e,w = CheckConstraints([,
# [{'2:0:Scale': 1.0, '5:0:Scale': 1.0, '10:0:Scale': 1.0, '6:0:Scale': 1.0, '9:0:Scale': 1.0, '8:0:Scale': 1.0,# '3:0:Scale': 1.0, '4:0:Scale': 1.0, '7:0:Scale': 1.0, '1:0:Scale': 1.0, '0:0:Scale': 1.0}],
# ['1.0'])
# if e: print 'error=',e
# if w: print 'error=',w
# 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']
# constrDict = [
# {'0::Afrac:4': 24.0, '0::Afrac:1': 16.0, '0::Afrac:3': 24.0, '0::Afrac:2': 16.0},
# {'0::Afrac:1': 1.0, '0::Afrac:2': 1.0},
# {'0::Afrac:4': 1.0, '0::Afrac:3': 1.0}]
# fixedList = ['40.0', '1.0', '1.0']
errmsg, warnmsg = CheckConstraints(varylist,constrDict,fixedList)
if errmsg:
print "*** Error ********************"
print errmsg
if warnmsg:
print "*** Warning ********************"
print warnmsg
if errmsg or warnmsg:
sys.exit()
groups,parmlist = GroupConstraints(constrDict)
GenerateConstraints(groups,parmlist,varylist,constrDict,fixedList)
print VarRemapShow(varylist)
parmdict.update( {
'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,
'0:0:eA': 0.0,
'2::C(10,6,1)': 0.2, '1::C(10,6,1)': 0.3,
'1::C(10,0,1)': 0.2, '2::C(10,0,1)': 0.3,
'1::AUiso:0': 0.02, '0::AUiso:0': 0.03,
'0::A0': 0.0,
'2::atomx:3':0.23,'2::atomy:3':-.23, '2::atomz:3':-0.11,
})
print 'parmdict start',parmdict
print 'varylist start',varylist
before = parmdict.copy()
Map2Dict(parmdict,varylist)
print 'parmdict before and after Map2Dict'
print ' key / before / after'
for key in sorted(parmdict.keys()):
print ' '+key,'\t',before.get(key),'\t',parmdict[key]
print 'varylist after',varylist
before = parmdict.copy()
Dict2Map(parmdict,varylist)
print 'after Dict2Map'
print ' key / before / after'
for key in sorted(parmdict.keys()):
print ' '+key,'\t',before.get(key),'\t',parmdict[key]
# dMdv = len(varylist)*[0]
# deriv = {}
# for i,v in enumerate(parmdict.keys()): deriv[v]=i
# Dict2Deriv(varylist,deriv,dMdv)
