source: trunk/GSASIImapvars.py @ 5145

Last change on this file since 5145 was 5145, checked in by toby, 9 months ago

fix testDeriv to handle derivatives of parameters created by constraints; fix their derivative computation

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