source: trunk/GSASIImapvars.py

Last change on this file was 5305, checked in by toby, 3 months ago

constraints fix

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