source: trunk/GSASIImapvars.py @ 5148

Last change on this file since 5148 was 5148, checked in by toby, 7 months ago

add/fix comments

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