Context Navigation

source: trunk/GSASIImapvars.py @ 5051

Last change on this file since 5051 was 5040, checked in by toby, 4 years ago
fix sym-gen constraint display
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 111.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: