source: trunk/GSASIImapvars.py @ 5149

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

more doc formatting cleanup

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