source: trunk/GSASIImapvars.py @ 3713

Last change on this file since 3713 was 3713, checked in by toby, 3 years ago

constraint bug fix; comments

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 78.0 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2018-10-31 04:53:13 +0000 (Wed, 31 Oct 2018) $
4# $Author: toby $
5# $Revision: 3713 $
6# $URL: trunk/GSASIImapvars.py $
7# $Id: GSASIImapvars.py 3713 2018-10-31 04:53:13Z toby $
8########### SVN repository information ###################
9"""
10*GSASIImapvars: Parameter constraints*
11======================================
12
13Module to implements algebraic contraints, parameter redefinition
14and parameter simplification contraints.
15
16Types of constraints
17--------------------
18
19There are four ways to specify constraints, as listed below.
20Note that parameters are initially stored in the
21main section of the GSAS-II data tree under heading ``Constraints``.
22This dict has four keys, 'Hist', 'HAP', 'Global', and 'Phase',
23each containing a list of constraints. An additional set of constraints
24are generated for each phase based on symmetry considerations by calling
25:func:`GSASIIstrIO.GetPhaseData`.
26
27Note that in the constraints, as stored in the GSAS-II data tree, parameters
28are stored as :class:`GSASIIobj.G2VarObj` objects, as these objects allow for
29changes in numbering of phases, histograms and atoms.
30When they are interpreted (in :func:`GSASIIstrIO.ProcessConstraints`),
31references
32to numbered objects are resolved using the appropriate random ids and the
33variable object is expressed as a string of form ``ph:hst:VARNAM:at``.
34
35Alternate parameters (New Var)
36^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
37
38Parameter redefinition ("New Var" constraints)
39is done by creating an expression that relates several
40parameters:
41
42::
43
44   Mx1 * Px + My1 * Py +...
45   Mx2 * Px + Mz2 * Pz + ...
46
47where Pj is a GSAS-II parameter name and Mjk is a constant (float) multiplier.
48Alternately, multipliers Mjk can contain a formula (str) that will be evaluated prior
49to the start of the refinement. In a formula, GSAS-II parameters will be replaced by the
50value of the parameter before the formula is evaluated, so ``'np.cos(0::Ax:2)'`` is a valid
51multiplier. At present, only phase (atom/cell) parameters are available for use in
52a formula, but this can be expanded if needed.
53
54This type of constraint describes an alternate
55degree of freedom where parameter Px and Py, etc. are varied to keep
56their ratio
57fixed according the expression. A new variable is assigned to each degree of
58freedom when refined. An example where this can be valuable is when
59two parameters, P1 and P2, have similar values and are highly correlated. It is often better to create a new variable, Ps = P1 + P2, and refine Ps.
60In the later stages of refinement, a second
61variable, Pd = P1 - P2, can be defined and it can be seen if refining Pd is
62supported by the data. Another use will be to define parameters that
63express "irrep modes" in terms of the fundamental structural parameters.
64
65These "New Var" constraints are stored as described for type "f" in the
66:ref:`constraint definitions table <Constraint_definitions_table>`.
67
68Constrained parameters (Const)
69^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
70
71A constraint on a set of variables can be supplied in the form of a
72linear algebraic equation: ::
73
74  Nx1 * Px + Ny1 * Py +... = C1
75
76where Cn is a constant (float), where Pj is a GSAS-II parameter name,
77and where Njk is a constant multiplier (float) or a formula (str) that will be evaluated prior
78to the start of the refinement. In a formula, GSAS-II parameters will be replaced by the
79value of the parameter before the formula is evaluated, so ``'np.cos(0::Ax:2)'`` is a valid
80multiplier. At present, only phase (atom/cell) parameters are available for use in
81a formula, but this can be expanded if needed.
82
83These equations set an interdependence between variables.
84Common uses of parameter constraints are to set rules that decrease the number of parameters,
85such as restricting the sum of fractional occupancies for atoms that share
86a site to sum to unity, thus reducing the effective number of variables by one.
87Likewise, the Uiso value for a H atom "riding" on a C, N or O atom
88can be related by a fixed offset (the so called B+1 "rule").
89
90A "Const" constraint is stored as described for type "c" in the
91:ref:`constraint definitions table <Constraint_definitions_table>`.
92
93Equivalenced parameters (Equiv)
94^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
95
96A simplifed way to set up a constraint equation is to define an equivalence,
97which can be of form: ::
98
99  C1 * P1 = C2 * Py
100
101or::
102
103  C1 * P1 = C2 * P2 = C3 * P3 = ...
104
105where Cn is a constant (float), where Pj is a GSAS-II parameter name. This
106means that parameters Py (or P2 and P3) are determined from (or "slaved" to)
107parameter P1. Alternately, equivalences can be created with :func:`StoreEquivalence`
108where the multipliers can be a formula (str) that will be evaluated prior to the start of
109the refinement. In a formula, GSAS-II parameters will be replaced by the value of the
110parameter before the formula is evaluate, so ``'np.cos(0::Ax:2)'`` is a valid multiplier.
111At present, only phase (atom/cell) parameters are available for use in
112a formula, but this can be expanded if needed.
113Note that
114the latter constraint expression is conceptually identical to
115defining constraint equations. In practice, however,
116equivalenced parameters are processed in a
117different and more direct manner than constraint equations. The previous
118set of equalities could also be written in this way as a set of constraint
119equations: ::
120
121  C1 * P1 - C2 * P2 = 0
122  C1 * P1 - C3 * P3 = 0
123  ...
124
125
126The first parameter (P1 above)
127is considered the independent variable
128and the remaining parameters are dependent variables. The dependent variables
129are set from the independent variable.
130An example of how this may be used woul be if, for example,
131a material has a number of O atoms, all in fairly similar bonding environments
132and the diffraction data are sparse, one my reduce the complexity of the model
133by defining Uiso for the first O atoms to be identical to the remaining atoms.
134The results of this refinement will be simpler to understand than if a set of
135constraint equations is used because the refined parameter will be the independent variable, which will be as ph::Uiso:n, corresponding to the first O atom.
136
137A parameter can be used in multiple
138equivalences as independent variable,
139but if parameter is used as both a dependent and independent variable or
140a parameter is used in equivalences and in "New Var" or "Const" constraints,
141this create conflicts that cannot be resolved within the equivalences implementation
142but can be handled as constraint equations.
143The equivalences that violate this are discovered in :func:`CheckEquivalences`
144and then :func:`MoveConfEquiv` is used to change these equivalences to "Const" equations.
145
146Equivalenced parameters ("EQUIV" constraints), when defined by users,
147are stored as described for type "e" in the
148:ref:`constraint definitions table <Constraint_definitions_table>`.
149Other equvalences are generated by symmetry prior
150to display or refinement in :func:`GSASIIstrIO.GetPhaseData`.
151These are not stored.
152
153Fixed parameters (Hold)
154^^^^^^^^^^^^^^^^^^^^^^^^
155
156When variables are refined where a single refinement flag determines that several variables
157are refined at the same time (examples are: cell parameters, atom positions, anisotropic
158displacement parameters, magnetic moments,...) it can be useful to specify that a
159specific parameter should not be varied. These will most commonly be generated due to symmetry,
160but under specific conditions, there may be other good reasons to constrain a variable.
161
162A "Hold" constraint is stored as described for type "h" in the
163:ref:`constraint definitions table <Constraint_definitions_table>`.
164
165Constraint Processing
166---------------------
167
168When constraints will be used or edited, they are processed using a series of
169calls:
170
171* First all of the stored constraints are appended into a single list. They are
172  initially stored in separate lists only to improve their creation and display
173  in the GUI.
174
175* Then :func:`InitVars` is used to initialize the global variables in
176  this module (:mod:`GSASIImapvars`).
177
178* Then :func:`GSASIIstrIO.ProcessConstraints` is used to initially process the
179  constraints, as described below.
180
181* Symmetry-generated equivalences are then created in
182  :func:`GSASIIstrIO.GetPhaseData`, which also calls
183  :func:`GSASIIstrIO.cellVary` and for Pawley refinements
184  :func:`GSASIIstrIO.GetPawleyConstr`. These are entered directly into this
185  module's globals using :func:`StoreEquivalence`.
186* Constraints/equivalences are then checked for possible conflicts with
187  :func:`CheckConstraints`, this requires grouping the constraints,
188  as described below.
189* In refinements, :func:`GenerateConstraints` is then called to
190  create the constraints that will be used, see below for
191* For debugging constraints, :func:`VarRemapShow` can be called after
192  :func:`GenerateConstraints` to display the generated constraints.
193
194Constraint Reorganization (:func:`~GSASIIstrIO.ProcessConstraints`)
195^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
196
197:func:`GSASIIstrIO.ProcessConstraints` is used to initially process the
198constraints. This does these things:
199
2001. The "Hold", "Const" and "New Var" expressions are split between two paired lists,
201   :data:`constDictList` and :data:`fixedList` which are set:
202 
203   * For "Hold" entries a dict with a single entry is placed in constDictList where
204     the key is the parameter name (associated value is 0.0) and fixedList gets a
205     value of 0.0.
206   * For "Const" entries, a dict with multiple entries is placed in constDictList where
207     the key is the parameter name and the value is the multiplier for the parameter,
208     while fixedList gets a string value corresponding to the constant value for
209     the expression.
210   * For "New Var" entries, a dict with multiple entries is placed in constDictList
211     where the key is the parameter name and the value is the multiplier
212     for the parameter; an additional key "_vary" is given the value of True or False
213     depending on the refinement flag setting. The corresponding entry in
214     fixedList is None
215
216   The output from this will have this form where the first entry is a "Const", the
217   second is a "New Var" and the final is a "Hold".
218
219  .. code-block:: python
220
221    constrDict = [
222         {'0:12:Scale': 2.0, '0:14:Scale': 4.0, '0:13:Scale': 3.0, '0:0:Scale': 0.5},
223         {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0, '_vary':True},
224         {'0::A0': 0.0}]
225    fixedList = ['5.0', None, '0']
226
227
2282. Equivalences are stored using :func:`StoreEquivalence`
229   into this module's globals (:data:`~GSASIImapvars.arrayList`,
230   :data:`~GSASIImapvars.invarrayList`,    :data:`~GSASIImapvars.indParmList`,
231   :data:`~GSASIImapvars.dependentParmList` and  :data:`~GSASIImapvars.symGenList`)
232
233
234Parameter Grouping (:func:`GenerateConstraints`)
235^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
236
237Functions :func:`CheckConstraints` and :func:`GenerateConstraints` are used to
238process the parameter equivalences and constraint lists created in
239:func:`~GSASIIstrIO.ProcessConstraints`. The former is used to generate
240error messages and the latter to generate the internal information used to
241apply the constraints.
242
243Initially, in both a list of parameters that are fixed and those used in
244constraint relations are tabulated in :func:`CheckEquivalences`. The equivalence
245relations are the scanned for the following potential problems:
246
2471. a parameter is used as a dependent variable in more than one
248   equivalence relation
2492. a parameter is fixed and used in in an equivalence relation either
250   as a dependent or independent variable
2513. a parameter is used as a dependent variable in one equivalence relation and
252   as a independent variable in another
2534. a parameter is used in in an equivalence relation (either
254   as a dependent or independent variable) and is used in a constraint
255   expression
2565. a parameter is not defined in a particular refinement, but is used in an
257   equivalence relation
2586. a parameter uses a wildcard for the histogram number (sequential refinements)
259
260Cases 1 & 2 above cannot be corrected, and result in errors. Cases 3 & 4 are
261potentially corrected with :func:`MoveConfEquiv`, as described below.
262Case 5 causes the equivalence to
263be dropped. Case 6 causes the current histogram number to be substituted for
264the wildcard.
265
266For cases 3 & 4, :func:`MoveConfEquiv` is used to change these equivalences
267into "Const" equations. This can potentially mean that additional
268equivalences will be problematic, so if there are changes made by
269:func:`MoveConfEquiv`, :func:`CheckEquivalences` is repeated. If any problem
270cases are noted, the refinement cannot be performed.
271
272Constraint expressions ("Const" and "New Var") are sorted into
273groups so that each group contains the minimum number of entries that
274ensures each parameter is referenced in only one group in
275:func:`GroupConstraints`. This is done by scanning the
276list of dicts in :data:`constDictList` one by one and making a list
277of parameters used in that constraint expression. Any expression that contains
278a parameter in is in that list is added to the current group and those
279parameters are added to this list of parameters. The list of ungrouped
280expressions is then scanned again until no more expressions are added to the
281current group. This process is repeated until every expression has been
282placed in a group. Function :func:`GroupConstraints` returns two lists of lists.
283The first has, for each group, a list of the indices in :data:`constDictList`
284that comprise the group (there can be only one). The second list contains,
285for each group, the unique parameter names in that group.
286
287Each constraint group is then processed. First, wildcard parameters are
288renamed (in a sequential refinement). Any fixed parameters that are used
289in constraints are noted as errors. The number of refined parameters and
290the number of parameters that are not defined in the current refinement are
291also noted. It is fine if all parameters in a group are not defined or all are
292not varied, but if some are defined and others not or some are varied and
293others not, this constitutes an error.
294
295The contents of each group is then examined. Groups with a single
296parameter (holds) are ignored. Then for each group, the number
297of parameters in the group (Np) and the number of expressions in the
298group (Nc) are counted and for each expression. If Nc > Np, then the constraint
299is overdetermined, which also constitutes an error.
300
301The parameter multipliers for each expression are then assembled:
302
303::
304
305   M1a * P1 + M2a * P2 +... Mka * Pk
306   M1b * P1 + M2b * P2 +... Mkb * Pk
307   ...
308   M1j * P1 + M2j * P2 +... Mkj * Pk
309
310
311From this it becomes possible to create a Nc x Np matrix, which is
312called the constraint matrix:
313
314 .. math::
315
316   \\left( \\begin{matrix}
317   M_{1a}  & M_{2a} &... & M_{ka} \\\\
318   M_{1b}  & M_{2b} &... & M_{kb} \\\\
319   ...  \\\\
320   M_{1j}  & M_{2j}  &... & M_{kj}
321   \\end{matrix}\\right)
322
323When Nc<Np, then additional rows need to be added to the matrix and to
324the vector that contains the value for each row (:data:`fixedList`) where
325values are ``None`` for New Vars and a constant for fixed values.
326This then can describe a system of Np simultaneous equations:
327
328 .. math::
329
330   \\left( \\begin{matrix}
331   M_{1a}  & M_{2a} &... & M_{ka} \\\\
332   M_{1b}  & M_{2b} &... & M_{kb} \\\\
333   ...  \\\\
334   M_{1j}  & M_{2j}  &... & M_{kj}
335   \\end{matrix}\\right)
336   \\left( \\begin{matrix}
337   P_{1} \\\\
338   P_{2} \\\\
339   ...  \\\\
340   P_{k}
341   \\end{matrix}\\right)
342   =
343   \\left( \\begin{matrix}
344   C_{1} & C_{2} &  ... & C_{k}
345   \\end{matrix}\\right)
346
347This is done by creating a square matrix from the group using
348:func:`_FillArray` with parameter ``FillDiagonals=False`` (the default). Any
349unspecified rows are left as all zero. The first Nc rows in the
350array are then coverted to row-echelon form using :func:`_RowEchelon`. This
351will create an Exception if any two rows are linearly dependent (which means
352that no matter what values are used for the remaining rows, that the matrix
353will be singular). :func:`_FillArray` is then called with parameter
354``FillDiagonals=True``, which again creates a square matrix but where
355unspecified rows are zero except for the diagonal elements. The 
356`Gram-Schmidt process <http://en.wikipedia.org/wiki/Gram-Schmidt>`_,
357implemented  in :func:`GramSchmidtOrtho`, is used to find orthonormal unit
358vectors for the remaining Np-Nc rows of the matrix. This will fail with
359a ``ConstraintException`` if this is not possible (singular matrix) or
360the result is placed in :data:`constrArr` as a numpy array.
361
362Rows in the matrix corresponding to "New Var" constraints and those that
363were generated by the Gram-Schmidt process are provided with variable names
364(this can be specified if a "New Var" entry by using a ``"_name"`` element
365in the constraint dict, but at present this is not implemented.) Names are
366generated using :data:`paramPrefix` which is set to ``"::constr"``, plus a
367number to make the new parameter name unique. Global dict :data:`genVarLookup`
368provides a lookup table, where the names of the parameters related to this new
369parameter can be looked up easily.
370
371Finally the parameters used as input to the constraint are placed in
372this module's globals
373:data:`~GSASIImapvars.dependentParmList` and the constraint matrix is
374placed in into  :data:`~GSASIImapvars.arrayList`. This can be used to compute
375the initial values for "New Var" parameters. The inverse of the
376constraint matrix is placed in :data:`~GSASIImapvars.invarrayList` and a list
377of the "New Var" parameters and a list of the fixed values (as str's)
378is placed in :data:`~GSASIImapvars.indParmList`. A lookup table for
379fixed values as floats is placed in :data:`~GSASIImapvars.fixedDict`.
380Finally the appropriate entry in :data:`~GSASIImapvars.symGenList` is set to
381False to indicate that this is not a symmetry generated constraint.
382
383
384*Externally-Accessible Routines*
385---------------------------------
386
387To define a set of constrained and unconstrained relations, one
388defines a list of dictionary defining constraint parameters and their
389values, a list of fixed values for each constraint and a list of
390parameters to be varied. In addition, one uses
391:func:`StoreEquivalence` to define parameters that are equivalent.
392Use :func:`EvaluateMultipliers` to convert formula-based constraint/equivalence
393multipliers to numbers and then
394use :func:`CheckConstraints` to check that the input is
395internally consistent and finally :func:`GroupConstraints` and
396:func:`GenerateConstraints` to generate the internally used
397tables. Routine :func:`Map2Dict` is used to initialize the parameter
398dictionary and routine :func:`Dict2Map`, :func:`Dict2Deriv`, and
399:func:`ComputeDepESD` are used to apply constraints. Routine
400:func:`VarRemapShow` is used to print out the constraint information,
401as stored by :func:`GenerateConstraints`. Further information on each routine
402is below:
403
404:func:`InitVars`
405  This is optionally used to clear out all defined previously defined constraint information
406 
407:func:`StoreEquivalence`
408  To implement parameter redefinition, one calls StoreEquivalence. This should be called for every set of
409  equivalence relationships. There is no harm in using StoreEquivalence with the same independent variable:
410
411  .. code-block:: python
412
413       StoreEquivalence('x',('y',))
414       StoreEquivalence('x',('z',))
415
416  or equivalently
417
418  .. code-block:: python
419
420       StoreEquivalence('x',('y','z'))
421
422  The latter will run more efficiently. Note that mixing independent
423  and dependent variables would require assignments, such as
424
425  .. code-block:: python
426
427        StoreEquivalence('x',('y',))
428        StoreEquivalence('y',('z',))
429
430  would require that equivalences be applied in a particular order and
431  thus is implemented as a constraint equation rather than an equivalence.
432       
433  Use StoreEquivalence before calling GenerateConstraints or CheckConstraints
434
435:func:`CheckConstraints`
436   check that input in internally consistent
437
438:func:`GenerateConstraints`
439   generate the internally used tables from constraints and equivalences
440
441:func:`EvaluateMultipliers`
442   Convert any string-specified (formula-based) multipliers to numbers. Call this before
443   using :func:`CheckConstraints` or :func:`GenerateConstraints`.
444   At present, the code may pass only the dict for phase (atom/cell) parameters,
445   but this could be expanded if needed.
446
447:func:`Map2Dict`
448   To determine values for the parameters created in this module, one
449   calls Map2Dict. This will not apply contraints.
450
451:func:`Dict2Map`
452   To take values from the new independent parameters and constraints,
453   one calls Dict2Map and set the parameter array, thus appling contraints.
454
455:func:`Dict2Deriv`
456   Use Dict2Deriv to determine derivatives on independent parameters
457   from those on dependent ones.
458
459:func:`ComputeDepESD`     
460   Use ComputeDepESD to compute uncertainties on dependent variables.
461
462:func:`VarRemapShow`
463   To show a summary of the parameter remapping, one calls VarRemapShow.
464
465*Global Variables*
466------------------
467
468dependentParmList:
469   a list containing group of lists of
470   parameters used in the group. Note that parameters listed in
471   dependentParmList should not be refined as they will not affect
472   the model
473
474indParmList:
475     a list containing groups of Independent parameters defined in
476     each group. This contains both parameters used in parameter
477     redefinitions as well as names of generated new parameters.
478
479arrayList:
480     a list containing group of relationship matrices to relate
481     parameters in dependentParmList to those in indParmList. Unlikely
482     to be used externally.
483
484invarrayList:
485     a list containing group of relationship matrices to relate
486     parameters in indParmList to those in dependentParmList. Unlikely
487     to be used externally.
488
489fixedVarList:
490     a list of variables that have been 'fixed'
491     by defining them as equal to a constant (::var: = 0). Note that
492     the constant value is ignored at present. These variables are
493     later removed from varyList which prevents them from being refined.
494     Unlikely to be used externally.
495
496fixedDict:
497     a dictionary containing the fixed values corresponding
498     to parameter equations.  The dict key is an ascii string, but the
499     dict value is a float.  Unlikely to be used externally.
500
501symGenList:
502     a list of boolean values that will be True to indicate that a constraint
503     (only equivalences) is generated by symmetry (or Pawley overlap)
504     
505problemVars:
506     a list containing variables that show up in constraints producing errors
507
508
509
510*Routines/variables*
511---------------------
512
513Note that parameter names in GSAS-II are strings of form ``<ph#>:<hst#>:<nam>`` or ``<ph#>::<nam>:<at#>``.
514"""
515
516from __future__ import division, print_function
517import numpy as np
518import sys
519import GSASIIpath
520GSASIIpath.SetVersionNumber("$Revision: 3713 $")
521# data used for constraints;
522debug = False # turns on printing as constraint input is processed
523
524# note that constraints are stored listed by contraint groups,
525# where each constraint
526# group contains those parameters that must be handled together
527dependentParmList = []
528'''a list of lists where each item contains a list of parameters in each constraint group.
529note that parameters listed in dependentParmList should not be refined directly.'''
530indParmList = [] # a list of names for the new parameters
531'''a list of lists where each item contains a list for each constraint group with
532fixed values for constraint equations and names of generated (New Var) parameters.
533'''
534arrayList = []
535'''a list of of relationship matrices that map model parameters in each
536constraint group (in :data:`dependentParmList`) to
537generated (New Var) parameters.
538'''
539invarrayList = []
540'''a list of of inverse-relationship matrices that map constrained values and
541generated (New Var) parameters (in :data:`indParmList`) to model parameters
542(in :data:`dependentParmList`).
543'''
544fixedDict = {}
545'''A dict lookup-table containing the fixed values corresponding
546to defined parameter equations. Note the key is the original ascii string
547and the value in the dict is a float.
548'''
549fixedVarList = []
550'''List of variables that should not be refined.
551'''
552symGenList = []
553'''A list of flags that if True indicates a constraint was generated by symmetry
554'''
555problemVars = []
556'''a list of variables causing errors
557'''
558dependentVars = []
559'A list of dependent variables, taken from (:data:`dependentParmList`).'
560independentVars = []
561'A list of dependent variables, taken from (:data:`indParmList`).'
562genVarLookup = {}
563'provides a list of variables that are related to each generated variable'
564paramPrefix = "::constr"
565'A prefix for generated parameter names'
566consNum = 0
567'The number to be assigned to the next constraint to be created'
568
569class ConstraintException(Exception):
570    '''Defines an Exception that is used when an exception is raised processing constraints
571    '''
572    pass
573
574def InitVars():
575    '''Initializes all constraint information'''
576    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,consNum,symGenList
577    dependentParmList = [] # contains a list of parameters in each group
578    arrayList = [] # a list of of relationship matrices
579    invarrayList = [] # a list of inverse relationship matrices
580    indParmList = [] # a list of names for the new parameters
581    fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations
582    symGenList = [] # Flag if constraint is generated by symmetry
583    consNum = 0 # number of the next constraint to be created
584    global genVarLookup
585    genVarLookup = {}
586
587def VarKeys(constr):
588    """Finds the keys in a constraint that represent variables
589    e.g. eliminates any that start with '_'
590
591    :param dict constr: a single constraint entry of form::
592
593        {'var1': mult1, 'var2': mult2,... '_notVar': val,...}
594
595        (see :func:`GroupConstraints`)
596    :returns: a list of keys where any keys beginning with '_' are
597      removed.
598    """
599    return [i for i in constr.keys() if not i.startswith('_')]
600
601
602def GroupConstraints(constrDict):
603    """divide the constraints into groups that share no parameters.
604
605    :param dict constrDict: a list of dicts defining relationships/constraints
606
607    ::
608   
609       constrDict = [{<constr1>}, {<constr2>}, ...]
610
611    where {<constr1>} is {'var1': mult1, 'var2': mult2,... }
612
613    :returns: two lists of lists:
614   
615      * a list of grouped contraints where each constraint grouped containts a list
616        of indices for constraint constrDict entries
617      * a list containing lists of parameter names contained in each group
618     
619      """
620    assignedlist = [] # relationships that have been used
621    groups = [] # contains a list of grouplists
622    ParmList = []
623    for i,consi in enumerate(constrDict):
624        if i in assignedlist: continue # already in a group, skip
625        # starting a new group
626        grouplist = [i,]
627        assignedlist.append(i)
628        groupset = set(VarKeys(consi))
629        changes = True # always loop at least once
630        while(changes): # loop until we can't find anything to add to the current group
631            changes = False # but don't loop again unless we find something
632            for j,consj in enumerate(constrDict):
633                if j in assignedlist: continue # already in a group, skip
634                if len(set(VarKeys(consj)) & groupset) > 0: # true if this needs to be added
635                    changes = True
636                    grouplist.append(j)
637                    assignedlist.append(j)
638                    groupset = groupset | set(VarKeys(consj))
639        group = sorted(grouplist)
640        varlist = sorted(list(groupset))
641        groups.append(group)
642        ParmList.append(varlist)
643    return groups,ParmList
644
645def CheckConstraints(varyList,constrDict,fixedList):
646    '''Takes a list of relationship entries comprising a group of
647    constraints and checks for inconsistencies such as conflicts in
648    parameter/variable definitions and or inconsistently varied parameters.
649
650    :param list varyList: a list of parameters names that will be varied
651
652    :param dict constrDict: a list of dicts defining relationships/constraints
653      (as created in :func:`GSASIIstrIO.ProcessConstraints` and
654      documented in :func:`GroupConstraints`)
655
656    :param list fixedList: a list of values specifying a fixed value for each
657      dict in constrDict. Values are either strings that can be converted to
658      floats or ``None`` if the constraint defines a new parameter rather
659      than a constant.
660
661    :returns: two strings:
662
663      * the first lists conflicts internal to the specified constraints
664      * the second lists conflicts where the varyList specifies some
665        parameters in a constraint, but not all
666       
667      If there are no errors, both strings will be empty
668    '''
669    import re
670    global dependentParmList,arrayList,invarrayList,indParmList,consNum
671    global problemVars
672    # Process the equivalences
673    #    If there are conflicting parameters, move them into constraints. This
674    #    may create new conflicts, requiring movement of other parameters, so
675    #    loop until there are no more changes to make.
676    parmsChanged = True
677    while parmsChanged:
678        parmsChanged = 0
679        errmsg,warnmsg,fixVlist,dropVarList,translateTable = CheckEquivalences(
680            constrDict,varyList)
681        #print('debug: using MoveConfEquiv to address =',errmsg)
682        if problemVars: parmsChanged = MoveConfEquiv(constrDict,fixedList)
683#    GSASIIpath.IPyBreak()
684
685    groups,parmlist = GroupConstraints(constrDict)
686    # scan through parameters in each relationship. Are all varied? If only some are
687    # varied, create a warning message.
688    for group,varlist in zip(groups,parmlist):
689        if len(varlist) == 1: continue
690        for rel in group:
691            varied = 0
692            notvaried = ''
693            for var in constrDict[rel]:
694                if var.startswith('_'): continue
695                if not re.match('[0-9]*:[0-9\*]*:',var):
696                    warnmsg += "\nVariable "+str(var)+" does not begin with a ':'"
697                if var in varyList:
698                    varied += 1
699                else:
700                    if notvaried: notvaried += ', '
701                    notvaried += var
702                if var in fixVlist:
703                    errmsg += '\nParameter '+var+" is Fixed and used in a constraint:\n\t"
704                    if var not in problemVars: problemVars.append(var)
705                    errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
706            if varied > 0 and varied != len(VarKeys(constrDict[rel])):
707                warnmsg += "\nNot all variables refined in constraint:\n\t"
708                warnmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
709                warnmsg += '\nNot refined: ' + notvaried + '\n'
710    if errmsg or warnmsg:
711        return errmsg,warnmsg
712
713    # now look for process each group and create the relations that are needed to form
714    # non-singular square matrix
715    for group,varlist in zip(groups,parmlist):
716        if len(varlist) == 1: continue # a constraint group with a single variable can be ignored
717        if len(varlist) < len(group): # too many relationships -- no can do
718            errmsg += "\nOver-constrained input. "
719            errmsg += "There are more constraints " + str(len(group))
720            errmsg += "\n\tthan variables " + str(len(varlist)) + "\n"
721            for rel in group:
722                errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
723                errmsg += "\n"
724                continue
725        try:
726            multarr = _FillArray(group,constrDict,varlist)
727            _RowEchelon(len(group),multarr,varlist)
728        except:
729            errmsg += "\nSingular input. "
730            errmsg += "There are internal inconsistencies in these constraints\n"
731            for rel in group:
732                errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
733                errmsg += "\n"
734            continue
735        try:
736            multarr = _FillArray(group,constrDict,varlist,FillDiagonals=True)
737            GramSchmidtOrtho(multarr,len(group))
738        except:
739            errmsg += "\nUnexpected singularity with constraints (in Gram-Schmidt)\n"
740            for rel in group:
741                errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
742                errmsg += "\n"
743            continue
744        mapvar = []
745        group = group[:]
746        # scan through all generated and input variables
747        # Check again for inconsistent variable use
748        # for new variables -- where varied and unvaried parameters get grouped
749        # together. I don't think this can happen when not flagged before, but
750        # it does not hurt to check again.
751        for i in range(len(varlist)):
752            varied = 0
753            notvaried = ''
754            if len(group) > 0:
755                rel = group.pop(0)
756                fixedval = fixedList[rel]
757                for var in VarKeys(constrDict[rel]):
758                    if var in varyList:
759                        varied += 1
760                    else:
761                        if notvaried: notvaried += ', '
762                        notvaried += var
763            else:
764                fixedval = None
765            if fixedval is None:
766                varname = paramPrefix + str(consNum) # assign a name to a variable
767                mapvar.append(varname)
768                consNum += 1
769            else:
770                mapvar.append(fixedval)
771            if varied > 0 and notvaried != '':
772                warnmsg += "\nNot all variables refined in generated constraint"
773                warnmsg += '\nPlease report this unexpected error\n'
774                for rel in group:
775                    warnmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
776                    warnmsg += "\n"
777                warnmsg += '\n\tNot refined: ' + notvaried + '\n'
778        try:
779            np.linalg.inv(multarr)           
780        except:
781            errmsg += "\nSingular input. "
782            errmsg += "The following constraints are not "
783            errmsg += "linearly independent\n\tor do not "
784            errmsg += "allow for generation of a non-singular set\n"
785            errmsg += 'Please report this unexpected error\n'
786            for rel in group:
787                errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
788                errmsg += "\n"
789    _setVarLists([])
790    return errmsg,warnmsg
791
792def GenerateConstraints(varyList,constrDict,fixedList,parmDict=None,SeqHist=None):
793    '''Takes a list of relationship entries comprising a group of
794    constraints and builds the relationship lists and their inverse
795    and stores them in global variables Also checks for internal
796    conflicts or inconsistencies in parameter/variable definitions.
797
798    :param list varyList: a list of parameters names (strings of form
799      ``<ph>:<hst>:<nam>``) that will be varied. Note that this is changed here.
800   
801    :param dict constrDict: a list of dicts defining relationships/constraints
802      (as defined in :func:`GroupConstraints`)
803
804    :param list fixedList: a list of values specifying a fixed value for each
805      dict in constrDict. Values are either strings that can be converted to
806      floats, float values or None if the constraint defines a new parameter.
807     
808    :param dict parmDict: a dict containing all parameters defined in current
809      refinement.
810
811    :param int SeqHist: number of current histogram, when used in a sequential
812      refinement. None (default) otherwise. Wildcard variable names are
813      set to the current histogram, when found if not None.
814    '''
815    global dependentParmList,arrayList,invarrayList,indParmList,consNum
816    global genVarLookup
817    msg = ''
818    shortmsg = ''
819    changed = False
820    # Process the equivalences
821    #    If there are conflicting parameters, move them into constraints. This
822    #    may create new conflicts, requiring movement of other parameters, so
823    #    loop until there are no more changes to make.
824    parmsChanged = True
825    while parmsChanged:
826        parmsChanged = 0
827        errmsg,warnmsg,fixVlist,dropVarList,translateTable = CheckEquivalences(
828            constrDict,varyList,parmDict,SeqHist)
829        if problemVars:
830            parmsChanged = MoveConfEquiv(constrDict,fixedList)
831            changed = True
832    if errmsg:
833        msg = errmsg
834    if warnmsg:
835        if msg: msg += '\n'
836        msg += warnmsg
837
838    # scan through parameters in each relationship. Are all varied? If only some are
839    # varied, create an error message.
840    groups,parmlist = GroupConstraints(constrDict)
841    for group,varlist in zip(groups,parmlist):
842        if len(varlist) == 1: continue
843        for rel in group:
844            varied = 0
845            notvaried = ''
846            unused = 0
847            notused = ''
848            for var in constrDict[rel]:
849                if var.startswith('_'): continue
850                if var.split(':')[1] == '*' and SeqHist is not None:
851                    # convert wildcard var to reference current histogram; save translation in table
852                    sv = var.split(':')
853                    sv[1] = str(SeqHist)
854                    translateTable[var] = ':'.join(sv)
855                    var = translateTable[var]
856                if parmDict is not None and var not in parmDict:
857                    unused += 1
858                    if notvaried: notused += ', '
859                    notused += var
860                if var in varyList:
861                    varied += 1
862                else:
863                    if notvaried: notvaried += ', '
864                    notvaried += var
865                if var in fixedVarList:
866                    msg += '\nError: parameter '+var+" is Fixed and used in a constraint:\n\t"
867                    msg += _FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
868            #if unused > 0:# and unused != len(VarKeys(constrDict[rel])):
869            if unused > 0 and unused != len(VarKeys(constrDict[rel])):
870                #msg += "\nSome (but not all) variables in constraint are not defined:\n\t"
871                #msg += _FormatConstraint(constrDict[rel],fixedList[rel])
872                #msg += '\nNot used: ' + notused + '\n'
873                shortmsg += notused+" not used in constraint "+_FormatConstraint(constrDict[rel],fixedList[rel])
874            elif varied > 0 and varied != len(VarKeys(constrDict[rel])):
875                #msg += "\nNot all variables refined in constraint:\n\t"
876                #msg += _FormatConstraint(constrDict[rel],fixedList[rel])
877                #msg += '\nNot refined: ' + notvaried + '\n'
878                shortmsg += notvaried+" not varied in constraint "+_FormatConstraint(constrDict[rel],fixedList[rel])
879    # if there were errors found, go no farther
880    if shortmsg and SeqHist is not None:
881        if msg:
882            print (' *** ERROR in constraint definitions! ***')
883            print (msg)
884            raise ConstraintException
885        print ('*** Sequential refinement: ignoring constraint definition(s): ***')
886        print (shortmsg)
887        msg = ''
888    elif shortmsg:
889        msg += shortmsg
890    if msg:
891        print (' *** ERROR in constraint definitions! ***')
892        print (msg)
893        raise ConstraintException
894               
895    # now process each group and create the relations that are needed to form
896    # a non-singular square matrix
897    # If all are varied and this is a constraint equation, then set VaryFree flag
898    # so that the newly created relationships will be varied
899    for group,varlist in zip(groups,parmlist):
900        if len(varlist) == 1: continue
901        # for constraints, if all included variables are refined,
902        # set the VaryFree flag, and remaining degrees of freedom will be
903        # varied (since consistency was checked, if any one variable is
904        # refined, then assume that all are)
905        varsList = [] # make a list of all the referenced variables as well
906        VaryFree = False
907        for rel in group:
908            varied = 0
909            unused = 0
910            for var in VarKeys(constrDict[rel]):
911                var = translateTable.get(var,var) # replace wildcards
912                if parmDict is not None and var not in parmDict:
913                    unused += 1                   
914                if var not in varsList: varsList.append(var)
915                if var in varyList: varied += 1
916            if fixedList[rel] is not None and varied > 0:
917                VaryFree = True
918        if len(varlist) < len(group): # too many relationships -- no can do
919            msg = 'too many relationships'
920            break
921        # Since we checked before, if any variables are unused, then all must be.
922        # If so, this set of relationships can be ignored
923        if unused:
924            if debug: print('Constraint ignored (all variables undefined)')
925            if debug: print ('    '+_FormatConstraint(constrDict[rel],fixedList[rel]))
926            continue
927        # fill in additional degrees of freedom
928        try:
929            arr = _FillArray(group,constrDict,varlist)
930            _RowEchelon(len(group),arr,varlist)
931            constrArr = _FillArray(group,constrDict,varlist,FillDiagonals=True)
932            GramSchmidtOrtho(constrArr,len(group))
933        except:
934            msg = 'Singular relationships'
935            break
936        mapvar = []
937        group = group[:]
938        # scan through all generated and input relationships, we need to add to the varied list
939        # all the new parameters where VaryFree has been set or where a New Var is varied.
940        #
941        # If a group does not contain any fixed values (constraint equations)
942        # and nothing in the group is varied, drop this group, so that the
943        # dependent parameters can be refined individually.
944        unused = True
945        for i in range(len(varlist)):
946            if len(group) > 0: # get the original equation reference
947                rel = group.pop(0)
948                fixedval = fixedList[rel]
949                varyflag = constrDict[rel].get('_vary',False)
950                varname = constrDict[rel].get('_name','')
951            else: # this relationship has been generated
952                varyflag = False
953                varname = ''
954                fixedval = None
955            if fixedval is None: # this is a new variable, not a constraint
956                if not varname:
957                    varname = paramPrefix + str(consNum) # no assigned name, create one
958                    consNum += 1
959                mapvar.append(varname)
960                genVarLookup[varname] = varlist # save list of variables related to this new var
961                # vary the new relationship if it is a degree of freedom in
962                # a set of contraint equations or if a New Var is flagged to be varied.
963                if VaryFree or varyflag: 
964                    unused = False
965                    varyList.append(varname)
966                    # fix (prevent varying) of all the variables inside the constraint group
967                    # (dependent vars)
968                    for var in varsList:
969                        if var in varyList: varyList.remove(var)
970            else:
971                unused = False
972                mapvar.append(fixedval)
973        if unused: # skip over constraints that don't matter (w/o fixed value or any refined variables)
974            if debug: print('Constraint ignored (all variables unrefined)')
975            if debug: print ('   '+_FormatConstraint(constrDict[rel],fixedList[rel]))
976            continue 
977        dependentParmList.append([translateTable.get(var,var) for var in varlist])
978        arrayList.append(constrArr)
979        invarrayList.append(np.linalg.inv(constrArr))
980        indParmList.append(mapvar)
981        symGenList.append(False)
982    if msg:
983        print (' *** ERROR in constraint definitions! ***')
984        print (msg)
985        print (VarRemapShow(varyList))
986        raise ConstraintException
987    # setup dictionary containing the fixed values
988    global fixedDict
989    # key is original ascii string, value is float
990    for fixedval in fixedList:
991        if fixedval:
992            fixedDict[fixedval] = float(fixedval)
993    _setVarLists(dropVarList)
994    if changed:
995        print(60*'=')
996        print('Constraints were reclassified to avoid conflicts, as below:')
997        print (VarRemapShow(varyList,True))
998        print(60*'=')
999   
1000def _setVarLists(dropVarList):
1001    '''Make list of dependent and independent variables (after dropping unused vars in dropVarList)
1002    '''
1003    global dependentParmList,indParmList
1004    global dependentVars
1005    global independentVars
1006    dependentVars = []
1007    independentVars = []
1008    for varlist,mapvars in zip(dependentParmList,indParmList):  # process all constraints
1009        for mv in mapvars:
1010            if mv in dropVarList: continue
1011            if mv not in independentVars: independentVars.append(mv)
1012        for mv in varlist:
1013            if mv in dropVarList: continue
1014            if mv not in dependentVars: dependentVars.append(mv)
1015    if debug: # on debug, show what is parsed & generated, semi-readable
1016        print (50*'-')
1017        #print (VarRemapShow(varyList))
1018        #print ('Varied: ',varyList)
1019        print ('Not Varied: ',fixedVarList)
1020
1021# def CheckEquivalences(constrDict,varyList):
1022#     global dependentParmList,arrayList,invarrayList,indParmList,consNum
1023#     global problemVars
1024#     warnmsg = ''
1025#     errmsg = ''
1026#     problemVars = []
1027#     # process fixed variables (holds)
1028#     fixVlist = [] # list of Fixed vars
1029#     constrVars = [] # list of vars in constraint expressions
1030#     for cdict in constrDict:
1031#         # N.B. No "_" names in holds
1032#         if len(cdict) == 1:
1033#             fixVlist.append(list(cdict.keys())[0])
1034#         else:
1035#             constrVars += cdict.keys() # this will include _vary (not a problem)
1036#     # process equivalences: make a list of dependent and independent vars
1037#     #    and check for repeated uses (repetition of a parameter as an
1038#     #    independent var is OK)
1039#     indepVarList = []
1040#     depVarList = []
1041#     multdepVarList = []
1042#     for varlist,mapvars,multarr,invmultarr in zip(
1043#         dependentParmList,indParmList,arrayList,invarrayList):
1044#         if multarr is None: # an equivalence
1045#             zeromult = False
1046#             for mv in mapvars:
1047#                 varied = 0
1048#                 notvaried = ''
1049#                 if mv in varyList:
1050#                     varied += 1
1051#                 else:
1052#                     if notvaried: notvaried += ', '
1053#                     notvaried += mv
1054#                 if mv not in indepVarList: indepVarList.append(mv)
1055#                 for v,m in zip(varlist,invmultarr):
1056#                     if v in indepVarList:
1057#                         errmsg += '\nVariable '+v+' is used to set values in a constraint before its value is set in another constraint\n'
1058#                         if v not in problemVars: problemVars.append(v)
1059#                     if m == 0: zeromult = True
1060#                     if v in varyList:
1061#                         varied += 1
1062#                     else:
1063#                         if notvaried: notvaried += ', '
1064#                         notvaried += v
1065#                     if v in depVarList:
1066#                         multdepVarList.append(v)
1067#                     else:
1068#                         depVarList.append(v)
1069#             if varied > 0 and varied != len(varlist)+1:
1070#                 warnmsg += "\nNot all variables refined in equivalence:\n\t"
1071#                 s = ""
1072#                 for v in varlist:
1073#                     if s != "": s+= " & "
1074#                     s += str(v)           
1075#                 warnmsg += str(mv) + " => " + s
1076#                 warnmsg += '\nNot refined: ' + notvaried + '\n'
1077#             if zeromult:
1078#                 errmsg += "\nZero multiplier is invalid in equivalence:\n\t"
1079#                 s = ""
1080#                 for v in varlist:
1081#                     if s != "": s+= " & "
1082#                     s += str(v)           
1083#                 errmsg += str(mv) + " => " + s + '\n'
1084#     # check for errors:
1085#     if len(multdepVarList) > 0:
1086#         errmsg += "\nThe following parameters(s) are used in conflicting Equivalence relations as dependent variables:\n"
1087#         s = ''
1088#         for var in sorted(set(multdepVarList)):
1089#             if v not in problemVars: problemVars.append(v)
1090#             if s != "": s+= ", "
1091#             s += str(var)           
1092#         errmsg += '\t'+ s + '\n'
1093#     equivVarList = list(set(indepVarList).union(set(depVarList)))
1094#     if debug: print ('equivVarList',equivVarList)
1095#     # check for parameters that are both fixed and in an equivalence (not likely)
1096#     inboth = set(fixVlist).intersection(set(equivVarList))
1097#     if len(inboth) > 0:
1098#         errmsg += "\nThe following parameter(s) are used in both Equivalence and Fixed constraints:\n"
1099#         s = ''
1100#         for var in sorted(inboth):
1101#             if var not in problemVars: problemVars.append(var)
1102#             if s != "": s+= ", "
1103#             s += str(var)
1104#         errmsg += '\t'+ s + '\n'
1105#     # check for parameters that in both an equivalence and a constraint expression
1106#     inboth = set(constrVars).intersection(set(equivVarList))
1107#     if len(inboth) > 0:
1108#         errmsg += "\nThe following parameter(s) are used in both Equivalence and Equiv or new var constraints:\n"
1109#         s = ''
1110#         for var in sorted(inboth):
1111#             if var not in problemVars: problemVars.append(var)
1112#             if s != "": s+= ", "
1113#             s += str(var)
1114#         errmsg += '\t'+ s + '\n'
1115#     return errmsg,warnmsg,fixVlist
1116
1117def CheckEquivalences(constrDict,varyList,parmDict=None,SeqHist=None):
1118    '''Process equivalence constraints, looking for conflicts such as
1119    where a parameter is used in both an equivalence and a constraint expression
1120    or where chaining is done (A->B and B->C).
1121    When called during refinements, parmDict is defined, and for sequential refinement
1122    SeqHist ia also defined.
1123
1124      * parmDict is used to remove equivalences where a parameter is not present
1125        in a refinement
1126      * SeqHist is used to rename wild-card parameter names in sequential
1127        refinements to use the current histogram.
1128    '''
1129    global dependentParmList,arrayList,invarrayList,indParmList,consNum
1130    global problemVars
1131    warnmsg = ''
1132    errmsg = ''
1133    problemVars = []
1134    # process fixed variables (holds)
1135    fixVlist = [] # list of Fixed vars
1136    constrVars = [] # list of vars in constraint expressions
1137    for cdict in constrDict:
1138        # N.B. No "_" names in holds
1139        if len(cdict) == 1:
1140            fixVlist.append(list(cdict.keys())[0])
1141        else:
1142            constrVars += cdict.keys() # this will include _vary (not a problem)
1143    # process equivalences: make a list of dependent and independent vars
1144    #    and check for repeated uses (repetition of a parameter as an
1145    #    independent var is OK)
1146    indepVarList = []
1147    depVarList = []
1148    multdepVarList = []
1149    dropVarList = []
1150    translateTable = {} # lookup table for wildcard referenced variables
1151    for varlist,mapvars,multarr,invmultarr in zip(
1152        dependentParmList,indParmList,arrayList,invarrayList):
1153        if multarr is None: # an equivalence
1154            zeromult = False
1155            for i,mv in enumerate(mapvars):
1156                if mv.split(':')[1] == '*' and SeqHist is not None:
1157                    # convert wildcard var to reference current histogram; save translation in table
1158                    sv = mv.split(':')
1159                    sv[1] = str(SeqHist)
1160                    mv = translateTable[mv] = ':'.join(sv)
1161                    mapvars[i] = mv
1162                varied = 0
1163                notvaried = ''
1164                if mv in varyList:
1165                    varied += 1
1166                else:
1167                    if notvaried: notvaried += ', '
1168                    notvaried += mv
1169                if parmDict is not None and mv not in parmDict:
1170                    print ("Dropping equivalence for variable "+str(mv)+". Not defined in this refinement")
1171                    if mv not in dropVarList: dropVarList.append(mv)
1172                if mv not in indepVarList: indepVarList.append(mv)
1173            for i,(v,m) in enumerate(zip(varlist,invmultarr)):
1174                if v.split(':')[1] == '*' and SeqHist is not None:
1175                    # convert wildcard var to reference current histogram; save translation in table
1176                    sv = v.split(':')
1177                    sv[1] = str(SeqHist)
1178                    varlist[i] = v = translateTable[v] = ':'.join(sv)
1179                if parmDict is not None and v not in parmDict:
1180                    print ("Dropping equivalence for dep. variable "+str(v)+". Not defined in this refinement")
1181                    if v not in dropVarList: dropVarList.append(v)
1182                    continue
1183                if m == 0: zeromult = True
1184                if v in varyList:
1185                    varied += 1
1186                else:
1187                    if notvaried: notvaried += ', '
1188                    notvaried += v
1189                if v in indepVarList:
1190                    errmsg += '\nVariable '+v+' is used to set values in a constraint before its value is set in another constraint\n'
1191                    if v not in problemVars: problemVars.append(v)
1192                if v in depVarList:
1193                    multdepVarList.append(v)
1194                else:
1195                    depVarList.append(v)
1196            if varied > 0 and varied != len(varlist)+1:
1197                warnmsg += "\nNot all variables refined in equivalence:\n\t"
1198                s = ""
1199                for v in varlist:
1200                    if s != "": s+= " & "
1201                    s += str(v)           
1202                warnmsg += str(mv) + " => " + s
1203                warnmsg += '\nNot refined: ' + notvaried + '\n'
1204            if zeromult:
1205                errmsg += "\nZero multiplier is invalid in equivalence:\n\t"
1206                s = ""
1207                for v in varlist:
1208                    if s != "": s+= " & "
1209                    s += str(v)           
1210                errmsg += str(mv) + " => " + s + '\n'
1211    # check for errors:
1212    if len(multdepVarList) > 0:
1213        errmsg += "\nThe following parameters(s) are used in conflicting Equivalence relations as dependent variables:\n"
1214        s = ''
1215        for var in sorted(set(multdepVarList)):
1216            if v not in problemVars: problemVars.append(v)
1217            if s != "": s+= ", "
1218            s += str(var)           
1219        errmsg += '\t'+ s + '\n'
1220    equivVarList = list(set(indepVarList).union(set(depVarList)))
1221    if debug: print ('equivVarList',equivVarList)
1222    # check for parameters that are both fixed and in an equivalence (not likely)
1223    inboth = set(fixVlist).intersection(set(equivVarList))
1224    if len(inboth) > 0:
1225        errmsg += "\nThe following parameter(s) are used in both Equivalence and Fixed constraints:\n"
1226        s = ''
1227        for var in sorted(inboth):
1228            if var not in problemVars: problemVars.append(var)
1229            if s != "": s+= ", "
1230            s += str(var)
1231        errmsg += '\t'+ s + '\n'
1232    # check for parameters that in both an equivalence and a constraint expression
1233    inboth = set(constrVars).intersection(set(equivVarList))
1234    if len(inboth) > 0:
1235        errmsg += "\nThe following parameter(s) are used in both Equivalence and Equiv or new var constraints:\n"
1236        s = ''
1237        for var in sorted(inboth):
1238            if var not in problemVars: problemVars.append(var)
1239            if s != "": s+= ", "
1240            s += str(var)
1241        errmsg += '\t'+ s + '\n'
1242    return errmsg,warnmsg,fixVlist,dropVarList,translateTable
1243
1244def MoveConfEquiv(constrDict,fixedList):
1245    '''Address conflicts in Equivalence constraints by creating an constraint equation
1246    that has the same action as the equivalence and removing the Equivalence
1247    '''
1248    global dependentParmList,arrayList,invarrayList,indParmList,consNum
1249    global problemVars
1250    parmsChanged = 0
1251    for i,(varlist,mapvars) in enumerate(zip(dependentParmList,indParmList)):
1252        conf = False
1253        for mv in mapvars:
1254            if mv in problemVars:
1255                conf = True
1256                break
1257        for v in varlist:
1258            if v in problemVars:
1259                conf = True
1260                break
1261        if conf:
1262            parmsChanged += 1
1263            indvar = indParmList[i][0]
1264            for dep,mult in zip(dependentParmList[i],invarrayList[i]):
1265                #print('debug replacing equiv with constraint equation 0=',{indvar:-1.,dep:mult[0]})
1266                constrDict += [{indvar:-1.,dep:mult[0]}]
1267                fixedList += ['0.0']
1268            dependentParmList[i] = None
1269    if parmsChanged:
1270        for i in range(len(dependentParmList)-1,-1,-1):
1271            if dependentParmList[i] is None:
1272                del dependentParmList[i],indParmList[i],arrayList[i],invarrayList[i]
1273    return parmsChanged
1274
1275def StoreEquivalence(independentVar,dependentList,symGen=True):
1276    '''Takes a list of dependent parameter(s) and stores their
1277    relationship to a single independent parameter (independentVar)
1278
1279    :param str independentVar: name of master parameter that will be used to determine the value
1280      to set the dependent variables
1281
1282    :param list dependentList: a list of parameters that will set from
1283         independentVar. Each item in the list can be a string with the parameter
1284         name or a tuple containing a name and multiplier:
1285         ``['::parm1',('::parm2',.5),]``
1286
1287    '''
1288   
1289    global dependentParmList,arrayList,invarrayList,indParmList,symGenList
1290    mapList = []
1291    multlist = []
1292    allfloat = True
1293    for var in dependentList:
1294        if isinstance(var, str):
1295            mult = 1.0
1296        elif len(var) == 2:
1297            var,mult = var
1298        else:
1299            raise Exception("Cannot parse "+repr(var) + " as var or (var,multiplier)")
1300        mapList.append(var)
1301        try:           
1302            multlist.append(tuple((float(mult),)))
1303        except:
1304            allfloat = False
1305            multlist.append(tuple((mult,)))
1306    # added relationships to stored values
1307    arrayList.append(None)
1308    if allfloat:
1309        invarrayList.append(np.array(multlist))
1310    else:
1311        invarrayList.append(multlist)
1312    indParmList.append(list((independentVar,)))
1313    dependentParmList.append(mapList)
1314    symGenList.append(symGen)
1315    return
1316
1317def EvaluateMultipliers(constList,*dicts):
1318    '''Convert multipliers for constraints and equivalences that are specified
1319    as strings into values. The strings can specify values in the parameter dicts as
1320    well as normal Python functions, such as "2*np.cos(0::Ax:2/2.)"
1321   
1322    :param list constList: a list of dicts containing constraint expressions
1323    :param \*dict1: one or more dicts containing GSAS-II parameters and their values
1324       can be specified
1325    :returns: an empty string if there were no errors, or an error message listing
1326       the strings that could not be converted.
1327    '''
1328    def SubfromParmDict(s,prmDict):
1329        for key in prmDict:
1330            if key in s:
1331                s = s.replace(key,str(prmDict[key]))
1332        return eval(s)
1333    prmDict = {}
1334    for d in dicts: prmDict.update(d) # combine all passed parameter dicts
1335    problemList = ""
1336    # loop through multipliers in contraint expressions
1337    for const in constList:
1338        for key in const:
1339            try:
1340                1+const[key]
1341                continue
1342            except:
1343                pass
1344            try:
1345                newval = SubfromParmDict(const[key][:],prmDict)
1346                if GSASIIpath.GetConfigValue('debug'):
1347                    print('Changing ',const[key],'to',newval)
1348                const[key] = newval               
1349            except:
1350                if problemList: problemList += ", "
1351                problemList += const[key]
1352    # loop through multipliers in equivalences
1353    global arrayList,invarrayList
1354    for i,(a,valList) in enumerate(zip(arrayList,invarrayList)):
1355        if a is not None: continue # ignore if not equiv
1356        try:
1357            valList.shape
1358            continue # ignore if already a numpy array
1359        except:
1360            pass
1361        repList = []
1362        for v in valList:
1363            try:
1364                1+v[0]
1365                repList.append(tuple((v[0],)))
1366                continue
1367            except:
1368                pass
1369            try:
1370                newval = SubfromParmDict(v[0][:],prmDict)
1371                if GSASIIpath.GetConfigValue('debug'):
1372                    print('Changing ',v[0],'to',newval)
1373                repList.append(tuple((newval,)))
1374            except:
1375                if problemList: problemList += ", "
1376                problemList += v[0]
1377                repList.append(tuple(('error',)))
1378        invarrayList[i] = np.array(repList)
1379    return problemList
1380
1381def GetDependentVars():
1382    '''Return a list of dependent variables: e.g. variables that are
1383    constrained in terms of other variables
1384
1385    :returns: a list of variable names
1386
1387    '''
1388    global dependentVars
1389    return dependentVars
1390
1391def GetIndependentVars():
1392    '''Return a list of independent variables: e.g. variables that are
1393    created by constraints of other variables
1394
1395    :returns: a list of variable names
1396
1397    '''
1398    global independentVars
1399    return independentVars
1400
1401def PrintIndependentVars(parmDict,varyList,sigDict,PrintAll=False,pFile=None):
1402    '''Print the values and uncertainties on the independent variables'''
1403    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
1404    printlist = []
1405    mapvars = GetIndependentVars()
1406    for i,name in enumerate(mapvars):
1407        if name in fixedDict: continue
1408        if PrintAll or name in varyList:
1409            sig = sigDict.get(name)
1410            printlist.append([name,parmDict[name],sig])
1411    if len(printlist) == 0: return
1412    s1 = ''
1413    s2 = ''
1414    s3 = ''
1415    pFile.write(130*'-'+'\n')
1416    pFile.write("Variables generated by constraints\n")
1417    printlist.append(3*[None])
1418    for name,val,esd in printlist:
1419        if len(s1) > 120 or name is None:
1420            pFile.write(''+'\n')
1421            pFile.write(s1+'\n')
1422            pFile.write(s2+'\n')
1423            pFile.write(s3+'\n')
1424            s1 = ''
1425            if name is None: break
1426        if s1 == "":
1427            s1 = ' name  :'
1428            s2 = ' value :'
1429            s3 = ' sig   :'
1430        s1 += '%15s' % (name)
1431        s2 += '%15.4f' % (val)
1432        if esd is None:
1433            s3 += '%15s' % ('n/a')
1434        else:   
1435            s3 += '%15.4f' % (esd)
1436
1437def ComputeDepESD(covMatrix,varyList,parmDict):
1438    '''Compute uncertainties for dependent parameters from independent ones
1439    returns a dictionary containing the esd values for dependent parameters
1440    '''
1441    sigmaDict = {}
1442    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
1443        #if invmultarr is None: continue # probably not needed
1444#        try:
1445#            valuelist = [parmDict[var] for var in mapvars]
1446#        except KeyError:
1447#            continue
1448        # get the v-covar matrix for independent parameters
1449        vcov = np.zeros((len(mapvars),len(mapvars)))
1450        for i1,name1 in enumerate(mapvars):
1451            if name1 not in varyList: continue
1452            iv1 = varyList.index(name1)
1453            for i2,name2 in enumerate(mapvars):
1454                if name2 not in varyList: continue
1455                iv2 = varyList.index(name2)
1456                vcov[i1][i2] = covMatrix[iv1][iv2]
1457        # vec is the vector that multiplies each of the independent values
1458        for v,vec in zip(varlist,invmultarr):
1459            sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec)))
1460    return sigmaDict
1461
1462def _FormatConstraint(RelDict,RelVal):
1463    '''Formats a Constraint or Function for use in a convenient way'''
1464    linelen = 45
1465    s = [""]
1466    for var,val in RelDict.items():
1467        if var.startswith('_'): continue
1468        if len(s[-1]) > linelen: s.append(' ')
1469        m = val
1470        if s[-1] != "" and m >= 0:
1471            s[-1] += ' + '
1472        elif s[-1] != "":
1473            s[-1] += ' - '
1474            m = abs(m)
1475        s[-1] += '%.3f*%s '%(m,var)
1476    if len(s[-1]) > linelen: s.append(' ')
1477    if RelVal is None:
1478        s[-1] += ' = New variable'
1479    else:
1480        s[-1] += ' = ' + RelVal
1481    s1 = ''
1482    for s2 in s:
1483        if s1 != '': s1 += '\n\t'
1484        s1 += s2
1485    return s1
1486
1487def VarRemapShow(varyList,inputOnly=False):
1488    '''List out the saved relationships. This should be done after the constraints have been
1489    defined using :func:`StoreEquivalence`, :func:`GroupConstraints` and :func:`GenerateConstraints`.
1490
1491    :returns: a string containing the details of the contraint relationships
1492    '''
1493    s = ''
1494    if len(fixedVarList) > 0:
1495        s += 'Fixed Variables:\n'
1496        for v in fixedVarList:
1497            s += '    ' + v + '\n'
1498    if not inputOnly:
1499        s += 'User-supplied variable mapping relations:\n'
1500    symout = ''
1501    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,symGenList
1502
1503    for varlist,mapvars,multarr,invmultarr,symFlag in zip(
1504        dependentParmList,indParmList,arrayList,invarrayList,symGenList):
1505        for i,mv in enumerate(mapvars):
1506            if multarr is None:
1507#                s1 = '  ' + str(mv) + ' is equivalent to parameter(s): '
1508                if len(varlist) == 1:
1509                    s1 = '   ' + str(mv) + ' is equivalent to '
1510                else:
1511                    s1 = '   ' + str(mv) + ' is equivalent to parameters: '
1512                j = 0
1513                for v,m in zip(varlist,invmultarr):
1514                    if debug: print ('v,m[0]: ',v,m[0])
1515                    if len(s1.split('\n')[-1]) > 75: s1 += '\n        '
1516                    if j > 0: s1 += ' & '
1517                    j += 1
1518                    s1 += str(v)
1519                    if m != 1:
1520                        s1 += " / " + str(m[0])
1521                if symFlag:
1522                    symout += s1 + '\n'
1523                else:
1524                    s += s1 + '\n'
1525                continue
1526            s += %s = ' % mv
1527            j = 0
1528            for m,v in zip(multarr[i,:],varlist):
1529                if m == 0: continue
1530                if j > 0: s += ' + '
1531                j += 1
1532                s += '(%s * %s)' % (m,v)
1533            if mv in varyList: s += ' VARY'
1534            s += '\n'
1535    if symout:
1536        s += 'Symmetry-generated relations:\n' + symout
1537    if inputOnly: return s
1538    s += 'Inverse variable mapping relations:\n'
1539    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
1540        for i,mv in enumerate(varlist):
1541            s += %s = ' % mv
1542            j = 0
1543            for m,v in zip(invmultarr[i,:],mapvars):
1544                if m == 0: continue
1545                if j > 0: s += ' + '
1546                j += 1
1547                s += '(%s * %s)' % (m,v)
1548            s += '\n'
1549    return s
1550
1551def GetSymEquiv():
1552    '''Return the automatically generated (equivalence) relationships.
1553
1554    :returns: a list of strings containing the details of the contraint relationships
1555    '''
1556    symout = []
1557    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,symGenList
1558
1559    for varlist,mapvars,multarr,invmultarr,symFlag in zip(
1560        dependentParmList,indParmList,arrayList,invarrayList,symGenList):
1561        for i,mv in enumerate(mapvars):
1562            if not symFlag: continue
1563            if multarr is None:
1564                #s1 = str(mv) + ' = '
1565                s1 = ''
1566                s2 = ' = ' + str(mv)
1567                j = 0
1568                for v,m in zip(varlist,invmultarr):
1569                    if debug: print ('v,m[0]: ',v,m[0])
1570                    if len(s1.split('\n')[-1]) > 75: s1 += '\n        '
1571                    if j > 0: s1 += ' =  '
1572                    j += 1
1573                    s1 += str(v)
1574                    if m != 1:
1575                        s1 += " / " + str(m[0])
1576                    symout.append(s1+s2)
1577                continue
1578            else:
1579                s = %s = ' % mv
1580                j = 0
1581                for m,v in zip(multarr[i,:],varlist):
1582                    if m == 0: continue
1583                    if j > 0: s += ' + '
1584                    j += 1
1585                    s += '(%s * %s)' % (m,v)
1586                print ('unexpected sym op='+s)
1587    return symout
1588
1589def Dict2Deriv(varyList,derivDict,dMdv):
1590    '''Compute derivatives for Independent Parameters from the
1591    derivatives for the original parameters
1592
1593    :param list varyList: a list of parameters names that will be varied
1594
1595    :param dict derivDict: a dict containing derivatives for parameter values keyed by the
1596      parameter names.
1597
1598    :param list dMdv: a Jacobian, as a list of np.array containing derivatives for dependent
1599      parameter computed from derivDict
1600
1601    '''
1602    global dependentParmList,arrayList,invarrayList,indParmList,invarrayList
1603    for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
1604        for i,name in enumerate(mapvars):
1605            # grouped variables: need to add in the derv. w/r
1606            # dependent variables to the independent ones
1607            if name not in varyList: continue # skip if independent var not varied
1608            if multarr is None:
1609                for v,m in zip(varlist,invmultarr):
1610                    if debug: print ('start dMdv',dMdv[varyList.index(name)])
1611                    if debug: print ('add derv',v,'/',m[0],'to derv',name,'add=',derivDict[v] / m[0])
1612                    if m == 0: continue
1613                    dMdv[varyList.index(name)] += derivDict[v] / m[0]
1614            else:
1615                for v,m in zip(varlist,multarr[i,:]):
1616                    if debug: print ('start dMdv',dMdv[varyList.index(name)])
1617                    if debug: print ('add derv',v,'*',m,'to derv',name,'add=',m * derivDict[v])
1618                    if m == 0: continue
1619                    dMdv[varyList.index(name)] += m * derivDict[v]
1620
1621def Map2Dict(parmDict,varyList):
1622    '''Create (or update) the Independent Parameters from the original
1623    set of Parameters
1624
1625    Removes dependent variables from the varyList
1626
1627    This should be done once, after the constraints have been
1628    defined using :func:`StoreEquivalence`,
1629    :func:`GroupConstraints` and :func:`GenerateConstraints` and
1630    before any variable refinement is done. This completes the parameter
1631    dictionary by defining independent parameters and it satisfies the
1632    constraint equations in the initial parameters
1633
1634    :param dict parmDict: a dict containing parameter values keyed by the
1635      parameter names.
1636      This will contain updated values for both dependent and independent
1637      parameters after Dict2Map is called. It will also contain some
1638      unexpected entries of every constant value {'0':0.0} & {'1.0':1.0},
1639      which do not cause any problems.
1640
1641    :param list varyList: a list of parameters names that will be varied
1642   
1643
1644    '''
1645    # process the Independent vars: remove dependent ones from varylist
1646    # and then compute values for the Independent ones from their dependents
1647    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
1648    for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
1649        for item in varlist:
1650            try:
1651                varyList.remove(item)
1652            except ValueError:
1653                pass
1654        if multarr is None: continue
1655        valuelist = [parmDict[var] for var in varlist]
1656        parmDict.update(zip(mapvars,
1657                            np.dot(multarr,np.array(valuelist)))
1658                        )
1659    # now remove fixed variables from the varyList
1660    global fixedVarList
1661    for item in fixedVarList:
1662        try:
1663            varyList.remove(item)
1664        except ValueError:
1665            pass
1666    # Set constrained parameters to their fixed values
1667    parmDict.update(fixedDict)
1668
1669def Dict2Map(parmDict,varyList):
1670    '''Applies the constraints defined using :func:`StoreEquivalence`,
1671    :func:`GroupConstraints` and :func:`GenerateConstraints` by changing
1672    values in a dict containing the parameters. This should be
1673    done before the parameters are used for any computations
1674
1675    :param dict parmDict: a dict containing parameter values keyed by the
1676      parameter names.
1677      This will contain updated values for both dependent and independent
1678      parameters after Dict2Map is called. It will also contain some
1679      unexpected entries of every constant value {'0':0.0} & {'1.0':1.0},
1680      which do not cause any problems.
1681
1682    :param list varyList: a list of parameters names that will be varied
1683   
1684    '''
1685    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
1686    # reset fixed values (should not be needed, but very quick)
1687    # - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended
1688    # not needed, but also not a problem - BHT
1689    parmDict.update(fixedDict)
1690    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
1691        #if invmultarr is None: continue
1692        try: 
1693            valuelist = [parmDict[var] for var in mapvars]
1694        except KeyError:
1695            continue
1696        parmDict.update(zip(varlist,np.dot(invmultarr,np.array(valuelist))))
1697
1698#======================================================================
1699# internal routines follow (these routines are unlikely to be called
1700# from outside the module)
1701
1702def GramSchmidtOrtho(a,nkeep=0):
1703    '''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to
1704    find orthonormal unit vectors relative to first row.
1705
1706    If nkeep is non-zero, the first nkeep rows in the array are not changed
1707   
1708    input:
1709       arrayin: a 2-D non-singular square array
1710    returns:
1711       a orthonormal set of unit vectors as a square array
1712    '''
1713    def proj(a,b):
1714        'Projection operator'
1715        return a*(np.dot(a,b)/np.dot(a,a))
1716    for j in range(nkeep,len(a)):
1717        for i in range(j):
1718            a[j] -= proj(a[i],a[j])
1719        if np.allclose(np.linalg.norm(a[j]),0.0):
1720            raise ConstraintException("Singular input to GramSchmidtOrtho")
1721        a[j] /= np.linalg.norm(a[j])
1722    return a
1723
1724def _FillArray(sel,dict,collist,FillDiagonals=False):
1725    '''Construct a n by n matrix (n = len(collist)
1726    filling in the rows using the relationships defined
1727    in the dictionaries selected by sel
1728
1729    If FillDiagonals is True, diagonal elements in the
1730    array are set to 1.0
1731    '''
1732    n = len(collist)
1733    if FillDiagonals:
1734        arr = np.eye(n)
1735    else:
1736        arr = np.zeros(2*[n,])
1737    # fill the top rows
1738    for i,cnum in enumerate(sel):
1739        for j,var in enumerate(collist):
1740            arr[i,j] = dict[cnum].get(var,0)
1741    return arr
1742
1743def _SwapColumns(i,m,v):
1744    '''Swap columns in matrix m as well as the labels in v
1745    so that element (i,i) is replaced by the first non-zero element in row i after that element
1746
1747    Throws an exception if there are no non-zero elements in that row
1748    '''
1749    for j in range(i+1,len(v)):
1750        if not np.allclose(m[i,j],0):
1751            m[:,(i,j)] = m[:,(j,i)]
1752            v[i],v[j] = v[j],v[i]
1753            return
1754    else:
1755        raise ConstraintException('Singular input')
1756
1757def _RowEchelon(m,arr,collist):
1758    '''Convert the first m rows in Matrix arr to row-echelon form
1759    exchanging columns in the matrix and collist as needed.
1760
1761    throws an exception if the matrix is singular because
1762    the first m rows are not linearly independent
1763    '''
1764    for i in range(m):
1765        if np.allclose(arr[i,i],0):
1766            _SwapColumns(i,arr,collist)
1767        arr[i,:] /= arr[i,i] # normalize row
1768        # subtract current row from subsequent rows to set values to left of diagonal to 0
1769        for j in range(i+1,m):
1770            arr[j,:] -= arr[i,:] * arr[j,i]
1771
1772if __name__ == "__main__":
1773    parmdict = {}
1774    constrDict = [
1775        {'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},
1776        {'0:0:eA': 0.0},
1777        {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
1778        {'1::C(10,0,1)': 1.0, '2::C(10,0,1)': 1.0},
1779        {'1::AUiso:0': 1.0, '0::AUiso:0': 1.0},
1780        {'0::A0': 0.0}
1781        ]
1782    fixedList = ['5.0', '0', None, None, '1.0', '0']
1783    StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), ))
1784    #StoreEquivalence('1::atomx:3',('2::atomx:3', ('2::atomz:3',2,), )) # error: dependent & independent vars mixed
1785    #StoreEquivalence('1::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) # error: dependent vars repeated
1786    #StoreEquivalence('0:1:eA',('0:0:eA',)) # error: equiv & fixed
1787    #StoreEquivalence('0:99:Scale',('0:12:Scale',)) # error: equiv & constrained
1788    #StoreEquivalence('0:12:Scale',('0:99:Scale',)) # error: equiv & constrained
1789    varylist = ['2::atomx:3',
1790                '2::C(10,6,1)', '1::C(10,6,1)',
1791                '2::atomy:3', '2::atomz:3',
1792                '0:12:Scale', '0:11:Scale', '0:14:Scale', '0:13:Scale', '0:0:Scale']
1793#    e,w = CheckConstraints([,
1794#                     [{'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}],
1795#                     ['1.0'])
1796#    if e: print 'error=',e
1797#    if w: print 'error=',w
1798#    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']
1799#    constrDict = [
1800#        {'0::Afrac:4': 24.0, '0::Afrac:1': 16.0, '0::Afrac:3': 24.0, '0::Afrac:2': 16.0},
1801#        {'0::Afrac:1': 1.0, '0::Afrac:2': 1.0},
1802#        {'0::Afrac:4': 1.0, '0::Afrac:3': 1.0}]
1803#    fixedList = ['40.0', '1.0', '1.0']
1804
1805    errmsg, warnmsg = CheckConstraints(varylist,constrDict,fixedList)
1806    if errmsg:
1807        print ("*** Error ********************")
1808        print (errmsg)
1809    if warnmsg:
1810        print ("*** Warning ********************")
1811        print (warnmsg)
1812    if errmsg or warnmsg:
1813        sys.exit()
1814    groups,parmlist = GroupConstraints(constrDict)
1815    GenerateConstraints(groups,parmlist,varylist,constrDict,fixedList)
1816    print (VarRemapShow(varylist))
1817    parmdict.update( {
1818        '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,
1819        '0:0:eA': 0.0,
1820        '2::C(10,6,1)': 0.2, '1::C(10,6,1)': 0.3,
1821        '1::C(10,0,1)': 0.2, '2::C(10,0,1)': 0.3,
1822        '1::AUiso:0': 0.02, '0::AUiso:0': 0.03,
1823        '0::A0': 0.0,
1824        '2::atomx:3':0.23,'2::atomy:3':-.23, '2::atomz:3':-0.11,
1825        })
1826    print ('parmdict start',parmdict)
1827    print ('varylist start',varylist)
1828    before = parmdict.copy()
1829    Map2Dict(parmdict,varylist)
1830    print ('parmdict before and after Map2Dict')
1831    print ('  key / before / after')
1832    for key in sorted(list(parmdict.keys())):
1833        print ('  '+key,'\t',before.get(key),'\t',parmdict[key])
1834    print ('varylist after',varylist)
1835    before = parmdict.copy()
1836    Dict2Map(parmdict,varylist)
1837    print ('after Dict2Map')
1838    print ('  key / before / after')
1839    for key in sorted(list(parmdict.keys())):
1840        print ('  '+key,'\t',before.get(key),'\t',parmdict[key])
1841#    dMdv = len(varylist)*[0]
1842#    deriv = {}
1843#    for i,v in enumerate(parmdict.keys()): deriv[v]=i
1844#    Dict2Deriv(varylist,deriv,dMdv)
Note: See TracBrowser for help on using the repository browser.