source: trunk/GSASIImapvars.py

Last change on this file was 5038, checked in by toby, 27 hours ago

Major revision to constraints including an update to the Sequential Refinement tutorial

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 111.2 KB
Line 
1# TODO: revisit SeqRefine and :meth:`GSASIIdataGUI.GSASII.OnSeqRefine` and :func:`GSASIIseqGUI.UpdateSeqResults`
2
3# -*- coding: utf-8 -*-
4########### SVN repository information ###################
5# $Date: 2021-09-27 23:02:41 +0000 (Mon, 27 Sep 2021) $
6# $Author: toby $
7# $Revision: 5038 $
8# $URL: trunk/GSASIImapvars.py $
9# $Id: GSASIImapvars.py 5038 2021-09-27 23:02:41Z toby $
10########### SVN repository information ###################
11"""
12*GSASIImapvars: Parameter constraints*
13======================================
14
15Module to implements algebraic contraints, parameter redefinition
16and parameter simplification contraints.
17
18
19*Externally-Accessible Routines*
20---------------------------------
21
22To define a set of constrained and unconstrained relations, one
23defines a list of dictionary defining constraint parameters and their
24values, a list of fixed values for each constraint and a list of
25parameters to be varied. In addition, one uses
26:func:`StoreEquivalence` to define parameters that are equivalent.
27See the :ref:`Constraints Processing section<Constraints_processing>` for details on how
28processing of constraints is done.
29
30=============================  ===================================================================
31  routine                      explanation
32=============================  ===================================================================
33:func:`InitVars`               This is used to clear out all defined previously
34                               defined constraint information
35 
36:func:`StoreEquivalence`       Implements parameter redefinition.
37                               This should be called for every set of equivalence relationships.
38                               Use :func:`StoreEquivalence` before calling
39                               :func:`GenerateConstraints`
40
41:func:`ProcessConstraints`     Initially constraints of all types are maintained in lists of
42                               dict entries that are stored in the data tree,
43                               with parameters are stored as
44                               :class:`~GSASIIobj.G2VarObj` objects so that they can
45                               be resolved if the phase/histogram order changes.
46                               :func:`ProcessConstraints` processes this list of dict entries,
47                               separating the "Equivalence", "Hold", “Const” and “New Var”
48                               entries for subsequent use.
49                               See the :ref:`Constraint Reorganization <ProcessConstraints>`
50                               section for more details.
51
52:func:`EvaluateMultipliers`    Convert any string-specified (formula-based) multipliers to
53                               numbers. Call this before using :func:`GenerateConstraints`.
54                               At present only values in dict for phase (atom/cell) parameters
55                               are used to evaluate multipliers containint formulae,
56                               but this could be changed if needed.
57
58:func:`GenerateConstraints`    Generates the internally-used tables from constraints and
59                               equivalences. Checks for internal consistency and repairs
60                               problems where possible. See the
61                               :ref:`Constraint Checking and Grouping <GenerateConstraints>`
62                               and :ref:`Equivalence Checking<CheckEquivalences>`
63                               sections for more details.
64
65:func:`Map2Dict`               To determine values for any parameters created in this module,
66                               call Map2Dict. This will not apply contraints.
67
68:func:`Dict2Map`               To apply the constraints and equivalences, call this.
69                               It takes values from the new independent parameters and
70                               constraints, and applies them to the parameter dict.
71
72:func:`Dict2Deriv`             This determines derivatives on independent parameters
73                               from those on dependent ones.
74
75:func:`ComputeDepESD`          Use ComputeDepESD to compute uncertainties on dependent variables.
76
77:func:`VarRemapShow`           Use this to show a summary of the parameter remapping.
78                               Call after :func:`GenerateConstraints`.
79=============================  ===================================================================
80
81Types of constraints
82--------------------
83
84There are four ways to specify constraints, as listed below.
85Note that constraints are initially stored in the
86main section of the GSAS-II data tree under heading ``Constraints``.
87This dict has four keys, 'Hist', 'HAP', 'Global', and 'Phase',
88each containing a list of constraints. An additional set of constraints
89are generated for each phase based on symmetry considerations by calling
90:func:`GSASIIstrIO.GetPhaseData`.
91
92Note that in the constraints, as stored in the GSAS-II data tree, parameters
93are stored as :class:`GSASIIobj.G2VarObj` objects, as these objects allow for
94changes in numbering of phases, histograms and atoms since :class:`GSASIIobj.G2VarObj` objects
95use random Id's for references.
96When constraints are interpreted (in :func:`ProcessConstraints`),
97these references are resolved to the numbered objects by looking up random Id's
98so that the parameter object is converted to a string of form ``ph:hst:VARNAM:at``.
99
100Constraints are initially stored as described in the
101:ref:`constraint definitions <Constraint_definitions_table>`, where the last value in the
102list determines which type of constraint is defined.
103
104Alternate parameters (New Var)
105^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
106
107Parameter redefinition ("New Var" constraints)
108is done by creating an expression that relates several
109parameters: ::
110
111   Mx1 * Px + My1 * Py +... = ::newvar1
112   Mx2 * Px + Mz2 * Pz + ... = ::newvar2
113
114where Pj is a GSAS-II parameter name and Mjk is a constant (float) multiplier.
115Alternately, multipliers Mjk can contain a formula (str) that will be evaluated prior
116to the start of the refinement. In a formula, GSAS-II parameters will be replaced by the
117value of the parameter before the formula is evaluated, so ``'np.cos(0::Ax:2)'`` is a valid
118multiplier. At present, only phase (atom/cell) parameters are available for use in
119a formula, from the GUI but this can be expanded if needed.
120
121This type of constraint describes an alternate
122degree of freedom where parameter Px and Py, etc. are varied to keep
123their ratio
124fixed according the expression. A new variable parameter is assigned to each degree of
125freedom when refined. An example where this can be valuable is when
126two 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.
127In the later stages of refinement, a second
128variable, Pd = P1 - P2, can be defined and it can be seen if refining Pd is
129supported by the data. Another use will be to define parameters that
130express "irrep modes" in terms of the fundamental structural parameters.
131
132The "New Var" constraints are stored as a type "f"
133:ref:`constraint (see definitions)<Constraint_definitions_table>`.
134
135Constrained parameters (Const)
136^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
137
138A constraint on a set of variables can be supplied in the form of a
139linear algebraic equation: ::
140
141  Nx1 * Px + Ny1 * Py +... = C1
142
143where Cn is a constant (float), where Pj is a GSAS-II parameter name,
144and where Njk is a constant multiplier (float) or a formula (str) that will be evaluated prior
145to the start of the refinement. In a formula, GSAS-II parameters will be replaced by the
146value of the parameter before the formula is evaluated, so ``'np.cos(0::Ax:2)'`` is a valid
147multiplier. At present, only phase (atom/cell) parameters are available for use in
148a formula, but this can be expanded if needed.
149
150These equations set an interdependence between parameters.
151Common uses of parameter constraints are to set rules that decrease the number of parameters,
152such as restricting the sum of fractional occupancies for atoms that share
153a site to sum to unity, thus reducing the effective number of variables by one.
154Likewise, the Uiso value for a H atom "riding" on a C, N or O atom
155can be related by a fixed offset (the so called B+1 "rule").
156
157The "Const" constraints are stored as a type "c"
158:ref:`constraint (see definitions)<Constraint_definitions_table>`.
159
160Equivalenced parameters (Equiv)
161^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
162
163A simplifed way to set up a constraint equation is to define an equivalence,
164which can be of form: ::
165
166  C1 * P1 = C2 * Py
167
168or::
169
170  C1 * P1 = C2 * P2 = C3 * P3 = ...
171
172where Cn is a constant (float), where Pj is a GSAS-II parameter name. This
173means that parameters Py (or P2 and P3) are determined from (or "slaved" to)
174parameter P1. Alternately, equivalences can be created with :func:`StoreEquivalence`
175where the multipliers can be a formula (str) that will be evaluated prior to the start of
176the refinement. In a formula, GSAS-II parameters will be replaced by the value of the
177parameter before the formula is evaluated, so a multiplier can be specifed as ``'2*np.cos(0::Ax:2)'``.
178At present, only phase (atom/cell) parameters are available for use in
179such a formula, but this can be expanded if needed.
180
181The first parameter (P1 above)
182is considered the independent variable
183and the remaining parameters are dependent variables. The dependent variables
184are then set from the independent variable.
185
186Note that a constraint expression is conceptually identical to
187defining constraint equations.
188The previous set of equalities could also be written as a set of constraint
189equations in this way: ::
190
191  C1 * P1 - C2 * P2 = 0
192  C1 * P1 - C3 * P3 = 0
193  ...
194
195In practice, however,
196equivalenced parameters are processed in a
197different and more direct manner than constraint equations.
198
199A parameter can be used in multiple
200equivalences where it is an independent variable,
201but if a parameter were used both as a dependent and independent variable then the order that
202shifts are applied becomes potentially significant. As an example, in this case these two
203equivalences are "chained"::
204
205  C1 * P1 = C2 * P2
206  C2 * P2 = C3 * P3
207
208where P2 is both a dependent and independent variable. Likewise, if a parameter is used both in equivalences
209and in "New Var" or "Const" constraints, it also becomes unclear how this should be processed. It is
210possible to specify equivalences that conflict with constraints.
211Should parameter be used as both a dependent and an independent variable or if a parameter is used both in
212an the equivalence and in a "New Var" or "Const" constraints, the equivalence
213is converted to a constraint (Const) to
214avoid conflicts. The equivalences that require this are addressed in ::func:`GenerateConstraints` where
215:func:`CheckEquivalences` is used to locate problematic variables in equivalences
216and then change these equivalences to "Const" equations. Also, unneeded equivalences are removed.
217
218For an example of how equivalences may be used, consider
219a material that has **N** O atoms in the asymmetric unit, all in fairly similar bonding environments
220and where the diffraction data are sparse. One may wish to reduce the complexity of the model fit to
221these data by defining Uiso for all O atoms to be the same. This is done by selecting Uiso for any one O atom
222as an independent variable in a equivalence and setting the remaining **N-1** other O atom Uiso
223variables as dependent parameters with multipliers of 1. This will require that all O atom Uiso values
224be identical.
225The results of this refinement will be simpler to understand than if a set of
226constraint equations is used, because the refined parameter (named as ``ph::Uiso:n``) will be the
227independent variable, corresponding to the first O atom and all other variables would be
228expressed in terms of that variable with a single Equivalence expression.
229The alternate would require **N-1** constraint equations, leaving one degree of freedom with a
230variable would that is likely only indirectly related to the Uiso values.
231
232Equivalenced parameters ("EQUIV" constraints), when defined by users,
233or when created to relate phases, are stored as a type "e"
234:ref:`constraint (see definitions)<Constraint_definitions_table>`.
235Symmetry-generated equivalences are generated prior
236to display or refinement in :func:`GSASIIstrIO.GetPhaseData`.
237These are not stored in the data tree.
238
239Hold parameters (Fixed)
240^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
241
242When parameters are refined where a single refinement flag determines that several variables
243are refined at the same time (examples are: cell parameters, atom positions, anisotropic
244displacement parameters, magnetic moments,...) it can be useful to specify that a
245specific parameter should not be varied. These will most commonly be generated due to symmetry,
246but under specific conditions, there may be other good reasons to constrain a parameter.
247
248The "Hold" constraints are stored as a type "h"
249:ref:`constraint (see definitions)<Constraint_definitions_table>`.
250
251.. _Constraints_processing:
252
253Constraint Processing
254---------------------
255
256When constraints will be used or edited, they are processed using a series of
257calls. This is done in GSAS-II from several locations:
258
259* For error checking from the tree in :mod:`GSASIIconstrGUI`,
260  :func:`GSASIIconstrGUI.CheckConstraints` loads constraints from
261  the data tree.
262
263* When the table of refined parameters is shown, constraints are also
264  processed in function :func:`GSASIIdataGUI.GSASII.OnShowLSParms` using
265  :func:`GSASIIconstrGUI.CheckConstraints`
266
267* To write parameters in the Export sections of the program,
268  :func:`GSASIIIO.loadParmDict` loads results as well as constraints
269  from the tree. This works a bit differently from the above, so it
270  makes direct calls to the constraints routines.
271
272* For error checking from a GPX file
273  :func:`GSASIIstrIO.ReadCheckConstraints` loads constraints
274  (called in :mod:`GSASIIdataGUI` and :mod:`GSASIIscriptable`),
275  which is similar to :func:`GSASIIconstrGUI.CheckConstraints`.
276  :func:`~GSASIIstrIO.ReadCheckConstraints` is called by
277  :meth:`GSASIIdataGUI.GSASII.OnRefine` and
278  :meth:`GSASIIdataGUI.GSASII.OnSeqRefine`
279  before constraints are generated for use in refinements so they can
280  be shown in the GUI. This is also called to check for errors in
281  :class:`GSASIIscriptable.G2Project`.
282
283* To create the constraints for use in a refinement, in
284  :mod:`GSASIIstrMain`, functions :func:`GSASIIstrMain.Refine` and
285  :func:`GSASIIstrMain.SeqRefine` load and process the constraints again.
286  This is repeated here because :func:`~GSASIIstrMain.Refine` and
287  :func:`~GSASIIstrMain.SeqRefine` are intended to operate as stand-alone
288  routines that may be called directly.
289
290* After sequential fits have been completed, the previously processed
291  constraint info is read from the sequential results section of the
292  data tree. Function
293  :func:`GSASIIseqGUI.UpdateSeqResults` displays the sequential results
294  table
295also processes constraints.
296
297TODO: Note that G2stIO.makeTwinFrConstr is called only in one place. It probably needs to be included in all of the above.
298
299When constraints are processed, the following steps are used: 
300
301#. Constraints are stored in separate lists in the data tree to
302   simplify their creation and their GUI display.
303   In the initial processing, all of the stored constraints are appended
304   into a single list.
305
306#. Then :func:`InitVars` is used to initialize the global variables in
307   this module (:mod:`GSASIImapvars`). This may be done before the previous
308   step.
309
310#. Then :func:`ProcessConstraints` is used to initially process the
311   constraints user-supplied constraints (from the data tree),
312   as described in :ref:`Constraint Reorganization <ProcessConstraints>`.
313   When constraints are read from a GPX file, rather than the data tree, use
314   :func:`GSASIIstrIO.ReadConstraints` (which calls :func:`ProcessConstraints`).
315
316#. Symmetry-generated equivalences are then created in
317   :func:`GSASIIstrIO.GetPhaseData`, which also calls
318   :func:`GSASIIstrIO.cellVary` and for Pawley refinements
319   :func:`GSASIIstrIO.GetPawleyConstr`. These are entered directly into this
320   module's globals using :func:`StoreEquivalence`.
321
322#. Constraints/equivalences are then checked for possible conflicts with
323   :func:`GenerateConstraints`, this requires grouping the constraints,
324   as described below.
325
326#. :func:`GenerateConstraints` is then called to
327   create the constraints that will be used,
328   :ref:`see below <GenerateConstraints>` for more details.
329
330#. For debugging constraints, :func:`VarRemapShow` can be called after
331   :func:`GenerateConstraints` to display the generated constraints.
332
333.. _ProcessConstraints:
334
335Constraint Reorganization (:func:`ProcessConstraints`)
336^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
337
338:func:`ProcessConstraints` is used to initially process the
339constraints from the list of dict entries. The "Const" and "New Var" are placed into two
340lists (:data:`constrDict` and :data:`fixedList`) that are later used for parameter
341grouping (in :func:`GenerateConstraints`). "Hold" and "Equivalence" constraints are
342separated into separate storage.
343 
344   For "**Const**" entries,
345     a dict with multiple entries is placed in :data:`constrDict` where
346     each dict key is the parameter name and the value is the multiplier for the parameter,
347     while :data:`fixedList` gets a string value corresponding to the constant value for
348     the expression.
349
350   For "**New Var**" entries,
351     a dict with multiple entries defined identically to
352     to that used in "Const" entries. The differences between "New Var" and "Const" entries is
353     that for "Const" entries, a constant value (as a string) is placed in :data:`fixedList` while
354     for "New Var" entries corresponding entry in :data:`fixedList` is None.
355     Also, one or two additional entries are created in the dict for "New Var" constraints:
356     an entry with key "_vary" is given the value of True or False
357     depending on the refinement flag setting;
358     an entry with key "_name" will be created if the "New Var" parameter has a supplied name.
359
360   For "**Hold**" entries,
361     the “Hold” constraints are stored in global variables :data:`holdParmList` by calling
362     :func:`StoreHold`; holds that are generated by this code are stored in :data:`newHolds`.
363
364   **Equivalences** are stored using :func:`StoreEquivalence` into this module's globals
365     (:data:`dependentParmList`, :data:`arrayList`, :data:`invarrayList`, :data:`indParmList`,
366     and :data:`symGenList`).
367     For each equivalence:
368
369     * a list with one entry, the name of the independent parameter is placed in :data:`~GSASIImapvars.indParmList`;
370     * a list with one or more parameter name is placed in :data:`~GSASIImapvars.dependentParmList`;
371     * the value None is added to :data:`~GSASIImapvars.arrayList`;
372     * a list of multipliers for each dependent variable is placed in :data:`~GSASIImapvars.invarrayList`
373     * an entry of either True or False is placed in :data:`~GSASIImapvars.symGenList`, where True indicates that the entry has been generated from symmetry.
374
375The output from :func:`ProcessConstraints` will have the form as below,
376where the first entry is a "Const", the second is a "New Var" and the final is a "Hold".
377
378  .. code-block:: python
379
380    constrDict = [
381         {'0:12:Scale': 2.0, '0:14:Scale': 4.0, '0:13:Scale': 3.0, '0:0:Scale': 0.5},
382         {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0, '_vary':True},
383         {'0::A0': 0.0}]
384    fixedList = ['5.0', None, '0']
385
386.. _GenerateConstraints:
387
388Constraint Checking and Grouping (:func:`GenerateConstraints`)
389^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
390
391Function :func:`GenerateConstraints` is used to
392process the parameter equivalences and constraint lists created in
393:func:`ProcessConstraints` (``constrDict`` and ``fixedList``). :func:`GenerateConstraints`
394is used to generate error/warning messages, to set up lists that are used to show this
395information for the GUI (using :func:`getConstrError`) and to
396generate the information stored in :ref:`global arrays <GlobalVariables>` that are used later to
397apply the constraints.
398
399When a sequential refinement is in progress, the constraints are scanned for parameters
400that have a wildcard (*) for the histogram number, such as 1:*:Scale which would refer
401to the phase fraction for Phase ` in every histogram. The "*" will be replaced with
402the number of the current histogram.
403
404Equivalences are checked with :func:`CheckEquivalences` (described in detail
405:ref:`below <CheckEquivalences>`). This may result in the creation of additional "Hold"
406and "Constr" constraints being added to the ``constrDict`` and ``fixedList`` lists.
407
408The "Const" and "New Var" constraint expressions are then scanned for problems:
409
410Constraints cannot be processed without changes if any of the terms within have the following:
411
412* **Hold (Fixed) parameters**; **Unvaried parameters**; **Multiplier of zero**
413
414  If all parameters in a constraint are either not refined, or are marked as "Hold", or have a multiplier
415  that evaluates as 0.0, the constraint can be ignored.
416
417  If all but one parameter in a constraint are either not refined, or are marked as "Hold", or have a multiplier
418  that evaluates as 0.0, the constraint determines a static value for the parameter. The appropriate
419  constant values are substituted into the constraint equation and the one remaining parameter
420  is set as "Hold".
421
422  If two or more parameter remain in a constraint after removing parameters that are not refined, or
423  are marked as "Hold", or have a multiplier that evaluates as 0.0, the constraint can still be used.
424  The appropriate constant values are substituted into the constraint equation, any zero multiplier
425  terms set as "Hold" and the remaining terms in the constraint are used.
426
427* **Undefined parameters**
428
429  If all parameters in a constraint are undefined, the constraint can be ignored.
430
431  If some, but not all, parameters in a constraint are undefined, the constraint cannot be processed.
432  The remaining parameters will be set as "Hold" and the constraint will be ignored. One exception:
433  atom position constraints (p::dA[xyz]:#) will be assumed as zero.
434
435Constraint expressions ("Const" and "New Var") are sorted by routine :func:`GroupConstraints` into
436groups so that each group contains the minimum number of entries that
437ensures each parameter is referenced in only one group.
438This is done by scanning the
439list of dicts in :data:`constrDict` one by one and making a list
440of parameters used in that constraint expression. Any expression that contains
441a parameter in that list is added to the current group and those
442parameters are added to this list of parameters. The list of ungrouped
443expressions is then scanned again until no more expressions are added to the
444current group. This process is repeated until every expression has been
445placed in a group. Function :func:`GroupConstraints` returns two lists of lists.
446The first has, for each group, a list of the indices in :data:`constrDict`
447that comprise the group (there can be only one). The second list contains,
448for each group, the unique parameter names in that group.
449
450Each constraint group is then processed. First, wildcard parameters are
451renamed (in a sequential refinement). Any held parameters that are used
452in constraints are noted as errors. The number of refined parameters and
453the number of parameters that are not defined in the current refinement are
454also noted. It is fine if all parameters in a group are not defined or all are
455not varied, but if some are defined and others not or some are varied and
456others not, this constitutes an error.
457
458The contents of each group is then examined. Groups with a single
459parameter (holds) are ignored. Then for each group, the number
460of parameters in the group (Np) and the number of expressions in the
461group (Nc) are counted and for each expression. If Nc > Np, then the constraint
462is overdetermined, which also constitutes an error.
463
464The parameter multipliers for each expression are then assembled:
465
466::
467
468   M1a * P1 + M2a * P2 +... Mka * Pk
469   M1b * P1 + M2b * P2 +... Mkb * Pk
470   ...
471   M1j * P1 + M2j * P2 +... Mkj * Pk
472
473
474From this it becomes possible to create a Nc x Np matrix, which is
475called the constraint matrix:
476
477 .. math::
478
479   \\left( \\begin{matrix}
480   M_{1a}  & M_{2a} &... & M_{ka} \\\\
481   M_{1b}  & M_{2b} &... & M_{kb} \\\\
482   ...  \\\\
483   M_{1j}  & M_{2j}  &... & M_{kj}
484   \\end{matrix}\\right)
485
486When Nc<Np, then additional rows need to be added to the matrix and to
487the vector that contains the value for each row (:data:`fixedList`) where
488values are ``None`` for New Vars and a constant for fixed values.
489This then can describe a system of Np simultaneous equations:
490
491 .. math::
492
493   \\left( \\begin{matrix}
494   M_{1a}  & M_{2a} &... & M_{ka} \\\\
495   M_{1b}  & M_{2b} &... & M_{kb} \\\\
496   ...  \\\\
497   M_{1j}  & M_{2j}  &... & M_{kj}
498   \\end{matrix}\\right)
499   \\left( \\begin{matrix}
500   P_{1} \\\\
501   P_{2} \\\\
502   ...  \\\\
503   P_{k}
504   \\end{matrix}\\right)
505   =
506   \\left( \\begin{matrix}
507   C_{1} & C_{2} &  ... & C_{k}
508   \\end{matrix}\\right)
509
510This is done by creating a square matrix from the group using
511:func:`_FillArray` with parameter ``FillDiagonals=False`` (the default). Any
512unspecified rows are left as all zero. The first Nc rows in the
513array are then coverted to row-echelon form using :func:`_RowEchelon`. This
514will create an Exception if any two rows are linearly dependent (which means
515that no matter what values are used for the remaining rows, that the matrix
516will be singular). :func:`_FillArray` is then called with parameter
517``FillDiagonals=True``, which again creates a square matrix but where
518unspecified rows are zero except for the diagonal elements. The 
519`Gram-Schmidt process <http://en.wikipedia.org/wiki/Gram-Schmidt>`_,
520implemented  in :func:`GramSchmidtOrtho`, is used to find orthonormal unit
521vectors for the remaining Np-Nc rows of the matrix. This will fail with
522a ``ConstraintException`` if this is not possible (singular matrix) or
523the result is placed in :data:`constrArr` as a numpy array.
524
525Rows in the matrix corresponding to "New Var" constraints and those that
526were generated by the Gram-Schmidt process are provided with parameter names
527(this can be specified if a "New Var" entry by using a ``"_name"`` element
528in the constraint dict, but at present this is not implemented.) Names are
529generated using :data:`paramPrefix` which is set to ``"::constr"``, plus a
530number to make the new parameter name unique.
531
532Finally the parameters used as input to the constraint are placed in
533this module's globals
534:data:`~GSASIImapvars.dependentParmList` and the constraint matrix is
535placed in into  :data:`~GSASIImapvars.arrayList`. This can be used to compute
536the initial values for "New Var" parameters. The inverse of the
537constraint matrix is placed in :data:`~GSASIImapvars.invarrayList` and a list
538of the "New Var" parameters and a list of the fixed values (as str's)
539is placed in :data:`~GSASIImapvars.indParmList`.
540Finally the appropriate entry in :data:`~GSASIImapvars.symGenList` is set to
541False to indicate that this is not a symmetry generated constraint.
542
543.. _CheckEquivalences:
544
545Equivalence Checking and Reorganization (:func:`CheckEquivalences`)
546^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
547
548Equivalences need to be checked for usages that could be internally conflicted
549or have possible conflicts with other constraints.
550
551**Mixed parameter use:**
552
553 **Note** that cycling through the equivalences may be needed to find all mixed-use
554 parameters, see below.
555
556* A parameter should not show up as a dependent variable in two equivalence expressions,
557  such as::
558
559      ::x1 -> ::x3
560      ::x2 -> ::x3
561
562  This will be processed by turning the equivalences into two constraint equations::
563
564      ::x1 - ::x3 = 0
565      ::x2 - ::x3 = 0
566
567  which can be satisfied when ``::x1 = ::x2 = ::x3``. If  ``::x1`` and ``::x2`` had been
568  intended to be independent parameters, then the above equivalences would be conflict and
569  cannot be statisfied.
570
571* If a parameter is used both as an independent and as a dependent variable (*chaining*),
572  as is in these two equivalence expressions::
573
574      ::x1 -> ::x2 & ::x4
575      ::x2 -> ::x3
576
577  This can also be addressed by turning these equivalences into three constraint equations::
578
579   ::x1 - ::x2 = 0
580   ::x1 - ::x4 = 0
581   ::x2 - ::x3 = 0
582
583  which can be satisfied when ``::x1 = ::x2 = ::x3 = ::x4``
584
585*  Use of parameters in both equivalences and "Const" or "New Var" constraint expressions makes
586   logical sense::
587
588   ::x1 -> ::x2 & ::x4
589   ::x2 + ::x3 = 0
590
591   This can also be addressed by turning the equivalence into two constraint equations::
592
593   ::x1 - ::x2 = 0
594   ::x1 - ::x4 = 0
595
596   With the addition of the "Const" equation (``::x2 + ::x3 = 0``), the solution will require
597   ``::x1 = ::x2 = -1.0*::x3 = ::x4``
598
599* Cycling is needed to find all equivalences that must be converted.
600  Consider this set of constraints::
601
602   ::x2 + ::x3 = 0
603   ::x1 -> ::x2
604   ::x1 -> ::x4
605
606 In the first pass the equivalence with ``::x2`` would be converted to a "Const" constraint
607 and in the second pass
608 the other equivalence with ``::x1`` would be converted.
609
610
611**Mixing Hold (Fixed) parameters in equivalences**
612
613* If one parameter is designated as a "Hold" in an equivalence, then all parameters in that
614  equivalence cannot be varied. Considering this equivalence::
615
616      ::x1 -> ::x2 & ::x4
617
618  If any of the three parameters (``::x1``, ``::x2``, or `::x4`) are marked as Hold, then
619  the other two parameters may not be varied and will also be set with a "Hold".
620
621
622**Unvaried parameters in equivalences**
623
624* If no parameters in an equivalence are varied, then the equivalence is ignored.
625
626* If only some parameters are marked as varied then
627  *none of the parameters can be varied*; any varied parameters will be set with a "Hold".
628
629
630**Undefined parameters in equivalences**
631
632  Parameters may be placed in equivalences that are not actually defined in a project.
633  This can occur in two ways. If an equivalence is created in the GUI for a parameter that
634  is later supplanted with a different model (for example, changing from isotropic size
635  broadening to uniaxial broadening replaces the isotropic broadening term with two different
636  uniaxial terms) or symmetry may require restrictions on anisotropic ADPs that are not
637  in use).
638
639* If the independent parameter is undefined, then any dependent parameters that are defined
640  are set as "Hold" and the equivalence is ignored.
641
642* If all dependent parameters are undefined, then the equivalence is ignored.
643
644* If a dependent parameter is undefined, then that parameter is dropped from the equivalence.
645
646
647**Multiplier of zero in equivalences**
648
649  Any dependent parameter that has a multiplier of zero will be dropped from the equivalence.
650  If no terms remain, then the equivalence is ignored. (Independent parameters do not
651  have a multiplier).
652
653
654.. _GlobalVariables:
655
656*Global Variables*
657------------------
658
659This module uses a number of global variables. One set is used to store the
660constraints and equivalences after processing by :func:`StoreEquivalence` and
661:func:`GenerateConstraints`. 
662These globals are expected to be used only by this module's (:mod:`GSASIImapvars`) internal routines.
663
664
665.. tabularcolumns:: |l|p{4.5in}|
666
667=============================  ===================================================================
668  variable                      explanation
669=============================  ===================================================================
670:data:`dependentParmList`        a list containing group of lists of
671                                 parameters used in the group. Note that parameters listed in
672                                 dependentParmList will not be included in the Hessian as their
673                                 derivatives will not affect the model
674
675:data:`indParmList`              a list containing groups of Independent parameters defined in
676                                 each group. This contains both parameters used in parameter
677                                 redefinitions as well as names of generated new parameters.
678
679:data:`arrayList`                a list containing group of relationship matrices to relate
680                                 parameters in dependentParmList to those in indParmList.
681
682:data:`invarrayList`             a list containing group of relationship matrices to relate
683                                 parameters in indParmList to those in dependentParmList.
684                                 Unlikely to be used externally.
685
686:data:`holdParmList`             a list of parameters that have been marked as "Hold".
687                                 Unlikely to be used externally. Holds have a single
688                                 entry in the indParmList parameter list.
689                                 Set in :func:`StoreHold`. (Also see :data:`newHolds`)
690
691:data:`symGenList`               a list of boolean values that will be True to indicate
692                                 that an equivalence was generated internally GSAS-II
693                                 meaning it is generated based on symmetry, twining
694                                 or Pawley overlap.
695
696:data:`dependentVars`            a list of dependent variables in equivalences, compiled
697                                 from (:data:`dependentParmList`).
698                                 Used within :func:`GetDependentVars`.
699
700:data:`independentVars`          a list of dependent variables in equivalences, compiled
701                                 from (:data:`indParmList`).
702                                 Used within :func:`GetIndependentVars`.
703
704:data:`saveVaryList`             a list of the varied parameters used when constraints
705                                 were last processed.
706
707=============================  ===================================================================
708
709
710A second set of global variables are set in :func:`GenerateConstraints` with lists of parameter
711names from equivalences and constraints. Used in :func:`CheckEquivalences` and
712:func:`getConstrError`.
713
714.. tabularcolumns:: |l|p{4.5in}|
715
716=============================  ===================================================================
717  variable                      explanation
718=============================  ===================================================================
719:data:`depVarList`             a list of the parameters used in equivalences as dependent
720                               parameters for all equivalences initially specified (including
721                               those to be reclassified as "Constr" constraints.)
722:data:`indepVarList`           a list of the parameters used in equivalences as independent
723                               parameters for all equivalences initially specified (including
724                               those to be reclassified as "Constr" constraints.)
725:data:`constrVarList`          a list of the parameters that are used in "Constr" or
726                               "New Var" constraints. Does not include those in equivalences
727                               to be reclassified as "Constr" constraints.)
728=============================  ===================================================================
729
730
731.. tabularcolumns:: |l|p{4.5in}|
732
733A third set of global variables to store equivalence warning information.
734Set in :func:`CheckEquivalences` and :func:`GenerateConstraints`.
735Used in :func:`getConstrError` to display warning messages.
736
737=============================  ===================================================================
738  variable                      explanation
739=============================  ===================================================================
740:data:`convVarList`            parameters in equivalences that were converted to "Const"
741                               constraints
742:data:`multdepVarList`         parameters used as dependent parameters in equivalences
743                               multiple times
744:data:`newHolds`               parameters to be added as "Hold"s
745:data:`unvariedParmsList`      parameters used in equivalences and constraints
746                               that are not varied
747:data:`undefinedVars`          parameters used in equivalences
748                               that are not defined in the parameter dict (parmDict)
749:data:`groupErrors`            parameters in constraints that cause grouping errors
750=============================  ===================================================================
751
752
753
754*GSASIImapvars Routines/variables*
755---------------------------------------
756
757Note that parameter names in GSAS-II are strings of form ``<ph#>:<hst#>:<nam>`` or ``<ph#>::<nam>:<at#>``
758where ``<ph#>`` is a phase number, ``<hst#>`` is a histogram number and ``<at#>`` is an atom number.
759``<nam>`` is a name that determines the parameter type (see :func:`GSASIIobj.CompileVarDesc`). When
760stored in the data tree, parameters are saved as :class:`GSASIIobj.G2VarObj` objects
761so that they can be resolved if the phase/histogram order changes.
762
763"""
764
765from __future__ import division, print_function
766import copy
767import numpy as np
768import GSASIIpath
769GSASIIpath.SetVersionNumber("$Revision: 5038 $")
770import GSASIIobj as G2obj 
771# data used for constraints;
772debug = False # turns on printing as constraint input is processed
773
774#------------------------------------------------------------------------------------
775# Global vars used for storing constraints and equivalences after processing
776#   note that constraints are stored listed by contraint groups,
777#   where each constraint
778#   group contains those parameters that must be handled together
779
780dependentParmList = []
781'''a list of lists where each item contains a list of parameters in each constraint group.
782note that parameters listed in dependentParmList should not be refined directly.'''
783indParmList = [] # a list of names for the new parameters
784'''a list of lists where each item contains a list for each constraint group with
785fixed values for constraint equations and names of generated (New Var) parameters.
786'''
787arrayList = []
788'''a list of of relationship matrices that map model parameters in each
789constraint group (in :data:`dependentParmList`) to
790generated (New Var) parameters.
791'''
792invarrayList = []
793'''a list of of inverse-relationship matrices that map constrained values and
794generated (New Var) parameters (in :data:`indParmList`) to model parameters
795(in :data:`dependentParmList`).
796'''
797holdParmList = []
798'''List of parameters that should not be refined ("Hold"s).
799Set in :func:`StoreHold`. Initialized in :func:`InitVars`.
800'''
801symGenList = []
802'''A list of flags that if True indicates a constraint was generated by symmetry
803'''
804dependentVars = []
805'''A list of dependent variables in equivalences, compiled from (:data:`dependentParmList`).
806Used within :func:`GetDependentVars`
807'''
808independentVars = []
809'''A list of dependent variables in equivalences, compiled from (:data:`indParmList`).
810Used within :func:`GetIndependentVars`
811'''
812saveVaryList = []
813'''A list of the varied parameters that was last supplied when constraints were
814processed. This is set in :func:`GenerateConstraints` and updated in
815:func:`Map2Dict`. Used in :func:`VarRemapShow`
816'''
817#------------------------------------------------------------------------------------
818# Global vars set in :func:`GenerateConstraints`. Used for intermediate processing of
819# constraints.
820
821constrVarList = []
822'List of parameters used in "Constr" and "New Var" constraints'
823indepVarList = []
824'A list of all independent parameters in equivalences'
825depVarList = []
826'A list of all dependent parameters in equivalences'
827
828#------------------------------------------------------------------------------------
829# global variables is used in :func:`getConstrError` to store error and warning information.
830# set in CheckEquivalences and in GenerateConstraints
831convVarList = []
832'parameters in equivalences that were converted to "Const" constraints'
833multdepVarList = []
834'parameters used as dependents multiple times in equivalences'
835newHolds = []
836'parameters that will be added as "Hold"s based on use in equivalences and constraints' 
837unvariedParmsList = []
838'parameters used in equivalences that are not varied'
839undefinedVars = []
840'parameters used in equivalences that are not defined in the parameter dict'
841groupErrors = []
842'parameters in constraints where parameter grouping and matrix inversion fails'
843
844#------------------------------------------------------------------------------------
845paramPrefix = "::constr"
846'A prefix for generated parameter names'
847consNum = 0
848'The number to be assigned to the next constraint to be created'
849
850#------------------------------------------------------------------------------------
851class ConstraintException(Exception):
852    '''Defines an Exception that is used when an exception is raised processing constraints.
853    Raised in :func:`GenerateConstraints` during sequential fits. Possible (but highly unlikely)
854    to be raised in :func:`CheckEquivalences` (called by :func:`GenerateConstraints`) if an
855    infinite loop is detected.
856    Also raised in :func:`GramSchmidtOrtho` and :func:`_SwapColumns` but caught
857    within :func:`GenerateConstraints`.
858    '''
859    pass
860
861def InitVars():
862    '''Initializes all constraint information'''
863    global dependentParmList,arrayList,invarrayList,indParmList,consNum,symGenList
864    dependentParmList = [] # contains a list of parameters in each group
865    arrayList = [] # a list of of relationship matrices
866    invarrayList = [] # a list of inverse relationship matrices
867    indParmList = [] # a list of names for the new parameters
868    consNum = 0 # number of the next constraint to be created
869    symGenList = [] # Flag if constraint is generated by symmetry
870    global holdParmList
871    holdParmList = []
872
873def VarKeys(constr):
874    """Finds the keys in a constraint that represent parameters
875    e.g. eliminates any that start with '_'
876
877    :param dict constr: a single constraint entry of form::
878
879        {'var1': mult1, 'var2': mult2,... '_notVar': val,...}
880
881        (see :func:`GroupConstraints`)
882    :returns: a list of keys where any keys beginning with '_' are
883      removed.
884    """
885    return [i for i in constr.keys() if not i.startswith('_')]
886
887
888def GroupConstraints(constrDict):
889    """Divide the constraints into groups that share no parameters.
890
891    :param dict constrDict: a list of dicts defining relationships/constraints
892
893    ::
894   
895       constrDict = [{<constr1>}, {<constr2>}, ...]
896
897    where {<constr1>} is {'var1': mult1, 'var2': mult2,... }
898
899    :returns: two lists of lists:
900   
901      * a list of grouped contraints where each constraint grouped containts a list
902        of indices for constraint constrDict entries
903      * a list containing lists of parameter names contained in each group
904     
905    """
906    assignedlist = [] # relationships that have been used
907    groups = [] # contains a list of grouplists
908    ParmList = []
909    for i,consi in enumerate(constrDict):
910        if i in assignedlist: continue # already in a group, skip
911        # starting a new group
912        grouplist = [i,]
913        assignedlist.append(i)
914        groupset = set(VarKeys(consi))
915        changes = True # always loop at least once
916        while(changes): # loop until we can't find anything to add to the current group
917            changes = False # but don't loop again unless we find something
918            for j,consj in enumerate(constrDict):
919                if j in assignedlist: continue # already in a group, skip
920                if len(set(VarKeys(consj)) & groupset) > 0: # true if this needs to be added
921                    changes = True
922                    grouplist.append(j)
923                    assignedlist.append(j)
924                    groupset = groupset | set(VarKeys(consj))
925        group = sorted(grouplist)
926        varlist = sorted(list(groupset))
927        groups.append(group)
928        ParmList.append(varlist)
929    return groups,ParmList
930
931def GenerateConstraints(varyList,constrDict,fixedList,parmDict=None,SeqHist=None):
932    '''Takes a list of relationship entries that have been stored by
933    :func:`ProcessConstraints` into lists ``constrDict`` and ``fixedList``
934
935    This routine then calls :func:`CheckEquivalences`
936    for internal consistency. This includes converting equivalenced variables into
937    constraints when a variable is used in both.
938
939    Once checked, parameters are grouped so that any parameters that are used in
940    more than one constraint are grouped together. This allows checking for incompatible
941    logic (for example, when four constraints are specified for three variables).
942
943    If parmDict is not None, the parameter groups are checked for constraints where
944    some parameters are varied, but not others. If so, the value for that unvaried
945    parameter is subtracted from the constant in the constraint.
946
947    Once all checks are complete, the constraints are then
948    converted to the form used to apply them, saving them as global variables within
949    this module.
950
951    :param list varyList: a list of parameters names (strings of form
952      ``<ph>:<hst>:<nam>``) that will be varied. Note that this is changed here.
953   
954    :param dict constrDict: a list of dicts defining relationships/constraints
955      (as described in :func:`GroupConstraints`)
956
957    :param list fixedList: a list of values specifying a fixed value for each
958      dict in constrDict. Values are either strings that can be converted to
959      floats, float values or None if the constraint defines a new parameter.
960     
961    :param dict parmDict: a dict containing all parameters defined in current
962      refinement.
963
964    :param int SeqHist: number of current histogram, when used in a sequential
965      refinement. None (default) otherwise. Wildcard parameter names are
966      set to the current histogram, when found if not None.
967
968    :returns: errmsg,warning,groups,parmlist
969
970      **errmsg**
971        Is an error message or empty if no errors were found
972      **warning**
973        Is a warning message about constraints that have been ignored or changed
974      **groups**
975        Lists parameter groups
976      **parmlist**
977        Lists parameters in each parameter groups
978    '''
979    warninfo = {'msg':'', 'shown':{}}
980    def warn(msg,cdict=None,val=None):
981        if cdict is not None and cdict != warninfo['shown']:
982            warninfo['shown'] = cdict
983            if warninfo['msg']: warninfo['msg'] += '\n'
984            warninfo['msg'] += '\nProblem with constraint: ' + _FormatConstraint(cdict,val)
985        if warninfo['msg']: warninfo['msg'] += '\n'
986        warninfo['msg'] += '  ' + msg
987   
988    global dependentParmList,arrayList,invarrayList,indParmList,consNum
989    # lists of parameters used for error reporting
990    global undefinedVars # parameters that are used in equivalences but are not defined
991    undefinedVars = []
992    global newHolds # additional parameters that should be held
993    newHolds = []
994    global groupErrors # parameters in constraints that cause grouping errors
995    groupErrors = []
996    global saveVaryList
997    saveVaryList = copy.copy(varyList)
998
999    errmsg = ''  # save error messages here. If non-blank, constraints cannot be used.
1000    warning = '' # save informational text messages here.
1001
1002    # tabulate a list of all parameters as defined initially setting the following global variables:
1003    global depVarList # parameters used in equivalences as dependent parameters
1004    global indepVarList # parameters used in equivalences as independent parameters
1005    global constrVarList  # list of parameters used in other constraints
1006    constrVarList = []
1007    for cdict in constrDict:
1008        constrVarList += [i for i in cdict if i not in constrVarList and not i.startswith('_')]
1009    # tabulate all parameters in equivalences by type and look for repeated dependent vars
1010    global depVarList,indepVarList
1011    depVarList = []  # list of all dependent parameters in equivalences
1012    indepVarList = []  # list of all independent parameters in equivalences
1013    for cnum,(varlist,mapvars,multarr,invmultarr) in enumerate(zip(
1014            dependentParmList,indParmList,arrayList,invarrayList)):
1015        if multarr is not None: continue # equivalence
1016        indepVarList += [mv for mv in mapvars if mv not in indepVarList]
1017        depVarList += [v for v in varlist if v not in depVarList]
1018   
1019    # Translate & generate lookup table for wildcard parameter names (sequential fits only)
1020    # TODO: where should this be done in the sequence of events? Probably here, so that
1021    # variable names can be looked up in varyList properly. Also, should names be changed in place
1022    # or stored in translateTable and then looked up later?     
1023    #translateTable = {} # lookup table for wildcard-referenced parameters
1024    if SeqHist is not None:
1025        for varlist,mapvars,multarr,invmultarr in zip(
1026            dependentParmList,indParmList,arrayList,invarrayList):
1027            for i,mv in enumerate(mapvars):
1028                if mv.split(':')[1] == '*':
1029                    # convert wildcard var to reference current histogram; save translation in table
1030                    sv = mv.split(':')
1031                    sv[1] = str(SeqHist)
1032                    #mapvars[i][1] = translateTable[mv] = ':'.join(sv)
1033                    mapvars[i][1] = ':'.join(sv)
1034            for i,(v,m) in enumerate(zip(varlist,invmultarr)):
1035                if v.split(':')[1] == '*':
1036                    # convert wildcard var to reference current histogram; save translation in table
1037                    sv = v.split(':')
1038                    sv[1] = str(SeqHist)
1039                    #varlist[i] = translateTable[v] = ':'.join(sv)
1040                    varlist[i] = ':'.join(sv)
1041        for rel in constrDict:
1042            for i,var in enumerate(rel):
1043                if var.startswith('_'): continue
1044                if var.split(':')[1] == '*': # convert wildcard var to reference current histogram
1045                    sv = var.split(':')
1046                    sv[1] = str(SeqHist)
1047                    #rel[i][1] = translateTable[var] = ':'.join(sv)   # add to translation in table for later use
1048                    rel[i][1] = ':'.join(sv)
1049
1050    # Process the equivalences; If there are conflicting parameters, move them into constraints
1051    warning = CheckEquivalences(constrDict,varyList,fixedList,parmDict)
1052
1053    # look through "Constr" and "New Var" constraints looking for zero multipliers and
1054    # Hold, Unvaried & Undefined parameters
1055    skipList = []
1056    for cnum,(cdict,val) in enumerate(zip(constrDict,fixedList)):
1057        valid = 0          # count of good parameter
1058        # error reporting
1059        zeroList = []      # parameters with zero multipliers
1060        holdList = []      # parameters with "Hold"'s
1061        noVaryList = []    # parameters not varied
1062        noWildcardList = [] # wildcard parameters in non-sequential fit
1063        notDefList = []    # parameters not defined
1064        # processing to be done
1065        problem = False    # constraint must be dropped
1066        dropList = []      # parameters to remove
1067        setAsHoldList = [] # parameters that will be held if the constraint is invalid
1068        for i,var in enumerate(cdict):
1069            if var.startswith('_'): continue
1070            if cdict[var] == 0:    # zero multiplier
1071                if var not in zeroList: zeroList.append(var)
1072                setAsHoldList.append(var)
1073                dropList.append(var)
1074            elif var in holdParmList:   # already hold
1075                #if var not in newHolds: newHolds.append(var)
1076                holdList.append(var)
1077                dropList.append(var)
1078            elif ':*:' in var and SeqHist is None:  # unvaried
1079                noWildcardList.append(var)
1080                dropList.append(var)
1081            elif var not in varyList:  # unvaried
1082                if var not in unvariedParmsList: unvariedParmsList.append(var)
1083                noVaryList.append(var)
1084                setAsHoldList.append(var)
1085                dropList.append(var)
1086            elif parmDict is not None and var not in parmDict: # not defined, constraint will not be used
1087                if var not in undefinedVars: undefinedVars.append(var)
1088                notDefList.append(var)
1089                if ':dAx:' in var or ':dAy:' in var or ':dAz:' in var: # undefined atoms
1090                    dropList.append(var) # these can be ignored
1091                else:
1092                    problem = True
1093            else:
1094                valid += 1
1095        for l,m in ((zeroList,"have zero multipliers"), # show warning
1096                      (holdList,'set as "Hold"'),
1097                      (noVaryList,"not varied"),
1098                      (noWildcardList,"wildcard in non-sequential fit"),
1099                      (notDefList,"not defined")):
1100            if l:
1101                msg = "parameter(s) " + m + ': '
1102                for i,v in enumerate(l):
1103                    if i != 0: msg += ', '
1104                    msg += v
1105                warn(msg,cdict,val)
1106        if not valid and (problem or len(dropList) > 0): # no valid entries
1107            warn('Ignoring constraint',cdict,val)
1108            skipList.append(cnum)
1109        elif problem: # mix of valid & refined and undefined items, cannot use this
1110            warn('Holding varied items & Dropping constraint',cdict,val)
1111            skipList.append(cnum)
1112            newHolds += [i for i in holdList if i not in newHolds]
1113        elif len(dropList) > 0: # mix of valid and problematic items, drop problem vars, but keep rest
1114            if GSASIIpath.GetConfigValue('debug'): 
1115                msg = ''
1116                for v in dropList:
1117                    if msg: msg += ' ,'
1118                    msg += v
1119                warn('removing: '+msg,cdict,val)
1120            value = fixedList[cnum]
1121            for var in dropList:   # do cleanup
1122                # TODO: evaluate expressions in constraint multipliers here? I am assuming no
1123                # as I think G2mv.EvaluateMultipliers does that already. For now assuming
1124                # otherwise; note SubfromParmDict
1125                if ':dAx:' in var or ':dAy:' in var or ':dAz:' in var: # undefined atoms can be ignored
1126                    pass
1127                elif cdict[var] != 0: 
1128                    value = float(value) - cdict[var]*parmDict[var]
1129                del cdict[var]
1130            if float(value) != float(fixedList[cnum]): fixedList[cnum] = str(np.round(value,12))
1131            if GSASIIpath.GetConfigValue('debug'):
1132                warn('revised as: '+_FormatConstraint(constrDict[cnum],fixedList[cnum]))
1133    for i in list(range(len(constrDict)-1,-1,-1)): # remove the dropped constraints
1134        if i in skipList:
1135            del constrDict[i]
1136            del fixedList[i]
1137           
1138    if warning: warning += '\n'
1139    warning += warninfo['msg']
1140           
1141    groups,parmlist = GroupConstraints(constrDict)
1142
1143    # now process each group and create the relations that are needed to form
1144    # a non-singular square matrix
1145    # Now check that all parameters are varied (probably do not need to do this
1146    # any more). For constraint equations, if all are varied, set VaryFree to True
1147    # and all newly created relationships will be varied. For NewVar constraints,
1148    # vary if the vary flag was set.
1149    for group,varlist in zip(groups,parmlist):
1150        if len(varlist) < len(group): # too many relationships -- no can do
1151            if errmsg: errmsg += '\n'
1152            errmsg += "Over-constrained input. "
1153            errmsg += "There are more constraints (" + str(len(group))
1154            errmsg += ") than parameters (" + str(len(varlist)) + ")\nin these constraints:"
1155            for rel in group:
1156                errmsg += '\n\t'+ _FormatConstraint(constrDict[rel],fixedList[rel])
1157            groupErrors += varlist
1158            continue # go on to next group
1159
1160        VaryFree = False
1161        for rel in group:
1162            varied = 0
1163            for var in VarKeys(constrDict[rel]):
1164                #var = translateTable.get(var,var) # replace wildcards
1165                #if parmDict is not None and var not in parmDict:
1166                if var in varyList: varied += 1
1167            if fixedList[rel] is not None and varied > 0:
1168                VaryFree = True
1169               
1170        # fill in additional degrees of freedom
1171        try:
1172            arr = _FillArray(group,constrDict,varlist)
1173            _RowEchelon(len(group),arr,varlist)
1174        except:
1175            if errmsg: errmsg += '\n'
1176            errmsg += "\nSingular input. "
1177            errmsg += "There are internal inconsistencies in these constraints:"
1178            for rel in group:
1179                errmsg += '\n\t' + _FormatConstraint(constrDict[rel],fixedList[rel])
1180            groupErrors += varlist
1181            continue
1182
1183        try:
1184            constrArr = _FillArray(group,constrDict,varlist,FillDiagonals=True)
1185            GramSchmidtOrtho(constrArr,len(group))
1186        except:
1187            if errmsg: errmsg += '\n'
1188            errmsg += "\nUnexpected singularity with constraints group (in Gram-Schmidt)"
1189            for rel in group:
1190                errmsg += '\n\t' + _FormatConstraint(constrDict[rel],fixedList[rel])
1191            groupErrors += varlist
1192            continue
1193       
1194        try:
1195            np.linalg.inv(constrArr)
1196        except:
1197            if errmsg: errmsg += '\n'
1198            errmsg += "\nSingular input. "
1199            errmsg += "The following constraints are not "
1200            errmsg += "linearly independent\nor do not "
1201            errmsg += "allow for generation of a non-singular set.\n"
1202            errmsg += 'This is unexpected. Please report this (toby@anl.gov)'
1203            for rel in group:
1204                errmsg += '\n\t' + _FormatConstraint(constrDict[rel],fixedList[rel])
1205            groupErrors += varlist
1206            continue
1207       
1208        mapvar = []
1209        group = group[:]
1210        # scan through all generated and input relationships, we need to add to the varied list
1211        # all the new parameters where VaryFree has been set or where a New Var is varied.
1212        #
1213        # If a group does not contain any fixed values (constraint equations)
1214        # and nothing in the group is varied, drop this group, so that the
1215        # dependent parameters can be refined individually.
1216        unused = True
1217        for i in range(len(varlist)):
1218            if len(group) > 0: # get the original equation reference
1219                rel = group.pop(0)
1220                fixedval = fixedList[rel]
1221                varyflag = constrDict[rel].get('_vary',False)
1222                varname = constrDict[rel].get('_name','')
1223            else: # this relationship has been generated
1224                varyflag = False
1225                varname = ''
1226                fixedval = None
1227            if fixedval is None: # this is a new parameter, not a constraint
1228                if not varname:
1229                    varname = paramPrefix + str(consNum) # no assigned name, create one
1230                    consNum += 1
1231                mapvar.append(varname)
1232                #genVarLookup[varname] = varlist # save list of parameters related to this new var
1233                # vary the new relationship if it is a degree of freedom in
1234                # a set of contraint equations or if a New Var is flagged to be varied.
1235                if VaryFree or varyflag: 
1236                    unused = False
1237                    varyList.append(varname)
1238                    # fix (prevent varying) of all the parameters inside the constraint group
1239                    # (dependent vars)
1240                    for var in varlist:
1241                        if var in varyList: varyList.remove(var)
1242            else:
1243                unused = False
1244                mapvar.append(float(fixedval))
1245        if unused: # skip over constraints that don't matter (w/o fixed value or any refined parameters)
1246            if GSASIIpath.GetConfigValue('debug'): 
1247                print('Unexpected: Constraint ignored (all parameters unvaried)')
1248                print ('   '+_FormatConstraint(constrDict[rel],fixedList[rel]))
1249            continue 
1250        #dependentParmList.append([translateTable.get(var,var) for var in varlist])
1251        dependentParmList.append(varlist)
1252        arrayList.append(constrArr)
1253        invarrayList.append(np.linalg.inv(constrArr))
1254        indParmList.append(mapvar)
1255        symGenList.append(False)
1256    if errmsg and SeqHist is not None:
1257        print (' *** ERROR in constraint definitions! ***')
1258        print (errmsg)
1259        if warning:
1260            print (' also note warnings in constraint processing:')
1261            print (warning)
1262        raise ConstraintException
1263    elif errmsg:
1264        return errmsg,warning,None,None
1265
1266    # Make list of dependent and independent variables (after possibly dropping unused vars)
1267    global dependentVars
1268    global independentVars
1269    dependentVars = []
1270    independentVars = []
1271    for varlist,mapvars in zip(dependentParmList,indParmList):  # process all constraints
1272        for mv in mapvars:
1273            if type(mv) is float: continue
1274            if mv not in independentVars: independentVars.append(mv)
1275        for mv in varlist:
1276            if mv not in dependentVars: dependentVars.append(mv)
1277    saveVaryList = copy.copy(varyList)
1278
1279    # if equivMoved:
1280    #     print(60*'=')
1281    #     print('Constraints were reclassified to avoid conflicts, as below:')
1282    #     print(mvMsg)
1283    #     print('New constraints are:')
1284    #     print (VarRemapShow(varyList,True))
1285    #     print(60*'=')
1286    return errmsg,warning,groups,parmlist # saved for sequential fits
1287   
1288def CheckEquivalences(constrDict,varyList,fixedList,parmDict=None):
1289    '''Process equivalence constraints, looking for conflicts such as
1290    where a parameter is used in both an equivalence and a constraint expression
1291    or where chaining is done (A->B and B->C).
1292
1293    Removes equivalences or parameters from equivalences or converts equivalences to
1294    constraints as described for :ref:`Equivalence Checking and Reorganization <CheckEquivalences>`.
1295
1296    :param dict constrDict: a list of dicts defining relationships/constraints
1297    :param list varyList: list of varied parameters (defined during refinements only)
1298    :param list fixedList: a list of values specifying a fixed value for each
1299       dict in constrDict. Values are either strings that can be converted to
1300       floats or ``None`` if the constraint defines a new parameter rather
1301       than a constant.
1302    :param dict parmDict: a dict containing defined parameters and their values. Used to find
1303       equivalences where a parameter is has been removed from a refinement.
1304
1305    :returns: warning messages about changes that need to be made to equivalences
1306    '''
1307   
1308    warninfo = {'msg':'', 'shown':-1}
1309    def warn(msg,cnum=None):
1310        if cnum is not None and cnum != warninfo['shown']:
1311            warninfo['shown'] = cnum
1312            if warninfo['msg']: warninfo['msg'] += '\n'
1313            warninfo['msg'] += '\nProblem with equivalence: ' + _showEquiv(
1314                dependentParmList[cnum],indParmList[cnum],invarrayList[cnum])
1315        if warninfo['msg']: warninfo['msg'] += '\n'
1316        warninfo['msg'] += '  ' + msg
1317
1318    global holdParmList # parameters set as "Hold"
1319    global depVarList # parameters used in equivalences as dependent parameters
1320    global indepVarList # parameters used in equivalences as independent parameters
1321    global constrVarList  # parameters used in other constraints
1322             
1323    # lists of parameters used for error reporting
1324    global undefinedVars # parameters that are used in equivalences but are not defined
1325    global newHolds # additional parameters that should be held
1326    global convVarList # parameters in equivalences that will be converted to constraints
1327    convVarList = [] # parameters in equivalences to be made into constraints
1328    global multdepVarList
1329    multdepVarList = [] # list of dependent parameters used in more than one equivalence
1330
1331    # local vars
1332    dropVarList = []  # parameters that can be removed from equivalences
1333    removeList = []  # equivalences that are not needed
1334    convertList = [] # equivalences that should be converted to "Const" constraints
1335
1336    # process equivalences: make a list of dependent and independent vars
1337    #    and check for repeated uses (repetition of a parameter as an
1338    #    independent var is OK)
1339    # look for parameters in equivalences that are used more than once as dependent parameters
1340    seenOnce = []
1341    for cnum,(varlist,multarr) in enumerate(zip(dependentParmList,arrayList)):
1342        if multarr is not None: continue # equivalences only
1343        for v in varlist:
1344            if v not in seenOnce:
1345                seenOnce.append(v)
1346            elif v not in multdepVarList:
1347                multdepVarList.append(v)               
1348
1349    # scan through equivalences looking for other "dual uses". Stop when no new ones are found
1350    changed = True
1351    count = 0
1352    while changed:
1353        changed = False
1354        count += 1
1355        if count > 1000:
1356            raise ConstraintException("Too many loops in CheckEquivalences")
1357       
1358        # look for repeated dependent vars
1359        convVarList = [] # parameters in equivalences to be made into constraints
1360        for cnum,(varlist,mapvars,multarr,invmultarr) in enumerate(zip(
1361            dependentParmList,indParmList,arrayList,invarrayList)):
1362            if multarr is not None: continue # equivalences only
1363            if cnum in convertList:
1364                convVarList += [v for v in mapvars+varlist if v not in convVarList and type(v) is not float]
1365                continue
1366           
1367            # identify equivalences that need to be converted to constraints.
1368            #  Where parameters:
1369            #     are used both in equivalences & constraints,
1370            #     are used as dependent multiple times or
1371            #     where are used as both dependent and independent (chained)
1372            msg = False
1373            for v in mapvars:
1374                if v in constrVarList+convVarList:
1375                    changed = True
1376                    msg = True
1377                    warn("Independent parameter "+str(v)+' used in constraint',cnum)
1378                    if cnum not in convertList: convertList.append(cnum)
1379            for v in varlist:
1380                if v in multdepVarList:
1381                    changed = True
1382                    msg = True
1383                    warn("Dependent parameter "+str(v)+' repeated',cnum)
1384                    if cnum not in convertList: convertList.append(cnum)
1385                elif v in indepVarList:
1386                    changed = True
1387                    msg = True
1388                    warn("Dependent parameter "+str(v)+' used elsewhere as independent',cnum)
1389                    if cnum not in convertList: convertList.append(cnum)
1390                elif v in constrVarList+convVarList:
1391                    changed = True
1392                    msg = True
1393                    warn("Dependent parameter "+str(v)+' used in constraint',cnum)
1394                    if cnum not in convertList: convertList.append(cnum)
1395            if msg:
1396                warn('Converting to "Constr"',cnum)
1397           
1398    global unvariedParmsList
1399    unvariedParmsList = []  # parameters in equivalences that are not varied
1400    # scan equivalences: look for holds
1401    for cnum,(varlist,mapvars,multarr,invmultarr) in enumerate(zip(
1402        dependentParmList,indParmList,arrayList,invarrayList)):
1403        if multarr is not None: continue # not an equivalence
1404        if cnum in convertList: continue
1405
1406        # look for holds
1407        gotHold = False
1408        holdList = []
1409        for v in varlist+mapvars:
1410            if v in holdParmList:
1411                gotHold = True
1412            elif type(v) is not float:
1413                holdList.append(v)
1414        if gotHold:
1415            if holdList:
1416                msg = 'Some parameters set as "Hold"; setting remainder as "Hold": '
1417                for i,var in enumerate(holdList):
1418                    if i != 0: msg += ", "
1419                    msg += var
1420                newHolds += [i for i in holdList if i not in newHolds]
1421                msg += ". "
1422            else:
1423                msg = 'All parameters set as "Hold". '
1424            msg += " Ignoring equivalence"
1425            warn(msg,cnum)
1426            removeList.append(cnum)
1427            continue
1428       
1429        # look for unvaried parameters
1430        gotVary = False
1431        gotNotVary = False
1432        holdList = []
1433        for v in varlist+mapvars:
1434            if v in varyList:
1435                gotVary = True
1436                holdList.append(v)
1437            elif type(v) is not float:
1438                gotNotVary = True
1439                if v not in unvariedParmsList: unvariedParmsList.append(v)
1440        if gotNotVary:  # at least some unvaried parameters
1441            if gotVary:  # mix of varied and unvaried parameters
1442                msg = 'Some parameters not varied; setting remainder as "Hold": '
1443                for i,var in enumerate(holdList):
1444                    if i != 0: msg += ", "
1445                    msg += var
1446                newHolds += [i for i in holdList if i not in newHolds]
1447                msg += ". "
1448            else:
1449                msg = 'No parameters varied. '               
1450            msg += " Ignoring equivalence"
1451            warn(msg,cnum)
1452            removeList.append(cnum)
1453            continue
1454           
1455        # look for undefined or zero multipliers
1456        holdList = []
1457        drop = 0
1458        for v,m in zip(varlist,invmultarr):   
1459            if parmDict is not None and v not in parmDict:
1460                if v not in undefinedVars: undefinedVars.append(v)
1461                if v not in dropVarList: dropVarList.append(v)
1462                drop += 1
1463            elif m == 0:
1464                warn("Parameter "+str(v)+" has a zero multiplier, dropping",cnum)
1465                if v not in dropVarList: dropVarList.append(v)
1466                drop += 1
1467            else:
1468                holdList.append(v)
1469        if drop == len(varlist):
1470            warn("No dependent parameters defined, ignoring equivalence",cnum)
1471            removeList.append(cnum)
1472            continue
1473        for mv in mapvars:
1474            if type(mv) is float: continue
1475            if parmDict is not None and mv not in parmDict:
1476                # independent parameter is undefined, but some dependent parameters are defined
1477                # hold them
1478                if mv not in undefinedVars: undefinedVars.append(mv)
1479                msg = "Parameter(s) "+str(mv)
1480                for v in varlist:
1481                    if v in dropVarList:
1482                        msg += ', ' + v
1483                msg += "not defined in this refinement\n"
1484                msg = "Setting holds for: "
1485                for i,var in enumerate(holdList):
1486                    if i != 0: msg += ", "
1487                    msg += var
1488                warn(msg,cnum)
1489                drop += 1
1490        if drop: # independent var and at least one dependent variable is defined
1491            msg = "Dropping undefined parameter(s) "
1492            i = 0
1493            for v in varlist:
1494                if v in dropVarList:
1495                    if i != 0: msg += ', '
1496                    i += 1
1497                    msg += v
1498            warn(msg,cnum)
1499            msg = "Some parameters not defined. Setting holds for: "
1500            for i,var in enumerate(holdList):
1501                if i != 0: msg += ", "
1502                msg += var
1503            warn(msg,cnum)
1504            newHolds += [i for i in holdList if i not in newHolds]
1505   
1506    # Convert equivalences where noted
1507    for cnum,varlist in enumerate(dependentParmList):
1508        if cnum not in convertList: continue
1509        indvar = indParmList[cnum][0]
1510        # msg = '\nChanging equivalence:\n    ' + _showEquiv(
1511        #     dependentParmList[cnum],indParmList[cnum],invarrayList[cnum])
1512        for dep,mult in zip(dependentParmList[cnum],invarrayList[cnum]):
1513            constrDict += [{indvar:-1.,dep:mult[0]}]
1514            fixedList += ['0.0']
1515        #msg += '\n  to constraint(s):'
1516        #msg += '\n    ' + _FormatConstraint(constrDict[-1],fixedList[-1])
1517        removeList.append(cnum)
1518    # Drop equivalences where noted
1519    if removeList:
1520        for i in sorted(set(removeList),reverse=True):
1521            del dependentParmList[i],indParmList[i],arrayList[i],invarrayList[i],symGenList[i]
1522    # Drop variables from remaining equivalences
1523    for cnum,varlist in enumerate(dependentParmList):
1524        for j,v in enumerate(varlist):
1525            drop = []
1526            if v in dropVarList:
1527                drop.append(j)
1528        if drop:
1529            for j in sorted(drop,reverse=True):
1530                del indParmList[cnum][j]
1531    return warninfo['msg']
1532
1533def ProcessConstraints(constList,seqmode='use-all',seqhst=None):
1534    """Interpret the constraints in the constList input into a dictionary, etc.
1535    All :class:`GSASIIobj.G2VarObj` objects are mapped to the appropriate
1536    phase/hist/atoms based on the object internals (random Ids). If this can't be
1537    done (if a phase has been deleted, etc.), the variable is ignored.
1538    If the constraint cannot be used due to too many dropped variables,
1539    it is counted as ignored. In the case of sequential refinements,
1540    the current histogram number is substituted for a histogram number of "*".
1541
1542    NB: this processing does not include symmetry imposed constraints
1543   
1544    :param list constList: a list of lists where each item in the outer list
1545      specifies a constraint of some form, as described in the :mod:`GSASIIobj`
1546      :ref:`Constraint definitions <Constraint_definitions_table>`.
1547    :param str seqmode: one of 'use-all', 'wildcards-only' or 'auto-wildcard'.
1548       When seqmode=='wildcards-only' then any constraint with a numerical
1549       constraint number is skipped. With seqmode=='auto-wildcard',
1550       any non-null constraint number is set to the selected histogram.
1551    :param int seqhst: number for current histogram (used for
1552      'wildcards-only' or 'auto-wildcard' only). Should be None for
1553      non-sequential fits.
1554
1555    :returns:  a tuple of (constrDict,fixedList,ignored) where:
1556     
1557      * constrDict (list of dicts) contains the constraint relationships
1558      * fixedList (list) contains the fixed values for each type
1559        of constraint.
1560      * ignored (int) counts the number of invalid constraint items
1561        (should always be zero!)
1562    """
1563    constrDict = []
1564    fixedList = []
1565    ignored = 0
1566    namedVarList = []
1567    for constr in constList:
1568        terms = copy.deepcopy(constr[:-3])
1569        if seqmode == 'wildcards-only' and seqhst is not None:
1570            skip = False
1571            for term in terms:
1572                if term[1].histogram == '*':
1573                    term[1] = term[1].varname(seqhst)
1574                elif term[1].histogram:
1575                    skip = True
1576            if skip: continue
1577        elif seqmode == 'auto-wildcard' and seqhst is not None:
1578            for term in terms:
1579                term[1] = term[1].varname(seqhst)
1580        else:
1581            for term in terms:
1582                if term[1].histogram == '*' and seqhst is not None:
1583                    term[1] = term[1].varname(seqhst)
1584                else:
1585                    term[1] = term[1].varname()
1586               
1587        if constr[-1] == 'h':
1588            # process a hold
1589            var = str(constr[0][1])
1590            if '?' not in var:
1591                StoreHold(var)
1592            else:
1593                ignored += 1
1594        elif constr[-1] == 'f':
1595            # process a new variable
1596            fixedList.append(None)
1597            D = {}
1598            varyFlag = constr[-2]
1599            varname = constr[-3]
1600            for term in terms:
1601                var = str(term[1])
1602                if '?' not in var:
1603                    D[var] = term[0]
1604            if len(D) > 1:
1605                # add extra dict terms for input variable name and vary flag
1606                if varname is not None:
1607                    varname = str(varname) # in case this is a G2VarObj
1608                    if varname.startswith(':'):
1609                        D['_name'] = varname
1610                    else:
1611                        D['_name'] = '::nv-' + varname
1612                    D['_name'] = G2obj.MakeUniqueLabel(D['_name'],namedVarList)
1613                D['_vary'] = varyFlag == True # force to bool
1614                constrDict.append(D)
1615            else:
1616                ignored += 1
1617            #constFlag[-1] = ['Vary']
1618        elif constr[-1] == 'c': 
1619            # process a contraint relationship
1620            D = {}
1621            for term in terms:
1622                var = str(term[1])
1623                if '?' not in var:
1624                    D[var] = term[0]
1625            if len(D) >= 1:
1626                fixedList.append(str(constr[-3]))
1627                constrDict.append(D)
1628            else:
1629                ignored += 1
1630        elif constr[-1] == 'e':
1631            # process an equivalence
1632            firstmult = None
1633            eqlist = []
1634            for term in terms:
1635                if term[0] == 0: term[0] = 1.0
1636                var = str(term[1])
1637                if '?' in var: continue
1638                if firstmult is None:
1639                    firstmult = term[0]
1640                    firstvar = var
1641                else:
1642                    eqlist.append([var,firstmult/term[0]])
1643            if len(eqlist) > 0:
1644                StoreEquivalence(firstvar,eqlist,False)
1645            else:
1646                ignored += 1
1647        else:
1648            ignored += 1
1649    return constrDict,fixedList,ignored
1650
1651def StoreHold(var,symGen=False):
1652    '''Takes a variable name and prepares it to be removed from the
1653    refined variables.
1654
1655    Called with user-supplied constraints by :func:`ProcessConstraints`.
1656    At present symGen is not used, but could be set up to track Holds generated
1657    by symmetry.
1658    '''
1659    global holdParmList
1660    if var not in holdParmList:
1661        holdParmList.append(var)
1662
1663   
1664def StoreEquivalence(independentVar,dependentList,symGen=True):
1665    '''Takes a list of dependent parameter(s) and stores their
1666    relationship to a single independent parameter (independentVar).
1667
1668    Called with user-supplied constraints by :func:`ProcessConstraints`,
1669    with Pawley constraints from :func:`GSASIIstrIO.GetPawleyConstr`,
1670    with Unit Cell constraints from :func:`GSASIIstrIO.cellVary`
1671    with symmetry-generated atom constraints from :func:`GSASIIstrIO.GetPhaseData`
1672
1673      There is no harm in using StoreEquivalence with the same independent variable::
1674
1675       StoreEquivalence('x',('y',))
1676       StoreEquivalence('x',('z',))
1677
1678      but the same outcome can be obtained with a single call::
1679
1680       StoreEquivalence('x',('y','z'))
1681
1682      The latter will run more efficiently.
1683
1684
1685      Note that mixing independent and dependent variables, such as::
1686
1687        StoreEquivalence('x',('y',))
1688        StoreEquivalence('y',('z',))
1689
1690      is a poor choice. The module will attempt to fix this by transforming the equivalence to a
1691      "Const" constraint.
1692
1693    :param str independentVar: name of master parameter that will be used to determine the value
1694      to set the dependent variables
1695
1696    :param list dependentList: a list of parameters that will set from
1697         independentVar. Each item in the list can be a string with the parameter
1698         name or a tuple containing a name and multiplier:
1699         ``['::parm1',('::parm2',.5),]``
1700
1701    '''
1702   
1703    global dependentParmList,arrayList,invarrayList,indParmList,symGenList
1704    mapList = []
1705    multlist = []
1706    allfloat = True
1707    for var in dependentList:
1708        if isinstance(var, str):
1709            mult = 1.0
1710        elif len(var) == 2:
1711            var,mult = var
1712        else:
1713            raise Exception("Cannot parse "+repr(var) + " as var or (var,multiplier)")
1714        mapList.append(var)
1715        try:           
1716            multlist.append(tuple((float(mult),)))
1717        except:
1718            allfloat = False
1719            multlist.append(tuple((mult,)))
1720    # added relationships to stored values
1721    arrayList.append(None)
1722    if allfloat:
1723        invarrayList.append(np.array(multlist))
1724    else:
1725        invarrayList.append(multlist)
1726    indParmList.append(list((independentVar,)))
1727    dependentParmList.append(mapList)
1728    symGenList.append(symGen)
1729    return
1730
1731def SubfromParmDict(s,prmDict):
1732    for key in prmDict:
1733        if key in s:
1734            s = s.replace(key,str(prmDict[key]))
1735    return eval(s)
1736
1737def EvaluateMultipliers(constList,*dicts):
1738    '''Convert multipliers for constraints and equivalences that are specified
1739    as strings into values. The strings can specify values in the parameter dicts as
1740    well as normal Python functions, such as "2*np.cos(0::Ax:2/2.)"
1741   
1742    :param list constList: a list of dicts containing constraint expressions
1743    :param \*dicts: one or more dicts containing GSAS-II parameters and their values
1744       can be specified
1745    :returns: an empty string if there were no errors, or an error message listing
1746       the strings that could not be converted.
1747    '''
1748    prmDict = {}
1749    for d in dicts: prmDict.update(d) # combine all passed parameter dicts
1750    problemList = ""
1751    # loop through multipliers in contraint expressions
1752    for const in constList:
1753        for key in const:
1754            if key.startswith('_'): continue
1755            try: # is this already a float, etc?
1756                1+const[key]
1757                continue
1758            except:
1759                pass
1760            try:
1761                newval = SubfromParmDict(const[key][:],prmDict)
1762                if GSASIIpath.GetConfigValue('debug'):
1763                    print('Changing ',const[key],'to',newval)
1764                const[key] = newval               
1765            except:
1766                if problemList: problemList += ", "
1767                problemList += const[key]
1768    # loop through multipliers in equivalences
1769    global arrayList,invarrayList
1770    for i,(a,valList) in enumerate(zip(arrayList,invarrayList)):
1771        if a is not None: continue # ignore if not equiv
1772        try:
1773            valList.shape
1774            continue # ignore if already a numpy array
1775        except:
1776            pass
1777        repList = []
1778        for v in valList:
1779            try:
1780                1+v[0]
1781                repList.append(tuple((v[0],)))
1782                continue
1783            except:
1784                pass
1785            try:
1786                newval = SubfromParmDict(v[0][:],prmDict)
1787                if GSASIIpath.GetConfigValue('debug'):
1788                    print('Changing ',v[0],'to',newval)
1789                repList.append(tuple((newval,)))
1790            except:
1791                if problemList: problemList += ", "
1792                problemList += v[0]
1793                repList.append(tuple(('error',)))
1794        invarrayList[i] = np.array(repList)
1795    return problemList
1796
1797def GetDependentVars():
1798    '''Return a list of dependent variables: e.g. parameters that are
1799    constrained in terms of other parameters
1800
1801    :returns: a list of parameter names
1802
1803    '''
1804    global dependentVars
1805    return dependentVars
1806
1807def GetIndependentVars():
1808    '''Return a list of independent variables: e.g. parameters that are
1809    slaved to other parameters by constraints
1810
1811    :returns: a list of parameter names
1812
1813    '''
1814    global independentVars
1815    return independentVars
1816
1817def PrintIndependentVars(parmDict,varyList,sigDict,PrintAll=False,pFile=None):
1818    '''Print the values and uncertainties on the independent parameters'''
1819    global dependentParmList,arrayList,invarrayList,indParmList
1820    printlist = []
1821    mvs = GetIndependentVars()
1822    for i,name in enumerate(mvs):
1823        if PrintAll or name in varyList:
1824            sig = sigDict.get(name)
1825            printlist.append([name,parmDict[name],sig])
1826    if len(printlist) == 0: return
1827    s1 = ''
1828    s2 = ''
1829    s3 = ''
1830    pFile.write(130*'-'+'\n')
1831    pFile.write("Parameters generated by constraints\n")
1832    printlist.append(3*[None])
1833    for name,val,esd in printlist:
1834        if len(s1) > 120 or name is None:
1835            pFile.write(''+'\n')
1836            pFile.write(s1+'\n')
1837            pFile.write(s2+'\n')
1838            pFile.write(s3+'\n')
1839            s1 = ''
1840            if name is None: break
1841        if s1 == "":
1842            s1 = ' name  :'
1843            s2 = ' value :'
1844            s3 = ' sig   :'
1845        s1 += '%15s' % (name)
1846        s2 += '%15.5f' % (val)
1847        if esd is None:
1848            s3 += '%15s' % ('n/a')
1849        else:
1850            s3 += '%15.5f' % (esd)
1851
1852def getConstrError(constrLst,seqmode,seqhst):
1853    '''This is used to display error messages for constraints and
1854    equivalence relations
1855
1856    :parm list constrLst: a single constraint or equivalence as saved
1857      in the data tree
1858      (see :ref:`constraint definitions <Constraint_definitions_table>`).
1859    :param str seqmode: one of 'use-all', 'wildcards-only' or 'auto-wildcard'
1860    :param int seqhst: number for current histogram (used for
1861      'wildcards-only' or 'auto-wildcard' only). Should be None for
1862      non-sequential fits.
1863
1864    :returns: error, msg where error (bool) is True if the
1865      constraint/equivalence creates an error, msg (str) can be a warning
1866      or an error
1867    '''
1868    msg = ''
1869    note = ''
1870    terms = copy.deepcopy(constrLst[:-3])
1871    if seqmode == 'wildcards-only' and seqhst is not None:
1872        if constrLst[-1] == 'e':
1873            msg = 'equivalence'
1874        else:
1875            msg = 'constraint'
1876        for term in terms:
1877            if term[1].histogram == '*':
1878                term[1] = term[1].varname(seqhst)
1879            elif term[1].histogram:
1880                return False,"Ignoring non-wildcard "+msg, "Ignore"
1881    elif seqmode == 'auto-wildcard' and seqhst is not None:
1882        for term in terms:
1883            term[1] = term[1].varname(seqhst)
1884    else:
1885        for term in terms:
1886            if term[1].histogram == '*':
1887                if seqhst is None:
1888                    msg = "Parameter "+str(terms[0][1])+" contains a wildcard, which are used only sequential refinements. Constraint ignored."
1889                    return False,msg, "Ignored"
1890                else:
1891                    term[1] = term[1].varname(seqhst)
1892            else:
1893                term[1] = term[1].varname()
1894    if constrLst[-1] == 'e':
1895        # conflicting uses
1896        if terms[0][1] in constrVarList+convVarList:
1897            msg = "Parameter "+str(terms[0][1])+" used in constraint. To be recast as constraint."
1898            return False,msg, "Recast as constraint."
1899        varList = []
1900        for m,v in terms:
1901            if v in constrVarList+convVarList:
1902                varList.append(str(v))
1903        if varList:
1904            msg = "Parameter(s) used in constraint: "
1905            for i,v in enumerate(varList):
1906                if i != 0: msg += ', '
1907                msg += v
1908            return False,msg,"Recast as constraint."
1909        varList = []
1910        for m,v in terms[1:]:
1911            if v in indepVarList:
1912                varList.append(str(v))
1913        if varList:
1914            msg = "Parameter(s) used elsewhere as independent: "
1915            for i,v in enumerate(varList):
1916                if i != 0: msg += ', '
1917                msg += v
1918            return False,msg,"Recast as constraint."
1919        varList = []
1920        for m,v in terms[1:]:
1921            if v in multdepVarList:
1922                varList.append(str(v))
1923        if varList:
1924            msg += "Parameter(s) repeated as dependent: "
1925            for i,v in enumerate(varList):
1926                if i != 0: msg += ', '
1927                msg += v
1928            return False,msg,"Recast as constraint."
1929
1930        # zero multiplier
1931        varList = []
1932        valid = 0
1933        for m,v in terms:
1934            if m == 0:
1935                varList.append(str(v))
1936            else:
1937                valid += 1
1938        if varList and valid > 1:
1939            msg += "Parameter(s) with zero multipliers: "
1940            for i,v in enumerate(varList):
1941                if i != 0: msg += ', '
1942                msg += v
1943            msg += " will be ignored"
1944        elif varList:
1945            msg += "Parameter(s) with zero multipliers:"
1946            for i,v in enumerate(varList):
1947                if i != 0: msg += ', '
1948                msg += v
1949            return False,msg,"Equivalence Ignored."
1950       
1951        # hold parameters
1952        s = ''
1953        for m,v in terms:
1954            if v in holdParmList:
1955                if s: s += ', '
1956                s += str(v)
1957        if s:
1958            if msg: msg += '; '
1959            msg += "Parameters set as Hold: "+s
1960            return False,msg,'Has holds: Nothing varied.'
1961               
1962        # unrefined parameters
1963        gotVary = False
1964        gotNotVary = False
1965        s = ''
1966        for m,v in terms:
1967            if v in unvariedParmsList:
1968                gotNotVary = True
1969                if s: s += ', '
1970                s += str(v)
1971            else:
1972                gotVary = True
1973        if gotNotVary and gotVary:  # mix of varied and unvaried parameters
1974            if msg: msg += '. '
1975            msg += 'Unvaried parameter(s): '+s+"; remainder set Hold. All parameters fixed."
1976            return False,msg, 'Equivalence Ignored.'
1977        elif gotNotVary:
1978            if msg: msg += '. '
1979            msg += 'All parameters not varied. Equivalence Ignored.'
1980            return False,msg, 'Equivalence Ignored.'
1981           
1982        # undefined parameters
1983        undef = 0
1984        s = ''
1985        for m,v in terms[1:]:
1986            if v in undefinedVars:
1987                undef += 1
1988                if s: s += ', '
1989                s += str(v)
1990        if undef == len(terms[1:]):
1991            msg += 'Ignored: None of the dependent parameters are defined'
1992        elif terms[0][1] in undefinedVars:
1993            if s:
1994                s = terms[0][1] + ', ' + s
1995            else:
1996                s = terms[0][1]
1997            msg += 'Undefined parameter(s): '+s+'. Remainder will be fixed'
1998        elif undef:
1999            msg += 'Undefined parameter(s): '+s+' will be dropped'
2000    elif constrLst[-1] == 'h':
2001        v = terms[0][1]
2002        if v in undefinedVars: return False,"Parameter is undefined","Ignored"
2003        if v in unvariedParmsList: return False,"Parameter is not refined","Ignored"
2004    else:
2005        # check for post-grouping errors
2006        for m,v in terms:
2007            if v in groupErrors:
2008                return True,'Constraint singularity: see error listing','Singular'
2009        zeroList = []
2010        tobeHold = []
2011        undef = []
2012        unvar = []
2013        hold = []
2014        for m,v in terms: # check for zero multiplier, undefined, unvaried or hold
2015            if m == 0:
2016                zeroList.append(str(v))
2017            elif v in undefinedVars:
2018                undef.append(str(v))
2019            elif v in unvariedParmsList:
2020                unvar.append(str(v))
2021            elif v in holdParmList:
2022                hold.append(str(v))
2023            else:
2024                tobeHold.append(str(v))
2025        s = ''
2026        for v in zeroList:
2027            if s: s += ', '
2028            s += str(v)
2029        if s:
2030            if msg: msg += '; '
2031            msg += "Parameter(s) with zero multipliers: "+s
2032        s = ''
2033        for v in undef:
2034            if s: s += ', '
2035            s += str(v)
2036        if s:
2037            if msg: msg += '; '
2038            msg += "Undefined parameter(s): "+s
2039        s = ''
2040        for v in unvar:
2041            if s: s += ', '
2042            s += str(v)
2043        if s:
2044            if msg: msg += '; '
2045            msg += "Unrefined parameter(s): "+s
2046        s = ''
2047        for v in hold:
2048            if s: s += ', '
2049            s += str(v)
2050        if s:
2051            if msg: msg += '; '
2052            msg += '"Hold" parameter(s): '+s
2053        if undef and tobeHold:
2054            s = ''
2055            for v in tobeHold:
2056                if s: s += ', '
2057                s += str(v)
2058            if msg: msg += '; '
2059            msg += "Adding Holds on "+s+"; Constraint Ignored."
2060            note = 'Ignored'
2061        elif undef or (len(tobeHold) == 0 and (zeroList or unvar or hold)):
2062            if msg: msg += '; '
2063            msg += "Constraint Ignored."
2064            note = 'Ignored'
2065        elif len(tobeHold) == 1 and (zeroList or unvar or hold):
2066            if msg: msg += '; '
2067            msg += "One parameter is retained; converted to fixed value."
2068            note = 'Converted'
2069        elif zeroList or unvar or hold:
2070            if msg: msg += ': '
2071            msg += 'Will be dropped. Constraint retained.'
2072            note = 'Parameters removed'
2073    return False,msg,note
2074       
2075def ComputeDepESD(covMatrix,varyList,parmDict):
2076    '''Compute uncertainties for dependent parameters from independent ones
2077    returns a dictionary containing the esd values for dependent parameters
2078    '''
2079    sigmaDict = {}
2080    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
2081        #if invmultarr is None: continue # probably not needed
2082#        try:
2083#            valuelist = [parmDict[var] for var in mapvars]
2084#        except KeyError:
2085#            continue
2086        # get the v-covar matrix for independent parameters
2087        vcov = np.zeros((len(mapvars),len(mapvars)))
2088        for i1,name1 in enumerate(mapvars):
2089            if name1 not in varyList: continue
2090            iv1 = varyList.index(name1)
2091            for i2,name2 in enumerate(mapvars):
2092                if name2 not in varyList: continue
2093                iv2 = varyList.index(name2)
2094                vcov[i1][i2] = covMatrix[iv1][iv2]
2095        # vec is the vector that multiplies each of the independent values
2096        for v,vec in zip(varlist,invmultarr):
2097            sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec)))
2098    return sigmaDict
2099
2100def _FormatConstraint(RelDict,RelVal):
2101    '''Formats a Constraint or Function for use in a convenient way'''
2102    linelen = 65
2103    s = [""]
2104    for var,val in RelDict.items():
2105        if var.startswith('_'): continue
2106        if len(s[-1]) > linelen: s.append(' ')
2107        m = val
2108        if s[-1] != "" and m >= 0:
2109            s[-1] += ' + '
2110        elif s[-1] != "":
2111            s[-1] += ' - '
2112            m = abs(m)
2113        if m == 1:
2114            s[-1] += '%s '%var
2115        else:
2116            s[-1] += '%.3f*%s '%(m,var)
2117    if len(s[-1]) > linelen: s.append(' ')
2118    if RelVal is None:
2119        s[-1] += ' = New variable'
2120    else:
2121        s[-1] += ' = ' + str(RelVal)
2122    s1 = ''
2123    for s2 in s:
2124        if s1 != '': s1 += '\n\t'
2125        s1 += s2
2126    return s1
2127
2128def _showEquiv(varlist,mapvars,invmultarr,longmsg=False):
2129    '''Format an equivalence relationship
2130    note that
2131    varlist,           mapvars,     invmultarr
2132    are elements of
2133    dependentParmList, indParmList, invarrayList
2134    '''
2135    for i,mv in enumerate(mapvars):
2136        s1 = str(mv)
2137        if not longmsg:
2138            s1 += ' ==> '
2139        elif len(varlist) == 1:
2140            s1 += ' is equivalent to '
2141        else:
2142            s1 += ' is equivalent to parameters: '
2143        j = 0
2144        for v,m in zip(varlist,invmultarr):
2145            if debug: print ('v,m[0]: ',v,m[0])
2146            if len(s1.split('\n')[-1]) > 75: s1 += '\n        '
2147            if j > 0: s1 += ' & '
2148            j += 1
2149            s1 += str(v)
2150            if m != 1:
2151                s1 += " / " + str(m[0])
2152    return s1
2153
2154def VarRemapShow(varyList=None,inputOnly=False):
2155    '''List out the saved relationships. This should be done after the constraints have been
2156    defined using :func:`StoreEquivalence`, :func:`GroupConstraints` and :func:`GenerateConstraints`.
2157
2158    :returns: a string containing the details of the contraint relationships
2159    '''
2160    if varyList is None:
2161        varyList = saveVaryList
2162    s = ''
2163    if len(holdParmList) > 0:
2164        s += 'User-supplied Fixed Parameters:\n'
2165        for v in holdParmList:
2166            s += '    ' + v + '\n'
2167    if len(newHolds) > 0:
2168        s += 'Additional Fixed Parameters:\n'
2169        for v in newHolds:
2170            s += '    ' + v + '\n'
2171    userOut = ''
2172    symOut = ''
2173    consOut = ''
2174    varOut = ''
2175    global dependentParmList,arrayList,invarrayList,indParmList,symGenList
2176
2177    for varlist,mapvars,multarr,invmultarr,symFlag in zip(
2178        dependentParmList,indParmList,arrayList,invarrayList,symGenList):
2179        for i,mv in enumerate(mapvars):
2180            if multarr is None:
2181#                s1 = '  ' + str(mv) + ' is equivalent to parameter(s): '
2182                if len(varlist) == 1:
2183                    s1 = '   ' + str(mv) + ' is equivalent to '
2184                else:
2185                    s1 = '   ' + str(mv) + ' is equivalent to parameters: '
2186                j = 0
2187                for v,m in zip(varlist,invmultarr):
2188                    if debug: print ('v,m[0]: ',v,m[0])
2189                    if len(s1.split('\n')[-1]) > 75: s1 += '\n        '
2190                    if j > 0: s1 += ' & '
2191                    j += 1
2192                    s1 += str(v)
2193                    if m != 1:
2194                        s1 += " / " + str(m[0])
2195                if symFlag:
2196                    symOut += s1 + '\n'
2197                else:
2198                    userOut += s1 + '\n'
2199                continue
2200            if mv in varyList: 
2201                lineOut = {} (**Varied**) = '.format(mv)
2202            else:
2203                lineOut = {} = '.format(mv)
2204            j = 0 
2205            for (m,v) in zip(multarr[i,:],varlist):
2206                if m == 0: continue
2207                if m < 0:
2208                    lineOut += ' - '
2209                    m *= -1
2210                elif j != 0:
2211                    lineOut += ' + '
2212                j += 1
2213                if len(lineOut) > 60 and type(mv) is float:
2214                    consOut += lineOut
2215                    lineOut = '\n\t'
2216                elif len(lineOut) > 60:
2217                    varOut += lineOut
2218                    lineOut = '\n\t'
2219                if m == 1:
2220                    lineOut += '{}'.format(v)
2221                else:
2222                    lineOut += '({:.4g} * {})'.format(m,v)
2223            if type(mv) is not float:
2224                varOut += lineOut + '\n'
2225            else:
2226                consOut += lineOut + '\n'
2227    if userOut:
2228        s += '\nUser-supplied equivalences:\n' + userOut
2229    if consOut:
2230        s += '\nUser-supplied Constraints:\n' + consOut
2231    if varOut:
2232        s += '\nUser-input generated New Variables:\n' + varOut
2233    if symOut:
2234        s += '\nSymmetry-generated equivalences:\n' + symOut
2235    if not (userOut or consOut or varOut or symOut):
2236        return s + '\nNo constraints or equivalences in use'
2237    elif inputOnly:
2238        return s
2239       
2240    s += '\nInverse parameter mapping relations:\n'
2241    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
2242        for i,mv in enumerate(varlist):
2243            lineOut = {} = '.format(mv)
2244            j = 0 
2245            for m,v in zip(invmultarr[i,:],mapvars):
2246                if m == 0: continue
2247                if v == 0: continue
2248                if m < 0:
2249                    lineOut += ' - '
2250                    m *= -1
2251                elif j != 0:
2252                    lineOut += ' + '
2253                j += 1
2254                if len(lineOut) > 60:
2255                    s += lineOut
2256                    lineOut = '\n\t'
2257                if m == 1:
2258                    lineOut += '{}'.format(v)
2259                else:
2260                    try:
2261                        lineOut += '{:.4g}'.format(m*v)
2262                    except:
2263                        lineOut += '({:.4g} * {})'.format(m,v)
2264            if j == 0: lineOut += '0'
2265            s += lineOut + '\n'
2266    return s
2267
2268def GetSymEquiv():
2269    '''Return the automatically generated (equivalence) relationships.
2270
2271    :returns: a list of strings containing the details of the contraint relationships
2272    '''
2273    symout = []
2274    symerr = []
2275    symhelp = []
2276    global dependentParmList,arrayList,invarrayList,indParmList,symGenList
2277
2278    for varlist,mapvars,multarr,invmultarr,symFlag in zip(
2279        dependentParmList,indParmList,arrayList,invarrayList,symGenList):
2280        if not symFlag: continue
2281        for i,mv in enumerate(mapvars):
2282            cnstr = [[1,mv]]
2283            if multarr is None:
2284                s1 = ''
2285                s2 = ' = ' + str(mv)
2286                j = 0
2287                helptext = 'Variable {:} '.format(mv) + " ("+ G2obj.fmtVarDescr(mv) + ")"
2288                if len(varlist) == 1:
2289                    cnstr.append([invmultarr[0][0],varlist[0]])
2290                    # format the way Bob prefers
2291                    if invmultarr[0][0] == 1: 
2292                        s1 = str(varlist[0]) + ' = ' + str(mv)
2293                    else:
2294                        s1 = str(varlist[0]) + ' = ' + str(
2295                            invmultarr[0][0]) + ' * '+ str(mv)
2296                    s2 = ''
2297                   
2298                    m = 1./invmultarr[0][0]
2299                    var1 = str(varlist[0])
2300                    helptext += "\n\nis equivalent to "
2301                    if m == 1:
2302                        helptext += '\n  {:} '.format(var1) + " ("+ G2obj.fmtVarDescr(var1) + ")"
2303                    else:
2304                        helptext += '\n  {:3g} * {:} '.format(m,var1) + " ("+ G2obj.fmtVarDescr(var1) + ")"
2305                else:
2306                    helptext += "\n\nis equivalent to the following:"
2307                    for v,m in zip(varlist,invmultarr):
2308                        cnstr.append([m,v])
2309                        #if debug: print ('v,m[0]: ',v,m[0])
2310                        if len(s1.split('\n')[-1]) > 75: s1 += '\n        '
2311                        if j > 0: s1 += ' =  '
2312                        j += 1
2313                        s1 += str(v)
2314                        if m != 1:
2315                            s1 += " / " + str(m[0])
2316                            helptext += '\n  {:3g} * {:} '.format(m,v) + " ("+ G2obj.fmtVarDescr(v) + ")"
2317                        else:
2318                            helptext += '\n  {:} '.format(v) + " ("+ G2obj.fmtVarDescr(v) + ")"
2319                err,msg,note = getConstrError(cnstr+[None,None,'e'])
2320                symerr.append([msg,note])
2321                symout.append(s1+s2)
2322                symhelp.append(helptext)
2323            else:
2324                s = %s = ' % mv
2325                j = 0
2326                for m,v in zip(multarr[i,:],varlist):
2327                    if m == 0: continue
2328                    if j > 0: s += ' + '
2329                    j += 1
2330                    s += '(%s * %s)' % (m,v)
2331                print ('unexpected sym op='+s)
2332    return symout,symerr,symhelp
2333
2334def Dict2Deriv(varyList,derivDict,dMdv):
2335    '''Compute derivatives for Independent Parameters from the
2336    derivatives for the original parameters
2337
2338    :param list varyList: a list of parameters names that will be varied
2339
2340    :param dict derivDict: a dict containing derivatives for parameter values keyed by the
2341      parameter names.
2342
2343    :param list dMdv: a Jacobian, as a list of np.array containing derivatives for dependent
2344      parameter computed from derivDict
2345
2346    '''
2347    global dependentParmList,arrayList,invarrayList,indParmList,invarrayList
2348    for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
2349        for i,name in enumerate(mapvars):
2350            # grouped parameters: need to add in the derv. w/r
2351            # dependent variables to the independent ones
2352            if name not in varyList: continue # skip if independent var not varied
2353            if multarr is None:
2354                for v,m in zip(varlist,invmultarr):
2355                    if debug: print ('start dMdv',dMdv[varyList.index(name)])
2356                    if debug: print ('add derv',v,'/',m[0],'to derv',name,'add=',derivDict[v] / m[0])
2357                    if m == 0: continue
2358                    dMdv[varyList.index(name)] += derivDict[v] / m[0]
2359            else:
2360                for v,m in zip(varlist,multarr[i,:]):
2361                    if debug: print ('start dMdv',dMdv[varyList.index(name)])
2362                    if debug: print ('add derv',v,'*',m,'to derv',name,'add=',m * derivDict[v])
2363                    if m == 0: continue
2364                    dMdv[varyList.index(name)] += m * derivDict[v]
2365
2366def Map2Dict(parmDict,varyList):
2367    '''Updates the parameter dictionary and the varyList using the
2368    equivalence and constraint input. This should be called at least once, after
2369    the constraints have been defined using :func:`StoreEquivalence`,
2370    :func:`GroupConstraints` and :func:`GenerateConstraints` and before any
2371    parameter refinement is done.
2372
2373    This completes the parameter dictionary by defining values for parameters
2374    created by constraints based on the constraints that define them
2375    using the values for the current parameters. It also removes all dependent
2376    variables from the varyList
2377
2378    :param dict parmDict: a dict containing parameter values keyed by the
2379      parameter names. For new variables created by constraints, entries
2380      will be added to the dictionary, if not alreay present, or the
2381      values will be recomputed.
2382
2383    :param list varyList: a list of parameters names. Will be modified.
2384    '''
2385    # remove fixed parameters from the varyList
2386    for item in holdParmList+newHolds:
2387        if item in varyList: varyList.remove(item)
2388           
2389    # process the independent parameters:
2390    # * remove dependent ones from varylist
2391    # * for equivalences apply the independent parameters onto dependent variables
2392    global dependentParmList,arrayList,invarrayList,indParmList
2393    for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
2394        for item in varlist:
2395            if item in varyList: varyList.remove(item)
2396        if multarr is None:
2397            #for v,val in zip(  # shows values to be set
2398            #    varlist,
2399            #    np.dot(invmultarr,np.array([parmDict[var] for var in mapvars]))
2400            #    ): print('parmDict set',v,':',val)
2401            parmDict.update(zip(
2402                varlist,
2403                np.dot(invmultarr,np.array([parmDict[var] for var in mapvars]))
2404                ))
2405
2406    # * for the created parameters, compute them from their dependents
2407    for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
2408        if multarr is None: continue
2409        # evaluate all constraints in the forward direction
2410        for key,value in zip(mapvars,np.dot(multarr,np.array([parmDict[var] for var in varlist]))):
2411            if type(key) is float: continue
2412            #print('parmDict set',key,':',value)
2413            parmDict[key] = value
2414    global saveVaryList
2415    saveVaryList = copy.copy(varyList)
2416           
2417def Dict2Map(parmDict):
2418    '''Applies the constraints defined using :func:`StoreEquivalence`,
2419    :func:`GroupConstraints` and :func:`GenerateConstraints` by changing
2420    values in a dict containing the parameters. This should be
2421    done after refinement and before the parameters are used for
2422    any computations
2423
2424    :param dict parmDict: a dict containing parameter values keyed by the
2425      parameter names. After this is called, all the dependent variables
2426      will be updated based on constraints and equivalences.
2427    '''
2428    global dependentParmList,arrayList,invarrayList,indParmList
2429    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
2430        if invmultarr is None:  # is this needed?
2431            if GSASIIpath.GetConfigValue('debug'): 
2432                print('Why does this constraint have None for invmultarr?',varlist,mapvars)
2433            continue
2434        valslist = [parmDict.get(var,var) for var in mapvars]
2435        #for v,val in zip(varlist,np.dot(invmultarr,np.array(valslist))): print(v,val) # shows what is being set
2436        parmDict.update(zip(varlist,np.dot(invmultarr,np.array(valslist))))
2437       
2438#======================================================================
2439# internal routines follow (these routines are unlikely to be called
2440# from outside the module)
2441
2442def GramSchmidtOrtho(a,nkeep=0):
2443    '''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to
2444    find orthonormal unit vectors relative to first row.
2445
2446    If nkeep is non-zero, the first nkeep rows in the array are not changed
2447   
2448    input:
2449       arrayin: a 2-D non-singular square array
2450    returns:
2451       a orthonormal set of unit vectors as a square array
2452    '''
2453    def proj(a,b):
2454        'Projection operator'
2455        return a*(np.dot(a,b)/np.dot(a,a))
2456    for j in range(nkeep,len(a)):
2457        for i in range(j):
2458            a[j] -= proj(a[i],a[j])
2459        if np.allclose(np.linalg.norm(a[j]),0.0):
2460            raise ConstraintException("Singular input to GramSchmidtOrtho")
2461        a[j] /= np.linalg.norm(a[j])
2462    return a
2463
2464def _FillArray(sel,dict,collist,FillDiagonals=False):
2465    '''Construct a n by n matrix (n = len(collist)
2466    filling in the rows using the relationships defined
2467    in the dictionaries selected by sel
2468
2469    If FillDiagonals is True, diagonal elements in the
2470    array are set to 1.0
2471    '''
2472    n = len(collist)
2473    if FillDiagonals:
2474        arr = np.eye(n)
2475    else:
2476        arr = np.zeros(2*[n,])
2477    # fill the top rows
2478    for i,cnum in enumerate(sel):
2479        for j,var in enumerate(collist):
2480            arr[i,j] = dict[cnum].get(var,0)
2481    return arr
2482
2483def _SwapColumns(i,m,v):
2484    '''Swap columns in matrix m as well as the labels in v
2485    so that element (i,i) is replaced by the first non-zero element in row i after that element
2486
2487    Throws an exception if there are no non-zero elements in that row
2488    '''
2489    for j in range(i+1,len(v)):
2490        if not np.allclose(m[i,j],0):
2491            m[:,(i,j)] = m[:,(j,i)]
2492            v[i],v[j] = v[j],v[i]
2493            return
2494    else:
2495        raise ConstraintException('Singular input')
2496
2497def _RowEchelon(m,arr,collist):
2498    '''Convert the first m rows in Matrix arr to row-echelon form
2499    exchanging columns in the matrix and collist as needed.
2500
2501    throws an exception if the matrix is singular because
2502    the first m rows are not linearly independent
2503    '''
2504    for i in range(m):
2505        if np.allclose(arr[i,i],0):
2506            _SwapColumns(i,arr,collist)
2507        arr[i,:] /= arr[i,i] # normalize row
2508        # subtract current row from subsequent rows to set values to left of diagonal to 0
2509        for j in range(i+1,m):
2510            arr[j,:] -= arr[i,:] * arr[j,i]
2511
2512if __name__ == "__main__":
2513    parmdict = {}
2514    constrDict = [
2515        {'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},
2516        {'0:0:eA': 0.0},
2517        {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
2518        {'1::C(10,0,1)': 1.0, '2::C(10,0,1)': 1.0},
2519        {'1::AUiso:0': 1.0, '0::AUiso:0': 1.0},
2520        {'0::A0': 0.0}
2521        ]
2522    fixedList = ['5.0', '0', None, None, '1.0', '0']
2523    StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), ))
2524    #StoreEquivalence('1::atomx:3',('2::atomx:3', ('2::atomz:3',2,), )) # error: dependent & independent vars mixed
2525    #StoreEquivalence('1::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) # error: dependent vars repeated
2526    #StoreEquivalence('0:1:eA',('0:0:eA',)) # error: equiv & fixed
2527    #StoreEquivalence('0:99:Scale',('0:12:Scale',)) # error: equiv & constrained
2528    #StoreEquivalence('0:12:Scale',('0:99:Scale',)) # error: equiv & constrained
2529    varyList = ['2::atomx:3',
2530                '2::C(10,6,1)', '1::C(10,6,1)',
2531                '2::atomy:3', '2::atomz:3',
2532                '0:12:Scale', '0:11:Scale', '0:14:Scale', '0:13:Scale', '0:0:Scale']
2533#    varyList = ['0::A0', '0::AUiso:0', '0::Afrac:1', '0::Afrac:2', '0::Afrac:3', '0::Afrac:4',
2534#       '0::dAx:5', '0::dAy:5', '0::dAz:5', '0::AUiso:5', ':0:Back;0', ':0:Back;1', ':0:Back;2', ':0:Back;3',
2535#       ':0:Back;4', ':0:Back;5', ':0:Back;6', ':0:Back;7', ':0:Back;8', ':0:Back;9', ':0:Back;10', ':0:Back;11'
2536#       :0:U', ':0:V', ':0:W', ':0:X', ':0:Y', ':0:Scale', ':0:DisplaceX', ':0:DisplaceY']
2537#    constrDict = [
2538#        {'0::Afrac:4': 24.0, '0::Afrac:1': 16.0, '0::Afrac:3': 24.0, '0::Afrac:2': 16.0},
2539#        {'0::Afrac:1': 1.0, '0::Afrac:2': 1.0},
2540#        {'0::Afrac:4': 1.0, '0::Afrac:3': 1.0}]
2541#    fixedList = ['40.0', '1.0', '1.0']
2542
2543    msg,warning,groups,parmlist = GenerateConstraints(varyList,constrDict,fixedList,parmdict)
2544    print (VarRemapShow(varyList))
2545    parmdict.update( {
2546        '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,
2547        '0:0:eA': 0.0,
2548        '2::C(10,6,1)': 0.2, '1::C(10,6,1)': 0.3,
2549        '1::C(10,0,1)': 0.2, '2::C(10,0,1)': 0.3,
2550        '1::AUiso:0': 0.02, '0::AUiso:0': 0.03,
2551        '0::A0': 0.0,
2552        '2::atomx:3':0.23,'2::atomy:3':-.23, '2::atomz:3':-0.11,
2553        })
2554    print ('parmdict start',parmdict)
2555    print ('varylist start',varyList)
2556    before = parmdict.copy()
2557    Map2Dict(parmdict,varyList)
2558    print ('parmdict before and after Map2Dict')
2559    print ('  key / before / after')
2560    for key in sorted(list(parmdict.keys())):
2561        print ('  '+key,'\t',before.get(key),'\t',parmdict[key])
2562    print ('varylist after',varyList)
2563    before = parmdict.copy()
2564    Dict2Map(parmdict)
2565    print ('after Dict2Map')
2566    print ('  key / before / after')
2567    for key in sorted(list(parmdict.keys())):
2568        print ('  '+key,'\t',before.get(key),'\t',parmdict[key])
2569#    dMdv = len(varylist)*[0]
2570#    deriv = {}
2571#    for i,v in enumerate(parmdict.keys()): deriv[v]=i
2572#    Dict2Deriv(varylist,deriv,dMdv)
Note: See TracBrowser for help on using the repository browser.