# TODO: revisit SeqRefine and :meth:`GSASIIdataGUI.GSASII.OnSeqRefine` and :func:`GSASIIseqGUI.UpdateSeqResults`
# -*- coding: utf-8 -*-
########### SVN repository information ###################
# $Date: 2024-02-07 14:07:07 -0600 (Wed, 07 Feb 2024) $
# $Author: toby $
# $Revision: 5723 $
# $URL: https://subversion.xray.aps.anl.gov/pyGSAS/trunk/GSASIImapvars.py $
# $Id: GSASIImapvars.py 5723 2024-02-07 20:07:07Z toby $
########### SVN repository information ###################
"""
Classes and routines defined in :mod:`GSASIImapvars` follow.
Note that parameter names in GSAS-II are strings of form ``<ph#>:<hst#>:<nam>`` or ``<ph#>::<nam>:<at#>``
where ``<ph#>`` is a phase number, ``<hst#>`` is a histogram number and ``<at#>`` is an atom number.
``<nam>`` is a name that determines the parameter type (see :func:`GSASIIobj.CompileVarDesc`). When
stored in the data tree, parameters are saved as :class:`GSASIIobj.G2VarObj` objects
so that they can be resolved if the phase/histogram order changes.
"""
# Note that documentation for GSASIImapvars.py has been moved
# to file docs/source/GSASIImapvars.rst
from __future__ import division, print_function
import copy
import numpy as np
import GSASIIpath
GSASIIpath.SetVersionNumber("$Revision: 5723 $")
import GSASIIobj as G2obj
# data used for constraints;
debug = False # turns on printing as constraint input is processed
#------------------------------------------------------------------------------------
# Global vars used for storing constraints and equivalences after processing
# note that new var and constraint equations are stored together in groups,
# where each constraint group contains those parameters that must be handled
# together. Equivalences are also stored in these
dependentParmList = []
'''a list of lists where each item contains a list of parameters in each constraint group.
note that parameters listed in dependentParmList should not be refined directly.'''
indParmList = [] # a list of names for the new parameters
'''a list of lists where each item contains a list for each constraint group with
fixed values for constraint equations and names of generated/New Var parameters.
In the case of equivalences, the name of a single independent parameter is stored.
'''
arrayList = []
'''a list of of relationship matrices that map model parameters in each
constraint group (in :data:`dependentParmList`) to
generated (New Var) parameters.
'''
invarrayList = []
'''a list of of inverse-relationship matrices that map constrained values and
generated (New Var) parameters (in :data:`indParmList`) to model parameters
(in :data:`dependentParmList`).
'''
symGenList = []
'''A list of flags that if True indicates a constraint was generated by symmetry
'''
holdParmList = []
'''List of parameters that should not be refined ("Hold"s).
Set in :func:`StoreHold`. Initialized in :func:`InitVars`.
'''
holdParmType = {}
'''The reason why a parameter has been marked as "Hold".
Initialized in :func:`InitVars`; set in :func:`StoreHold`.
'''
constrParms = {'indep-equiv':[], 'indep-constr':[], 'dep-equiv':[], 'dep-constr':[]}
'''A dict with parameters in equivalences, compiled from
(:data:`dependentParmList`) and (:data:`indParmList`).
Used within :func:`GetIndependentVars` and :func:`GetDependentVars`.
'''
saveVaryList = []
'''A list of the varied parameters that was last supplied when constraints were
processed. This is set in :func:`GenerateConstraints` and updated in
:func:`Map2Dict`. Used in :func:`VarRemapShow`
'''
droppedSym = []
'''A list of symmetry generated equivalences that have been converted to
constraint equations in :func:`CheckEquivalences`
'''
#------------------------------------------------------------------------------------
# Global vars set in :func:`GenerateConstraints`. Used for intermediate processing of
# constraints.
constrVarList = []
'List of parameters used in "Constr" and "New Var" constraints'
indepVarList = []
'A list of all independent parameters in equivalences'
depVarList = []
'A list of all dependent parameters in equivalences'
#------------------------------------------------------------------------------------
# global variables used in :func:`getConstrError` to store error and warning information.
# set in CheckEquivalences and in GenerateConstraints
convVarList = []
'parameters in equivalences that were converted to "Const" constraints'
multdepVarList = []
'parameters used as dependents multiple times in equivalences'
unvariedParmsList = []
'parameters used in equivalences that are not varied'
undefinedVars = []
'parameters used in equivalences that are not defined in the parameter dict'
groupErrors = []
'parameters in constraints where parameter grouping and matrix inversion fails'
#------------------------------------------------------------------------------------
paramPrefix = "::constr"
'A prefix for generated parameter names'
consNum = 0
'The number to be assigned to the next constraint to be created'
#------------------------------------------------------------------------------------
[docs]
class ConstraintException(Exception):
'''Defines an Exception that is used when an exception is raised processing constraints.
Raised in :func:`GenerateConstraints` during sequential fits. Possible (but highly unlikely)
to be raised in :func:`CheckEquivalences` (called by :func:`GenerateConstraints`) if an
infinite loop is detected.
Also raised in :func:`GramSchmidtOrtho` and :func:`_SwapColumns` but caught
within :func:`GenerateConstraints`.
'''
pass
[docs]
def InitVars():
'''Initializes all constraint information'''
global dependentParmList,arrayList,invarrayList,indParmList,consNum,symGenList
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
consNum = 0 # number of the next constraint to be created
symGenList = [] # Flag if constraint is generated by symmetry
global holdParmList,holdParmType
holdParmList = []
holdParmType = {}
global droppedSym
droppedSym = []
global constrVarList,indepVarList,depVarList
constrVarList = []
indepVarList = []
depVarList = []
[docs]
def VarKeys(constr):
"""Finds the keys in a constraint that represent parameters
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,constrI 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(constrI))
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,constrJ in enumerate(constrDict):
if j in assignedlist: continue # already in a group, skip
if len(set(VarKeys(constrJ)) & groupset) > 0: # true if this needs to be added
changes = True
grouplist.append(j)
assignedlist.append(j)
groupset = groupset | set(VarKeys(constrJ))
group = sorted(grouplist)
varlist = sorted(list(groupset))
groups.append(group)
ParmList.append(varlist)
return groups,ParmList
[docs]
def GenerateConstraints(varyList,constrDict,fixedList,parmDict=None,
seqHistNum=None,raiseException=False):
'''Takes a list of relationship entries that have been stored by
:func:`ProcessConstraints` into lists ``constrDict`` and ``fixedList``
This routine then calls :func:`CheckEquivalences` for internal
consistency. This includes converting equivalenced variables into
constraints when a variable is used in both.
Once checked, parameters are grouped so that any parameters that are used in
more than one constraint are grouped together. This allows checking for incompatible
logic (for example, when four constraints are specified for three variables).
If parmDict is not None, the parameter groups are checked for constraints where
some parameters are varied, but not others. If so, the value for that unvaried
parameter is subtracted from the constant in the constraint.
Once all checks are complete, the constraints are then
converted to the form used to apply them, saving them as global variables within
this module.
:param list varyList: a list of parameters names (strings of form
``<ph>:<hst>:<nam>``) that will be varied. Note that this is changed
here unless set to None. None is used to indicate that all constraints
should be generated.
:param dict constrDict: a list of dicts defining relationships/constraints
(as described 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 seqHistNum: the hId number of the current histogram in a sequential
fit. None (default) otherwise.
:param bool raiseException: When True, generation of an error causes
an exception to be raised (used in sequential fits)
:returns: errmsg,warning,groups,parmlist
**errmsg**
Is an error message or empty if no errors were found
**warning**
Is a warning message about constraints that have been ignored or changed
**groups**
Lists parameter groups
**parmlist**
Lists parameters in each parameter groups
'''
warninfo = {'msg':'', 'shown':{}}
def warn(msg,cdict=None,val=None):
if cdict is not None and cdict != warninfo['shown']:
warninfo['shown'] = cdict
if warninfo['msg']: warninfo['msg'] += '\n'
if '_vary' in cdict:
warninfo['msg'] += '\nProblem with new var expression: ' + _FormatConstraint(cdict,cdict.get('_name','New Var'))
else:
warninfo['msg'] += '\nProblem with constraint equation: ' + _FormatConstraint(cdict,val)
if warninfo['msg']: warninfo['msg'] += '\n'
warninfo['msg'] += ' ' + msg
global dependentParmList,arrayList,invarrayList,indParmList,consNum
# lists of parameters used for error reporting
global undefinedVars # parameters that are used in equivalences but are not defined
undefinedVars = []
global groupErrors # parameters in constraints that cause grouping errors
groupErrors = []
global saveVaryList
saveVaryList = copy.copy(varyList)
errmsg = '' # save error messages here. If non-blank, constraints cannot be used.
warning = '' # save informational text messages here.
# find parameters used in constraint equations & new var assignments (all are dependent)
global constrVarList
constrVarList = []
for cnum,(cdict,fixVal) in enumerate(zip(constrDict,fixedList)):
constrVarList += [i for i in cdict if i not in constrVarList and not i.startswith('_')]
# Process the equivalences; If there are conflicting parameters, move them into constraints
warning = CheckEquivalences(constrDict,varyList,fixedList,parmDict,seqHistNum=seqHistNum)
# look through "Constr" and "New Var" constraints looking for zero multipliers and
# Hold, Unvaried & Undefined parameters
skipList = []
invalidParms = []
for cnum,(cdict,fixVal) in enumerate(zip(constrDict,fixedList)):
#constrVarList += [i for i in cdict if i not in constrVarList and not i.startswith('_')]
valid = 0 # count of good parameters
# error reporting
zeroList = [] # parameters with zero multipliers
holdList = [] # parameters with "Hold"'s
noVaryList = [] # parameters not varied
noWildcardList = [] # wildcard parameters in non-sequential fit
notDefList = [] # parameters not defined
# processing to be done
problem = False # constraint must be dropped
dropList = [] # parameters to remove
for var in VarKeys(cdict): # assemble warning info
if cdict[var] == 0: # zero multiplier
if var not in zeroList: zeroList.append(var)
if var not in dropList: dropList.append(var)
elif var in holdParmList: # hold invalid in New Var, drop from constraint eqn
holdList.append(var)
if fixVal is None: # hold in a newvar is not allowed
problem = True
else:
if var not in dropList: dropList.append(var)
elif ':*:' in var : # wildcard still present should be treated as undefined
if var not in undefinedVars: undefinedVars.append(var)
noWildcardList.append(var)
problem = True
elif parmDict is not None and var not in parmDict: # not defined, constraint will not be used
if var not in undefinedVars: undefinedVars.append(var)
notDefList.append(var)
if seqHistNum is None:
if ':dAx:' in var or ':dAy:' in var or ':dAz:' in var: # coordinates from undefined atoms
if fixVal is None:
problem = True # invalid in New Var
else:
if var not in dropList: dropList.append(var) # ignore in constraint eqn
else:
problem = True
elif varyList is None:
valid += 1
elif var not in varyList and fixVal is not None: # unvaried, constraint eq. only
if var not in unvariedParmsList: unvariedParmsList.append(var)
noVaryList.append(var)
dropList.append(var)
else:
valid += 1
if seqHistNum is not None and len(notDefList) > 0 and valid == 0:
# for sequential ref can quietly ignore constraints with all undefined vars
notDefList = []
for l,m in ((zeroList,"have zero multipliers"), # show warning
(holdList,'set as "Hold"'),
(noVaryList,"not varied"),
(noWildcardList,"wildcard in non-sequential fit"),
(notDefList,"not defined")):
if l and cdict.get('_vary',True): # true for constr eq & varied New Var
msg = "parameter(s) " + m + ': '
for i,v in enumerate(l):
if i != 0: msg += ', '
msg += v
warn(msg,cdict,fixVal)
if valid == 0: # no valid entries
if seqHistNum is None: warn('Ignoring constraint, no valid entries',cdict)
skipList.append(cnum)
elif problem: # mix of valid & refined and undefined items, cannot use this
if cdict.get('_vary',True): # true for constr eq & varied New Var
warn('New Var constraint will be ignored',cdict)
skipList.append(cnum)
invalidParms += VarKeys(cdict)
elif len(dropList) > 0: # mix of valid and problematic items, drop problem vars, but keep rest
if GSASIIpath.GetConfigValue('debug'):
msg = ''
for v in dropList:
if msg: msg += ' ,'
msg += v
warn('removing: '+msg,cdict)
value = fixedList[cnum]
for var in dropList: # do cleanup
# NB expressions in constraint multipliers have already been evaluated
if ':dAx:' in var or ':dAy:' in var or ':dAz:' in var:
pass # treat delta coords as 0; no change in total is needed
elif cdict[var] != 0:
value = float(value) - cdict[var]*parmDict[var]
del cdict[var]
if float(value) != float(fixedList[cnum]): fixedList[cnum] = float(np.round(value,12))
if GSASIIpath.GetConfigValue('debug'):
warn('revised as: '+_FormatConstraint(constrDict[cnum],fixedList[cnum]))
for i in list(range(len(constrDict)-1,-1,-1)): # remove the dropped constraints
if i in skipList:
del constrDict[i]
del fixedList[i]
for i in invalidParms: StoreHold(i,"Used in invalid constraint")
if warning: warning += '\n'
warning += warninfo['msg']
groups,parmlist = GroupConstraints(constrDict)
# now process each group and create the relations that are needed to form
# a non-singular square matrix
# Now check that all parameters are varied (probably do not need to do this
# any more). For constraint equations, if all are varied, set VaryFree to True
# and all newly created relationships will be varied. For NewVar constraints,
# vary if the vary flag was set.
for group,depPrmList in zip(groups,parmlist):
if len(depPrmList) < len(group): # too many relationships -- no can do
if errmsg: errmsg += '\n'
errmsg += "Over-constrained input. "
errmsg += "There are more constraints (" + str(len(group))
errmsg += ") than parameters (" + str(len(depPrmList)) + ")\nin these constraints:"
for rel in group:
errmsg += '\n\t'+ _FormatConstraint(constrDict[rel],fixedList[rel])
groupErrors += depPrmList
continue # go on to next group
try:
constrArr = _FillArray(group,constrDict,depPrmList)
except Exception as err:
if errmsg: errmsg += '\n'
if 'Initial' in str(err):
errmsg += "\nSingular input. "
errmsg += "There are internal inconsistencies in these constraints:"
else:
errmsg += "\nError expanding matrix with these constraints:"
for rel in group:
errmsg += '\n\t' + _FormatConstraint(constrDict[rel],fixedList[rel])
groupErrors += depPrmList
continue
try:
GramSchmidtOrtho(constrArr,len(group))
except:
if errmsg: errmsg += '\n'
errmsg += "\nUnexpected singularity with constraints group (in Gram-Schmidt)"
for rel in group:
errmsg += '\n\t' + _FormatConstraint(constrDict[rel],fixedList[rel])
groupErrors += depPrmList
continue
try:
invConstrArr = np.linalg.inv(constrArr)
except:
if errmsg: errmsg += '\n'
errmsg += "\nSingular input. "
errmsg += "The following constraints are not "
errmsg += "linearly independent\nor do not "
errmsg += "allow for generation of a non-singular set.\n"
errmsg += 'This is unexpected. Please report this (toby@anl.gov)'
for rel in group:
errmsg += '\n\t' + _FormatConstraint(constrDict[rel],fixedList[rel])
groupErrors += depPrmList
continue
# scan through current group looking for new var assignments
hasNewVar = False
for rel in group:
if fixedList[rel] is None:
hasNewVar = True # there a New Var relationship in this group
break
else: # only constraint equations, check for unvaried parameters in one
unvaried = False
# this should not happen as they should have been removed
for var in depPrmList:
if varyList is not None and var not in varyList:
unvaried = True
break
if unvaried: # something is not varied: skip group & remove all parameters from varyList
for var in depPrmList: StoreHold(var,'Unexpected: mixed use')
if GSASIIpath.GetConfigValue('debug'):
print('Unexpected: Constraint group ignored (some parameters unvaried)')
for rel in group:
print (' '+_FormatConstraint(constrDict[rel],fixedList[rel]))
continue
maplist = [] # value or param name mapped to each row in expression matrix
if not hasNewVar: # constraint equations; all remaining entries varied, vary generated
for i in range(len(depPrmList)):
if i >= len(group): # tag generated degrees of freedom with names and vary them
varname = paramPrefix + str(consNum) # assign a unique name
consNum += 1
maplist.append(varname)
if varyList is not None: varyList.append(varname)
else:
rel = group[i]
maplist.append(fixedList[rel])
else: # ------------------------- groups with new var assignments, vary only NV's w/flags set
for i,rel in enumerate(group):
if fixedList[rel] is None:
varname = constrDict[rel].get('_name','::?????')
maplist.append(varname)
if constrDict[rel].get('_vary',False):
if varyList is not None: varyList.append(varname)
else:
maplist.append(fixedList[rel]) # constraint equation
for i in range(len(depPrmList)):
if i >= len(group): # name generated degrees of freedom
varname = paramPrefix + str(consNum) # assign a unique name
consNum += 1
maplist.append(varname)
for var in depPrmList: StoreHold(var,'New Var use')
# keep this group
dependentParmList.append(depPrmList)
arrayList.append(constrArr)
invarrayList.append(invConstrArr)
indParmList.append(maplist)
symGenList.append(False)
if errmsg and raiseException:
print (' *** ERROR in constraint definitions! ***')
print (errmsg)
if warning:
print (' also note warnings in constraint processing:')
print (warning)
raise ConstraintException
elif errmsg:
return errmsg,warning,None,None
# Make list of dependent and independent variables for all constraints
constrParms['dep-equiv'] = []
constrParms['dep-constr'] = []
constrParms['indep-equiv'] = []
constrParms['indep-constr'] = []
for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList): # process all constraints
for mv in mapvars:
if type(mv) is float or type(mv) is int: continue
if multarr is None and mv not in constrParms['indep-equiv']:
constrParms['indep-equiv'].append(mv)
elif mv not in constrParms['indep-constr']:
constrParms['indep-constr'].append(mv)
for mv in varlist:
if multarr is None and mv not in constrParms['dep-equiv']:
constrParms['dep-equiv'].append(mv)
elif mv not in constrParms['dep-constr']:
constrParms['dep-constr'].append(mv)
StoreHold(mv,'dependent param')
saveVaryList = copy.copy(varyList) # save varyList so it can be used within module
if varyList is None: saveVaryList = []
# if equivMoved:
# print(60*'=')
# print('Constraints were reclassified to avoid conflicts, as below:')
# print(mvMsg)
# print('New constraints are:')
# print (VarRemapShow(varyList,True))
# print(60*'=')
return errmsg,warning,groups,parmlist # saved for sequential fits
[docs]
def CheckEquivalences(constrDict,varyList,fixedList,parmDict=None,seqHistNum=None):
'''Process equivalence constraints, looking for conflicts such as
where a parameter is used in both an equivalence and a constraint expression
or where chaining is done (A->B and B->C).
Removes equivalences or parameters from equivalences or converts equivalences to
constraints as described for :ref:`Equivalence Checking and Reorganization <CheckEquivalences>`.
:param dict constrDict: a list of dicts defining relationships/constraints
:param list varyList: list of varied parameters (defined during refinements only)
: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.
:param dict parmDict: a dict containing defined parameters and their values. Used to find
equivalences where a parameter is has been removed from a refinement.
:param int seqHistNum: the hId number of the current histogram in a sequential
fit. None (default) otherwise.
:returns: warning messages about changes that need to be made to equivalences
'''
warninfo = {'msg':'', 'shown':-1}
def warnEqv(msg,cnum=None):
if cnum is not None and cnum != warninfo['shown']:
warninfo['shown'] = cnum
if warninfo['msg']: warninfo['msg'] += '\n'
warninfo['msg'] += '\nProblem with equivalence: ' + _showEquiv(
dependentParmList[cnum],indParmList[cnum],invarrayList[cnum])
if warninfo['msg']: warninfo['msg'] += '\n'
warninfo['msg'] += ' ' + msg
global depVarList # parameters used in equivalences as dependent parameters
global indepVarList # parameters used in equivalences as independent parameters
global constrVarList # parameters used in other constraints
# lists of parameters used for error reporting
global undefinedVars # parameters that are used in equivalences but are not defined
global convVarList # parameters in equivalences that will be converted to constraints
convVarList = [] # parameters in equivalences to be made into constraints
global multdepVarList
multdepVarList = [] # list of dependent parameters used in more than one equivalence
# local vars
dropVarList = [] # parameters that can be removed from equivalences
removeList = [] # equivalences that are not needed
convertList = [] # equivalences that should be converted to "Const" constraints
# tabulate parameters used in equivalences by type
depVarList = [] # list of all dependent parameters in equivalences
indepVarList = [] # list of all independent parameters in equivalences
for cnum,(varlist,mapvars,multarr,invmultarr) in enumerate(zip(
dependentParmList,indParmList,arrayList,invarrayList)):
#if multarr is not None: continue # equivalence
indepVarList += [mv for mv in mapvars if mv not in indepVarList]
depVarList += [v for v in varlist if v not in depVarList]
# 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)
# look for parameters in equivalences that are used more than once as dependent parameters
seenOnce = []
for cnum,(varlist,multarr) in enumerate(zip(dependentParmList,arrayList)):
if multarr is not None: continue # equivalences only
for v in varlist:
if v not in seenOnce:
seenOnce.append(v)
elif v not in multdepVarList:
multdepVarList.append(v)
# scan through equivalences looking for other "dual uses". Stop when no new ones are found
changed = True
count = 0
while changed:
changed = False
count += 1
if count > 1000:
raise ConstraintException("Too many loops in CheckEquivalences")
# look for repeated dependent vars
convVarList = [] # parameters in equivalences to be made into constraints
for cnum,(varlist,mapvars,multarr,invmultarr) in enumerate(zip(
dependentParmList,indParmList,arrayList,invarrayList)):
if multarr is not None: continue # equivalences only
if cnum in convertList:
convVarList += [v for v in mapvars+varlist if v not in convVarList and type(v) is not float]
continue
# identify equivalences that need to be converted to constraints.
# Where parameters:
# are used both in equivalences & constraints,
# are used as dependent multiple times or
# where are used as both dependent and independent (chained)
msg = False
for v in mapvars:
if v in constrVarList+convVarList:
changed = True
msg = True
warnEqv("Independent parameter "+str(v)+' used in constraint',cnum)
if cnum not in convertList: convertList.append(cnum)
for v in varlist:
if v in multdepVarList:
changed = True
msg = True
warnEqv("Dependent parameter "+str(v)+' repeated',cnum)
if cnum not in convertList: convertList.append(cnum)
elif v in indepVarList:
changed = True
msg = True
warnEqv("Dependent parameter "+str(v)+' used elsewhere as independent',cnum)
if cnum not in convertList: convertList.append(cnum)
elif v in constrVarList+convVarList:
changed = True
msg = True
warnEqv("Dependent parameter "+str(v)+' used in constraint',cnum)
if cnum not in convertList: convertList.append(cnum)
if msg:
warnEqv('Converting to "Constr"',cnum)
if warninfo['msg'] and seqHistNum is not None:
# sequential fit -- print the recasts but don't stop with a warning
print(warninfo['msg'])
warninfo = {'msg':'', 'shown':-1}
global unvariedParmsList
unvariedParmsList = [] # parameters in equivalences that are not varied
# scan equivalences: look for holds
for cnum,(varlist,mapvars,multarr,invmultarr) in enumerate(zip(
dependentParmList,indParmList,arrayList,invarrayList)):
if multarr is not None: continue # not an equivalence
if cnum in convertList: continue
# look for holds
gotHold = False
holdList = []
for v in varlist+mapvars:
if v in holdParmList:
gotHold = True
elif type(v) is not float:
holdList.append(v)
if gotHold:
if holdList:
msg = ' Some parameters set as "Hold"; setting remainder as "Hold": '
for i,var in enumerate(holdList):
if i != 0: msg += ", "
msg += var
StoreHold(var,'Equiv fixed')
else:
msg = ' All parameters set as "Hold" '
msg += "\n Will ignore equivalence"
warnEqv(msg,cnum)
removeList.append(cnum)
continue
# look for unvaried parameters
gotVary = False
gotNotVary = False
holdList = []
for v in varlist+mapvars:
if varyList is None or v in varyList:
gotVary = True
holdList.append(v)
elif type(v) is not float:
gotNotVary = True
if v not in unvariedParmsList: unvariedParmsList.append(v)
if gotNotVary: # at least some unvaried parameters
if gotVary: # mix of varied and unvaried parameters
msg = '\nSome parameters not varied; setting remainder as "Hold": '
for i,var in enumerate(holdList):
if i != 0: msg += ", "
msg += var
StoreHold(var,'Equiv fixed')
elif seqHistNum is not None: # don't need to warn for sequential fit
removeList.append(cnum)
continue
else:
msg = 'No parameters varied '
msg += "\n Will ignore equivalence"
warnEqv(msg,cnum)
removeList.append(cnum)
continue
# look for undefined or zero multipliers
holdList = []
drop = 0
for v,m in zip(varlist,invmultarr):
if parmDict is not None and v not in parmDict:
if v not in undefinedVars: undefinedVars.append(v)
if v not in dropVarList: dropVarList.append(v)
drop += 1
elif m == 0:
warnEqv("Parameter "+str(v)+" has a zero multiplier, dropping",cnum)
if v not in dropVarList: dropVarList.append(v)
drop += 1
else:
holdList.append(v)
if drop == len(varlist):
warnEqv("No dependent parameters defined, ignoring equivalence",cnum)
removeList.append(cnum)
continue
for mv in mapvars:
if type(mv) is float: continue
if parmDict is not None and mv not in parmDict:
# independent parameter is undefined, but some dependent parameters are defined
# hold them
if mv not in undefinedVars: undefinedVars.append(mv)
msg = "Parameter(s) "+str(mv)
for v in varlist:
if v in dropVarList:
msg += ', ' + v
msg += "not defined in this refinement\n"
msg = "Setting holds for: "
for i,var in enumerate(holdList):
if i != 0: msg += ", "
msg += var
warnEqv(msg,cnum)
drop += 1
if drop: # independent var and at least one dependent variable is defined
msg = "Dropping undefined parameter(s) "
i = 0
for v in varlist:
if v in dropVarList:
if i != 0: msg += ', '
i += 1
msg += v
warnEqv(msg,cnum)
msg = "Some parameters not defined. Setting holds for: "
for i,var in enumerate(holdList):
if i != 0: msg += ", "
msg += var
StoreHold(var,'Equiv fixed')
warnEqv(msg,cnum)
# Convert equivalences where noted
for cnum,varlist in enumerate(dependentParmList):
if cnum not in convertList: continue
indvar = indParmList[cnum][0]
# msg = '\nChanging equivalence:\n ' + _showEquiv(
# dependentParmList[cnum],indParmList[cnum],invarrayList[cnum])
for dep,mult in zip(dependentParmList[cnum],invarrayList[cnum]):
constrDict += [{indvar:-1.,dep:1./mult[0]}]
fixedList += ['0.0']
#msg += '\n to constraint(s):'
#msg += '\n ' + _FormatConstraint(constrDict[-1],fixedList[-1])
removeList.append(cnum)
# Drop equivalences where noted
global droppedSym
droppedSym = []
if removeList:
for i in set(removeList):
if symGenList[i]:
droppedSym.append((dependentParmList[i],indParmList[i],arrayList[i],invarrayList[i]))
for i in sorted(set(removeList),reverse=True):
del dependentParmList[i],indParmList[i],arrayList[i],invarrayList[i],symGenList[i]
# Drop variables from remaining equivalences
for cnum,varlist in enumerate(dependentParmList):
for j,v in enumerate(varlist):
drop = []
if v in dropVarList:
drop.append(j)
if drop:
for j in sorted(drop,reverse=True):
del indParmList[cnum][j]
return warninfo['msg']
[docs]
def ProcessConstraints(constList,seqmode='use-all',seqhst=None):
"""Interpret the constraints in the constList input into a dictionary, etc.
All :class:`GSASIIobj.G2VarObj` objects are mapped to the appropriate
phase/hist/atoms based on the object internals (random Ids). If this can't be
done (if a phase has been deleted, etc.), the variable is ignored.
If the constraint cannot be used due to too many dropped variables,
it is counted as ignored. In the case of sequential refinements,
the current histogram number is substituted for a histogram number of "*".
NB: this processing does not include symmetry imposed constraints
:param list constList: a list of lists where each item in the outer list
specifies a constraint of some form, as described in the :mod:`GSASIIobj`
:ref:`Constraint definitions <Constraint_definitions_table>`.
:param str seqmode: one of 'use-all', 'wildcards-only' or 'auto-wildcard'.
When seqmode=='wildcards-only' then any constraint with a numerical
histogram number is skipped. With seqmode=='auto-wildcard',
any non-null constraint number is set to the selected histogram.
:param int seqhst: number for current histogram (used for
'wildcards-only' or 'auto-wildcard' only). Should be None for
non-sequential fits.
:returns: a tuple of (constrDict,fixedList,ignored) where:
* constrDict (list of dicts) contains the constraint relationships
* fixedList (list) contains the fixed values for each type
of constraint.
* ignored (int) counts the number of invalid constraint items
(should always be zero!)
"""
constrDict = []
fixedList = []
ignored = 0
namedVarList = []
for constr in constList:
terms = copy.deepcopy(constr[:-3]) # don't change the tree contents
# deal with wildcards in sequential fits
if seqmode == 'wildcards-only' and seqhst is not None:
skip = False
for term in terms:
if term[1].histogram == '*':
term[1] = term[1].varname(seqhst)
elif term[1].histogram:
skip = True
if skip: continue
elif seqmode == 'auto-wildcard' and seqhst is not None:
for term in terms:
term[1] = term[1].varname(seqhst)
elif seqhst is not None:
for term in terms:
if term[1].histogram == '*':
term[1] = term[1].varname(seqhst)
# else:
# term[1] = term[1].varname() # does this change anything???
# separate processing by constraint type
if constr[-1] == 'h':
# process a hold
var = str(terms[0][1])
if '?' not in var:
StoreHold(var,'User supplied')
else:
ignored += 1
elif constr[-1] == 'f':
# process a new variable
fixedList.append(None)
D = {}
varyFlag = constr[-2]
varname = constr[-3]
for term in terms:
var = str(term[1])
if '?' not in var:
D[var] = term[0]
# add extra dict terms for input variable name and vary flag
if varname is None: # no assigned name, create one
global consNum
varname = str(consNum)
consNum += 1
else:
varname = str(varname) # in case this is a G2VarObj
if '::' in varname:
D['_name'] = varname.replace('::','::nv-')
else:
D['_name'] = '::nv-' + varname
D['_name'] = G2obj.MakeUniqueLabel(D['_name'],namedVarList)
D['_vary'] = varyFlag == True # force to bool
constrDict.append(D)
elif constr[-1] == 'c':
# process a constraint equation
D = {}
for term in terms:
var = str(term[1])
if '?' not in var:
D[var] = term[0]
if len(D) >= 1:
fixedList.append(float(constr[-3]))
constrDict.append(D)
else:
ignored += 1
elif constr[-1] == 'e':
# process an equivalence
firstmult = None
eqlist = []
for term in terms:
if term[0] == 0: term[0] = 1.0
var = str(term[1])
if '?' in var: continue
if firstmult is None:
firstmult = term[0]
firstvar = var
else:
eqlist.append([var,firstmult/term[0]])
if len(eqlist) > 0:
StoreEquivalence(firstvar,eqlist,False)
else:
ignored += 1
else:
ignored += 1
return constrDict,fixedList,ignored
[docs]
def StoreHold(var,holdType=None):
'''Takes a variable name and prepares it to be removed from the
refined variables.
Called with user-supplied constraints by :func:`ProcessConstraints`.
At present symGen is not used, but could be set up to track Holds generated
by symmetry.
'''
global holdParmList,holdParmType
if var not in holdParmList:
holdParmList.append(var)
if holdType: holdParmType[var] = holdType
[docs]
def StoreEquivalence(independentVar,dependentList,symGen=True):
'''Takes a list of dependent parameter(s) and stores their
relationship to a single independent parameter (independentVar).
Called with user-supplied constraints by :func:`ProcessConstraints`,
with Pawley constraints from :func:`GSASIIstrIO.GetPawleyConstr`,
with Unit Cell constraints from :func:`GSASIIstrIO.cellVary`
with symmetry-generated atom constraints from :func:`GSASIIstrIO.GetPhaseData`
There is no harm in using StoreEquivalence with the same independent variable::
StoreEquivalence('x',('y',))
StoreEquivalence('x',('z',))
but the same outcome can be obtained with a single call::
StoreEquivalence('x',('y','z'))
The latter will run more efficiently.
Note that mixing independent and dependent variables, such as::
StoreEquivalence('x',('y',))
StoreEquivalence('y',('z',))
is a poor choice. The module will attempt to fix this by transforming the equivalence to a
"Const" constraint.
: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,symGenList
mapList = []
multlist = []
allfloat = True
for var in dependentList:
if isinstance(var, str):
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)
try:
multlist.append(tuple((float(mult),)))
except:
allfloat = False
multlist.append(tuple((mult,)))
# added relationships to stored values
arrayList.append(None)
if allfloat:
invarrayList.append(np.array(multlist))
else:
invarrayList.append(multlist)
indParmList.append(list((independentVar,)))
dependentParmList.append(mapList)
symGenList.append(symGen)
return
[docs]
def SubfromParmDict(s,prmDict):
'''Process a string as a multiplier and convert it to a float value. This
is done by subsituting any GSAS-II parameter names that appear in the
string that have associated values in the parameter dict with the value
for that parameter.
:param str s: a string to be converted to a value
:param dict prmDict: a dictionary with keys as GSAS-II parameter names
and values the corresponding parameter value.
:returns: the evaluated expression as a float.
'''
# TODO: perhaps SubfromParmDict should be called to convert the
# fixed-val in constraint equations from strings to values.
for key in prmDict:
if key in s:
s = s.replace(key,str(prmDict[key]))
return eval(s)
[docs]
def EvaluateMultipliers(constList,*dicts):
'''Convert multipliers for constraints and equivalences that are specified
as strings into values. The strings can specify values in the parameter dicts as
well as normal Python functions, such as "2*np.cos(0::Ax:2/2.)"
:param list constList: a list of dicts containing constraint expressions
:param \\*dicts: one or more dicts containing GSAS-II parameters and their values
can be specified
:returns: an empty string if there were no errors, or an error message listing
the strings that could not be converted.
'''
prmDict = {}
for d in dicts: prmDict.update(d) # combine all passed parameter dicts
problemList = ""
# loop through multipliers in contraint expressions
for const in constList:
for key in const:
if key.startswith('_'): continue
try: # is this already a float, etc?
1+const[key]
continue
except:
pass
try:
newval = SubfromParmDict(const[key][:],prmDict)
if GSASIIpath.GetConfigValue('debug'):
print('Changing ',const[key],'to',newval)
const[key] = newval
except:
if problemList: problemList += ", "
problemList += const[key]
# loop through multipliers in equivalences
global arrayList,invarrayList
for i,(a,valList) in enumerate(zip(arrayList,invarrayList)):
if a is not None: continue # ignore if not equiv
try:
valList.shape
continue # ignore if already a numpy array
except:
pass
repList = []
for v in valList:
try:
1+v[0]
repList.append(tuple((v[0],)))
continue
except:
pass
try:
newval = SubfromParmDict(v[0][:],prmDict)
if GSASIIpath.GetConfigValue('debug'):
print('Changing ',v[0],'to',newval)
repList.append(tuple((newval,)))
except:
if problemList: problemList += ", "
problemList += v[0]
repList.append(tuple(('error',)))
invarrayList[i] = np.array(repList)
return problemList
[docs]
def GetDependentVars(opt=None):
'''Return a list of dependent variables: e.g. parameters that are
constrained in terms of other parameters
:param str opt: type of dependent variables.
'equiv': from equivalences,
'constr': from constraints
None (default): all
:returns: a list of parameter names
'''
if opt == 'equiv':
return constrParms['dep-equiv']
elif opt == 'constr':
return constrParms['dep-constr']
else:
return constrParms['dep-equiv']+constrParms['dep-constr']
[docs]
def GetIndependentVars():
'''Return a list of independent variables: e.g. parameters that are
slaved to other parameters by constraints
:returns: a list of parameter names
'''
return constrParms['indep-equiv']+constrParms['indep-constr']
[docs]
def getConstrError(constrLst,seqmode,seqhst):
'''This is used to display error messages for constraints and
equivalence relations
:parm list constrLst: a single constraint or equivalence as saved
in the data tree
(see :ref:`constraint definitions <Constraint_definitions_table>`).
:param str seqmode: one of 'use-all', 'wildcards-only' or 'auto-wildcard'
:param int seqhst: number for current histogram (used for
'wildcards-only' or 'auto-wildcard' only). Should be None for
non-sequential fits.
:returns: error, msg where error (bool) is True if the
constraint/equivalence creates an error, msg (str) can be a warning
or an error
'''
msg = ''
note = ''
terms = copy.deepcopy(constrLst[:-3])
if seqmode == 'wildcards-only' and seqhst is not None:
if constrLst[-1] == 'e':
msg = 'equivalence'
else:
msg = 'constraint'
for term in terms:
if term[1].histogram == '*':
term[1] = term[1].varname(seqhst)
elif term[1].histogram:
return False,"Ignoring non-wildcard "+msg, "Ignore"
elif seqmode == 'auto-wildcard' and seqhst is not None:
for term in terms:
term[1] = term[1].varname(seqhst)
else:
for term in terms:
if term[1].histogram == '*':
if seqhst is None:
msg = "Parameter "+str(terms[0][1])+" contains a wildcard, which are used only sequential refinements. Constraint ignored."
return False,msg, "Ignored"
else:
term[1] = term[1].varname(seqhst)
else:
term[1] = term[1].varname()
if constrLst[-1] == 'e':
# conflicting uses
if terms[0][1] in constrVarList+convVarList:
msg = "Parameter "+str(terms[0][1])+" used in constraint. To be recast as constraint."
return False,msg, "Recast as constraint."
varList = []
for m,v in terms:
if v in constrVarList+convVarList:
varList.append(str(v))
if varList:
msg = "Parameter(s) used in constraint: "
for i,v in enumerate(varList):
if i != 0: msg += ', '
msg += v
return False,msg,"Recast as constraint."
varList = []
for m,v in terms[1:]:
if v in indepVarList:
varList.append(str(v))
if varList:
msg = "Parameter(s) used elsewhere as independent: "
for i,v in enumerate(varList):
if i != 0: msg += ', '
msg += v
return False,msg,"Recast as constraint."
varList = []
for m,v in terms[1:]:
if v in multdepVarList:
varList.append(str(v))
if varList:
msg += "Parameter(s) repeated as dependent: "
for i,v in enumerate(varList):
if i != 0: msg += ', '
msg += v
return False,msg,"Recast as constraint."
# zero multiplier
varList = []
valid = 0
for m,v in terms:
if m == 0:
varList.append(str(v))
else:
valid += 1
if varList and valid > 1:
msg += "Parameter(s) with zero multipliers: "
for i,v in enumerate(varList):
if i != 0: msg += ', '
msg += v
msg += " will be ignored"
elif varList:
msg += "Parameter(s) with zero multipliers:"
for i,v in enumerate(varList):
if i != 0: msg += ', '
msg += v
return False,msg,"Ignored: not varied"
# hold parameters
s = ''
for m,v in terms:
if v in holdParmList and (
"User supplied" in holdParmType.get(v,'') or
"symmetry" in holdParmType.get(v,'') or
"rigid body" in holdParmType.get(v,'')):
if s: s += ', '
s += str(v)
if s:
if msg: msg += '; '
msg += "Parameter(s) set as Hold: "+s
msg += "\nThis equivalence will be ignored; all parameters are held."
return False,msg,'Has holds: Not varied.'
# unrefined parameters
gotVary = False
gotNotVary = False
s = ''
for m,v in terms:
if v in unvariedParmsList:
gotNotVary = True
if s: s += ', '
s += str(v)
else:
gotVary = True
if gotNotVary and gotVary: # mix of varied and unvaried parameters
if msg: msg += '. '
msg += 'Unvaried parameter(s): '+s+"; remainder set Hold. All parameters fixed."
return False,msg, 'Ignored: not all varied'
elif gotNotVary:
if msg: msg += '. '
msg += 'All parameters not varied. Equivalence Ignored.'
return False,msg, 'Ignored: not varied'
# undefined parameters
undef = 0
s = ''
for m,v in terms[1:]:
if v in undefinedVars:
undef += 1
if s: s += ', '
s += str(v)
if undef == len(terms[1:]):
msg += 'Ignored: None of the dependent parameters are defined'
elif terms[0][1] in undefinedVars:
if s:
s = terms[0][1] + ', ' + s
else:
s = terms[0][1]
msg += 'Undefined parameter(s): '+s+'. Remainder will be fixed'
elif undef:
msg += 'Undefined parameter(s): '+s+' will be dropped'
elif constrLst[-1] == 'h':
v = terms[0][1]
if v in undefinedVars: return False,"Parameter is undefined","Ignored"
if v in unvariedParmsList: return False,"Parameter is not refined","Ignored"
else:
# check for post-grouping errors in new var & constr. eq. groups
for m,v in terms:
if v in groupErrors:
return True,'Constraint singularity: see error listing','Singular'
zeroList = []
toBeUsed = []
undef = []
unvar = []
hold = []
for m,v in terms: # check for zero multiplier, undefined, unvaried or hold
if constrLst[-1] == 'f': # new var expressions;
if m == 0:
zeroList.append(str(v))
elif v in undefinedVars:
undef.append(str(v))
elif (v in holdParmList and constrLst[-2] and
"User supplied" in holdParmType.get(v,'') or
"symmetry" in holdParmType.get(v,'') or
"rigid body" in holdParmType.get(v,'')):
hold.append(str(v))
else: # constraint equation
if m == 0:
zeroList.append(str(v))
elif v in undefinedVars:
undef.append(str(v))
elif v in unvariedParmsList:
unvar.append(str(v))
elif v in holdParmList and holdParmType.get(v,'') != 'dependent param':
hold.append(str(v))
else:
toBeUsed.append(str(v))
s = ''
for v in zeroList:
if s: s += ', '
s += str(v)
if s:
if msg: msg += '; '
msg += "Parameter(s) with zero multipliers: "+s
s = ''
for v in undef:
if s: s += ', '
s += str(v)
if s:
if msg: msg += '; '
msg += "Undefined parameter(s): "+s
s = ''
for v in unvar:
if s: s += ', '
s += str(v)
if s:
if msg: msg += '; '
msg += "Unrefined parameter(s) will be dropped: "+s
s = ''
for v in hold:
if s: s += ', '
s += str(v)
if s:
if msg: msg += '; '
msg += '"Hold" parameter(s): '+s
if hold and constrLst[-1] == 'f':
if msg: msg += '; '
msg += "\nNew var with holds cannot be processed; constraint ignored."
note = 'Ignored'
elif undef and toBeUsed:
s = ''
for v in toBeUsed:
if s: s += ', '
s += str(v)
if msg: msg += '; '
msg += "Adding Holds on "+s+"; Constraint Ignored."
note = 'Ignored'
elif undef or (len(toBeUsed) == 0 and (zeroList or unvar or hold)):
if msg: msg += '; '
msg += "No parameters remain. Constraint Ignored."
note = 'Ignored'
elif len(toBeUsed) == 1 and (zeroList or unvar or hold):
if msg: msg += '; '
msg += "One parameter is retained; converted to fixed value."
note = 'Converted'
elif zeroList or unvar or hold:
if msg: msg += ': '
msg += '\nConstraint adjusted for fixed values and retained.'
note = 'Parameter(s) removed'
return False,msg,note
[docs]
def ComputeDepESD(covMatrix,varyList,noSym=False):
'''Compute uncertainties for dependent parameters from independent ones
returns a dictionary containing the esd values for dependent parameters
:param np.array covMatrix: the full covariance matrix
:param list varyList: the names of the variables matching the columns
and rows in covMatrix
:param bool noSym: When True symmetry generated parameters are
not included. Do this so that redundant s.u.'s eare not shown.
When False (default) s.u. values for all dependent
parameters are placed in the returned dict.
'''
sigmaDict = {}
for varlist,mapvars,multarr,invmultarr,symgen in zip(
dependentParmList,indParmList,arrayList,invarrayList,symGenList):
if symgen and noSym: continue # skip symmetry generted
varied = 0
# 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
varied += 1
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 i,(v,vec) in enumerate(zip(varlist,invmultarr)):
#if i == varied: break # this limits the number of generated params
# to match the number varied. Not sure why I did this.
sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec)))
return sigmaDict
[docs]
def _showEquiv(varlist,mapvars,invmultarr,longmsg=False):
'''Format an equivalence relationship, note that varlist, mapvars, invmultarr
are elements of dependentParmList, indParmList, invarrayList
'''
for i,mv in enumerate(mapvars):
s1 = str(mv)
if not longmsg:
s1 += ' ==> '
elif len(varlist) == 1:
s1 += ' is equivalent to '
else:
s1 += ' is equivalent to parameters: '
j = 0
for v,m in zip(varlist,invmultarr):
if debug: print ('v,m[0]: ',v,m[0])
if len(s1.split('\n')[-1]) > 60: s1 += '\n '
if j > 0: s1 += ' & '
j += 1
s1 += str(v)
if m != 1:
s1 += " / " + str(m[0])
return s1
[docs]
def VarRemapShow(varyList=None,inputOnly=False,linelen=60):
'''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
'''
if varyList is None:
varyList = saveVaryList
s = ''
depVars = constrParms['dep-equiv']+constrParms['dep-constr']
if len(depVars) > 0:
s += '\nDependent parameters from constraints & equivalences:\n'
for v in sorted(set(depVars)):
s += ' ' + v + '\n'
first = True
for v in sorted(holdParmList):
if v in depVars: continue
if v in holdParmType: v += '\t'+holdParmType[v]
if first:
s += '\nParameters set as Hold:\n'
first = False
s += ' ' + v + '\n'
userOut = ''
symOut = ''
consOut = ''
varOut = ''
freeOut = ''
global dependentParmList,arrayList,invarrayList,indParmList,symGenList
for varlist,mapvars,multarr,invmultarr,symFlag in zip(
dependentParmList,indParmList,arrayList,invarrayList,symGenList):
for i,mv in enumerate(mapvars):
if multarr is None:
# s1 = ' ' + str(mv) + ' is equivalent to parameter(s): '
if len(varlist) == 1:
s1 = ' ' + str(mv) + ' is equivalent to '
else:
s1 = ' ' + str(mv) + ' is equivalent to parameters: '
j = 0
for v,m in zip(varlist,invmultarr):
if debug: print ('v,m[0]: ',v,m[0])
if j > 0: s1 += ' & '
j += 1
s1 += str(v)
if m != 1:
s1 += " / " + '{:.4f}'.format(m[0])
#if len(s1.split('\n')[-1]) > 70:
# s1 = ' \n &'.join(s1.rsplit('&',1))
if symFlag:
symOut += s1 + '\n'
else:
userOut += s1 + '\n'
continue
ln = ''
# if mv in varyList:
# lineOut = ' (V)* {} = '.format(mv)
# else:
# lineOut = ' {} = '.format(mv)
j = 0
lineOut = ' {} = '.format(mv)
for (m,v) in zip(multarr[i,:],varlist):
if m == 0: continue
if m < 0:
lineOut += ' - '
m *= -1
elif j != 0:
lineOut += ' + '
j += 1
if len(lineOut) > linelen:
ln += lineOut
lineOut = '\n '
if m == 1:
lineOut += '{}'.format(v)
else:
lineOut += '({:.4g} * {})'.format(m,v)
if mv in varyList:
lineOut += '\t *VARIED*'
ln += lineOut
if type(mv) is float:
consOut += ln + '\n'
elif '::nv-' in mv:
varOut += ln + '\n'
else:
freeOut += ln + '\n'
if userOut:
s += '\nUser-supplied equivalences:\n' + userOut
if droppedSym:
for varlist,mapvars,multarr,invmultarr in droppedSym:
for i,mv in enumerate(mapvars):
if len(varlist) == 1:
s1 = ' ' + str(mv) + ' is equivalent to '
else:
s1 = ' ' + str(mv) + ' is equivalent to parameters: '
j = 0
for v,m in zip(varlist,invmultarr):
if debug: print ('v,m[0]: ',v,m[0])
if j > 0: s1 += ' & '
j += 1
s1 += str(v)
if m != 1:
s1 += " / " + '{:.4f}'.format(m[0])
symOut += s1 + '\t *recast as equation*\n'
if symOut:
s += '\nSymmetry-generated equivalences:\n' + symOut
if consOut:
s += '\nConstraint Equations:\n' + consOut
if varOut:
s += '\nNew Variable assignments:\n' + varOut
if freeOut:
s += '\nGenerated constraint equations:\n' + freeOut
if not (userOut or consOut or varOut or symOut):
return s + '\nNo constraints or equivalences in use'
elif inputOnly:
return s
s += '\nInverse parameter mapping relations:\n'
lineDict = {} # store so we can sort them
for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
varied = False
for i,mv in enumerate(varlist):
constrVal = 0 # offset to constraint from constant terms
lineOut = ' {} = '.format(mv)
j = 0 # number of terms shown
for m,v in zip(invmultarr[i,:],mapvars):
if np.isclose(m,0): continue
try:
if np.isclose(float(v),0): continue
except:
pass
if v in varyList: varied = True
try:
constrVal += m*v # v is a constant; add to offset
continue
except:
pass
if m < 0:
lineOut += ' - '
m *= -1
elif j != 0:
lineOut += ' + '
j += 1
if len(lineOut) > linelen:
s += lineOut
lineOut = '\n '
if m == 1:
lineOut += '{}'.format(v)
else:
lineOut += '({:.4g} * {})'.format(m,v)
if constrVal < 0:
lineOut += ' - {:.4g}'.format(-constrVal)
elif constrVal > 0:
lineOut += ' + {:.4g}'.format(constrVal)
elif j == 0:
lineOut += '0' # no terms, no constants: var fixed at zero
if varied: lineOut += '\t *VARIED*'
lineDict[mv] = lineOut
for key in sorted(lineDict):
s += lineDict[key] + '\n'
return s
[docs]
def getInvConstraintEq(var,varyList):
'''For a dependent variable, find the constraint that
defines the dependent variable in terms of varied independent variables.
This works for constraint equations (via new var or generated parameters)
or equivalences. For equivalences the result will lists of length 1
:param str var: named of refined variable (e.g. 0:0:Scale)
:param list varyList: list of refined variables
:returns: vList,mList where vList is a list of variables and
mList is a list of multipliers for that variable (floats)
'''
for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
if var not in varlist: continue
i = varlist.index(var)
vList = []
mList = []
for m,v in zip(invmultarr[i,:],mapvars):
if v not in varyList: continue
if m == 0: continue
if v == 0: continue
vList.append(v)
mList.append(m)
return vList,mList
if GSASIIpath.GetConfigValue('debug'): print('getInvConstraintEq: not found: ',var)
return [],[] # unexpected -- not an independent parameter
[docs]
def GetSymEquiv(seqmode,seqhistnum):
'''Return the automatically generated (equivalence) relationships.
:returns: a list of strings containing the details of the contraint relationships
'''
symout = []
symerr = []
symhelp = []
global dependentParmList,arrayList,invarrayList,indParmList,symGenList
for varlist,mapvars,multarr,invmultarr,symFlag in zip(
dependentParmList,indParmList,arrayList,invarrayList,symGenList):
if not symFlag: continue
for i,mv in enumerate(mapvars):
cnstr = [[1,G2obj.G2VarObj(mv)]]
if multarr is None:
s1 = ''
s2 = ' = ' + str(mv)
j = 0
helptext = 'Variable {:} '.format(mv) + " ("+ G2obj.fmtVarDescr(mv) + ")"
if len(varlist) == 1:
cnstr.append([invmultarr[0][0],G2obj.G2VarObj(varlist[0])])
# format the way Bob prefers
if invmultarr[0][0] == 1:
s1 = str(varlist[0]) + ' = ' + str(mv)
else:
s1 = str(varlist[0]) + ' = ' + str(
invmultarr[0][0]) + ' * '+ str(mv)
s2 = ''
m = 1./invmultarr[0][0]
var1 = str(varlist[0])
helptext += "\n\nis equivalent to "
if m == 1:
helptext += '\n {:} '.format(var1) + " ("+ G2obj.fmtVarDescr(var1) + ")"
else:
helptext += '\n {:3g} * {:} '.format(m,var1) + " ("+ G2obj.fmtVarDescr(var1) + ")"
else:
helptext += "\n\nis equivalent to the following:"
for v,m in zip(varlist,invmultarr):
cnstr.append([m,G2obj.G2VarObj(v)])
#if debug: print ('v,m[0]: ',v,m[0])
if len(s1.split('\n')[-1]) > 75: s1 += '\n '
if j > 0: s1 += ' = '
j += 1
s1 += str(v)
if m != 1:
s1 += " / " + str(m[0])
helptext += '\n {:3g} * {:} '.format(m,v) + " ("+ G2obj.fmtVarDescr(v) + ")"
else:
helptext += '\n {:} '.format(v) + " ("+ G2obj.fmtVarDescr(v) + ")"
err,msg,note = getConstrError(cnstr+[None,None,'e'],seqmode,seqhistnum)
symerr.append([msg,note])
symout.append(s1+s2)
symhelp.append(helptext)
else:
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)
print ('unexpected sym op='+s)
return symout,symerr,symhelp
[docs]
def GetDroppedSym(seqmode,seqhistnum):
'''Return automatically generated (equivalence) relationships that were
converted to constraint equations
:returns: a list of strings containing the details of the equivalences
'''
symout = []
symerr = []
symhelp = []
global droppedSym
for varlist,mapvars,multarr,invmultarr in droppedSym:
for i,mv in enumerate(mapvars):
cnstr = [[1,G2obj.G2VarObj(mv)]]
if multarr is None:
s1 = ''
s2 = ' = ' + str(mv)
j = 0
helptext = 'Variable {:} '.format(mv) + " ("+ G2obj.fmtVarDescr(mv) + ")"
if len(varlist) == 1:
cnstr.append([invmultarr[0][0],G2obj.G2VarObj(varlist[0])])
# format the way Bob prefers
if invmultarr[0][0] == 1:
s1 = str(varlist[0]) + ' = ' + str(mv)
else:
s1 = str(varlist[0]) + ' = ' + str(
invmultarr[0][0]) + ' * '+ str(mv)
s2 = ''
m = 1./invmultarr[0][0]
var1 = str(varlist[0])
helptext += "\n\nis equivalent to "
if m == 1:
helptext += '\n {:} '.format(var1) + " ("+ G2obj.fmtVarDescr(var1) + ")"
else:
helptext += '\n {:3g} * {:} '.format(m,var1) + " ("+ G2obj.fmtVarDescr(var1) + ")"
else:
helptext += "\n\nis equivalent to the following:"
for v,m in zip(varlist,invmultarr):
cnstr.append([m,G2obj.G2VarObj(v)])
#if debug: print ('v,m[0]: ',v,m[0])
if len(s1.split('\n')[-1]) > 75: s1 += '\n '
if j > 0: s1 += ' = '
j += 1
s1 += str(v)
if m != 1:
s1 += " / " + str(m[0])
helptext += '\n {:3g} * {:} '.format(m,v) + " ("+ G2obj.fmtVarDescr(v) + ")"
else:
helptext += '\n {:} '.format(v) + " ("+ G2obj.fmtVarDescr(v) + ")"
err,msg,note = getConstrError(cnstr+[None,None,'e'],seqmode,seqhistnum)
symerr.append([msg,note])
symout.append(s1+s2)
symhelp.append(helptext)
else:
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)
print ('unexpected sym op='+s)
return symout,symerr,symhelp
[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 parameters: 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:
if debug: print ('start dMdv for',name,dMdv[varyList.index(name)])
for v,m in zip(varlist,invmultarr):
if m[0] == 0: continue
dMdv[varyList.index(name)] += derivDict[v]/ m[0]
else:
for v,m in zip(varlist,invmultarr[:,i]):
if m == 0: continue
dMdv[varyList.index(name)] += m * derivDict[v]
[docs]
def Map2Dict(parmDict,varyList):
'''Updates the parameter dictionary and the varyList using the
equivalence and constraint input. This should be called at least once, after
the constraints have been defined using :func:`StoreEquivalence`,
:func:`GroupConstraints` and :func:`GenerateConstraints` and before any
parameter refinement is done.
This completes the parameter dictionary by defining values for parameters
created by constraints based on the constraints that define them
using the values for the current parameters. It also removes all dependent
variables from the varyList
:param dict parmDict: a dict containing parameter values keyed by the
parameter names. For new variables created by constraints, entries
will be added to the dictionary, if not alreay present, or the
values will be recomputed.
:param list varyList: a list of parameters names. Will be modified.
'''
# remove fixed parameters from the varyList
for item in holdParmList:
if varyList is not None and item in varyList: varyList.remove(item)
# process the independent parameters:
# * remove dependent ones from varylist
# * for equivalences apply the independent parameters onto dependent variables
global dependentParmList,arrayList,invarrayList,indParmList
for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
for item in varlist: # TODO: is this still needed?
if varyList is not None and item in varyList: varyList.remove(item)
if multarr is None:
#for v,val in zip( # shows values to be set
# varlist,
# np.dot(invmultarr,np.array([parmDict[var] for var in mapvars]))
# ): print('parmDict set',v,':',val)
parmDict.update(zip(varlist,np.dot(invmultarr,np.array([parmDict[var] for var in mapvars]))))
# * for the created parameters, compute them from their dependents
for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
if multarr is None: continue
# evaluate constraints in the forward direction
A = np.array([parmDict[var] for var in varlist])
z = zip(mapvars,np.dot(multarr,A))
# add/replace in parameter dict
parmDict.update([i for i in z if type(i[0]) is not float and ':' in i[0]])
# parmDict.update([i for i in zip(mapvars,np.dot(multarr,A)) if ':' in i[0]])
global saveVaryList
if varyList is not None:
saveVaryList = copy.copy(varyList)
else:
saveVaryList = []
[docs]
def normParms(parmDict):
'''Attempt to put parameters into the right ballpark by scaling to
enforce constraint equations
'''
for varlist,mapvars,multarr,invmultarr in zip(
dependentParmList,indParmList,arrayList,invarrayList):
if multarr is None or invmultarr is None: continue # unexpected
for i,s in enumerate(mapvars):
try:
s = float(s)
except:
continue
sumcons = 0.
for var,m in zip(varlist,multarr[i]):
if var not in parmDict:
print('normParms error: Parameter',var,'not in parmDict')
break
sumcons += parmDict[var] * m
if sumcons != s and sumcons != 0 and s != 0:
for var in varlist:
parmDict[var] *= s/sumcons
[docs]
def Dict2Map(parmDict):
'''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 after refinement and before the parameters are used for
any computations
:param dict parmDict: a dict containing parameter values keyed by the
parameter names. After this is called, all the dependent variables
will be updated based on constraints and equivalences.
'''
global dependentParmList,arrayList,invarrayList,indParmList
for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
if invmultarr is None: # is this needed?
if GSASIIpath.GetConfigValue('debug'):
print('Why does this constraint have None for invmultarr?',varlist,mapvars)
continue
valslist = np.array([float(parmDict.get(var,var)) for var in mapvars])
#for v,val in zip(varlist,np.dot(invmultarr,np.array(valslist))): print(v,val) # shows what is being set
parmDict.update(zip(varlist,np.dot(invmultarr,valslist)))
#======================================================================
# 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 ConstraintException("Singular input to GramSchmidtOrtho")
a[j] /= np.linalg.norm(a[j])
return a
[docs]
def _FillArray(sel,d,collist):
'''Construct a n by n matrix [n = len(collist)]
with the initial m rows [m = len(sel)] using the
relationships defined in the expressions dict, d.
Since m may be smaller than n, the remaining rows
are filled with rows that are tested to not create
a singular matrix.
:param list sel: a list of indices in dict d
:param list d: a list of dict's where each dict describes an
expression from a constraint equation or a new var
:param list collist: a list parameter names.
:returns: an n by n numpy.array matrix
'''
n = len(collist)
m = len(sel)
arr = np.zeros(2*[n,])
# fill the top rows
for i,cnum in enumerate(sel):
for j,var in enumerate(collist):
arr[i,j] = d[cnum].get(var,0)
try:
_RowEchelon(m,copy.copy(arr),collist[:])
except:
raise Exception('Initial constraints singular')
for i in range(m,n):
arr[i][i] = 1 # add a diagonal element
try:
_RowEchelon(i+1,copy.copy(arr),collist[:])
continue
except:
pass
for j in range(n):
if j == i: continue
arr[i][j] = 1 # add another element
try:
_RowEchelon(i+1,copy.copy(arr),collist[:])
break
except:
pass
arr[i][j] = -1 # try a different valuefor this element
try:
_RowEchelon(i+1,copy.copy(arr),collist[:])
break
except:
arr[i][j] = 0 # reset to add another element
else:
raise Exception('Unable to create non-singular matrix')
return arr
[docs]
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 ConstraintException('Singular input')
[docs]
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
'''
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__":
d = [{'a': 0, 'b': 1.5,'d':0}, {'d': -1}]
lbls = ['a','c','b','d']
sel = (0,1)
try:
arr2 = _FillArray(sel,d,lbls)
except Exception as err:
if 'Initial' in str(err):
print('initial error')
else:
print('unable to extend matrix error')
import sys; sys.exit()
print(arr2)
d = [{'a': 1, 'b': 1,}]
lbls = ['a','b']
sel = (0,)
arr1 = _FillArray(sel,d,lbls)
print(arr1)
print(GramSchmidtOrtho(arr1,1))
import sys; sys.exit()
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']
# 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']
msg,warning,groups,parmlist = GenerateConstraints(varyList,constrDict,fixedList,parmdict)
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(list(parmdict.keys())):
print (' '+key,'\t',before.get(key),'\t',parmdict[key])
print ('varylist after',varyList)
before = parmdict.copy()
Dict2Map(parmdict)
print ('after Dict2Map')
print (' key / before / after')
for key in sorted(list(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)