source: trunk/GSASIImapvars.py @ 4851

Last change on this file since 4851 was 4851, checked in by toby, 8 months ago

reformat equivalence pairs as Bob prefers

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 79.2 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2021-03-12 02:35:23 +0000 (Fri, 12 Mar 2021) $
4# $Author: toby $
5# $Revision: 4851 $
6# $URL: trunk/GSASIImapvars.py $
7# $Id: GSASIImapvars.py 4851 2021-03-12 02:35:23Z toby $
8########### SVN repository information ###################
9"""
10*GSASIImapvars: Parameter constraints*
11======================================
12
13Module to implements algebraic contraints, parameter redefinition
14and parameter simplification contraints.
15
16Types of constraints
17--------------------
18
19There are four ways to specify constraints, as listed below.
20Note that parameters are initially stored in the
21main section of the GSAS-II data tree under heading ``Constraints``.
22This dict has four keys, 'Hist', 'HAP', 'Global', and 'Phase',
23each containing a list of constraints. An additional set of constraints
24are generated for each phase based on symmetry considerations by calling
25:func:`GSASIIstrIO.GetPhaseData`.
26
27Note that in the constraints, as stored in the GSAS-II data tree, parameters
28are stored as :class:`GSASIIobj.G2VarObj` objects, as these objects allow for
29changes in numbering of phases, histograms and atoms.
30When they are interpreted (in :func:`GSASIIstrIO.ProcessConstraints`),
31references
32to numbered objects are resolved using the appropriate random ids and the
33parameter object is converted to a string of form ``ph:hst:VARNAM:at``.
34
35Alternate parameters (New Var)
36^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
37
38Parameter redefinition ("New Var" constraints)
39is done by creating an expression that relates several
40parameters:
41
42::
43
44   Mx1 * Px + My1 * Py +...
45   Mx2 * Px + Mz2 * Pz + ...
46
47where Pj is a GSAS-II parameter name and Mjk is a constant (float) multiplier.
48Alternately, multipliers Mjk can contain a formula (str) that will be evaluated prior
49to the start of the refinement. In a formula, GSAS-II parameters will be replaced by the
50value of the parameter before the formula is evaluated, so ``'np.cos(0::Ax:2)'`` is a valid
51multiplier. At present, only phase (atom/cell) parameters are available for use in
52a formula, but this can be expanded if needed.
53
54This type of constraint describes an alternate
55degree of freedom where parameter Px and Py, etc. are varied to keep
56their ratio
57fixed according the expression. A new variable parameter is assigned to each degree of
58freedom when refined. An example where this can be valuable is when
59two 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.
60In the later stages of refinement, a second
61variable, Pd = P1 - P2, can be defined and it can be seen if refining Pd is
62supported by the data. Another use will be to define parameters that
63express "irrep modes" in terms of the fundamental structural parameters.
64
65These "New Var" constraints are stored as described for type "f" in the
66:ref:`constraint definitions table <Constraint_definitions_table>`.
67
68Constrained parameters (Const)
69^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
70
71A constraint on a set of variables can be supplied in the form of a
72linear algebraic equation: ::
73
74  Nx1 * Px + Ny1 * Py +... = C1
75
76where Cn is a constant (float), where Pj is a GSAS-II parameter name,
77and where Njk is a constant multiplier (float) or a formula (str) that will be evaluated prior
78to the start of the refinement. In a formula, GSAS-II parameters will be replaced by the
79value of the parameter before the formula is evaluated, so ``'np.cos(0::Ax:2)'`` is a valid
80multiplier. At present, only phase (atom/cell) parameters are available for use in
81a formula, but this can be expanded if needed.
82
83These equations set an interdependence between parameters.
84Common uses of parameter constraints are to set rules that decrease the number of parameters,
85such as restricting the sum of fractional occupancies for atoms that share
86a site to sum to unity, thus reducing the effective number of variables by one.
87Likewise, the Uiso value for a H atom "riding" on a C, N or O atom
88can be related by a fixed offset (the so called B+1 "rule").
89
90A "Const" constraint is stored as described for type "c" in the
91:ref:`constraint definitions table <Constraint_definitions_table>`.
92
93Equivalenced parameters (Equiv)
94^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
95
96A simplifed way to set up a constraint equation is to define an equivalence,
97which can be of form: ::
98
99  C1 * P1 = C2 * Py
100
101or::
102
103  C1 * P1 = C2 * P2 = C3 * P3 = ...
104
105where Cn is a constant (float), where Pj is a GSAS-II parameter name. This
106means that parameters Py (or P2 and P3) are determined from (or "slaved" to)
107parameter P1. Alternately, equivalences can be created with :func:`StoreEquivalence`
108where the multipliers can be a formula (str) that will be evaluated prior to the start of
109the refinement. In a formula, GSAS-II parameters will be replaced by the value of the
110parameter before the formula is evaluate, so ``'np.cos(0::Ax:2)'`` is a valid multiplier.
111At present, only phase (atom/cell) parameters are available for use in
112a formula, but this can be expanded if needed.
113Note that
114the latter constraint expression is conceptually identical to
115defining constraint equations. In practice, however,
116equivalenced parameters are processed in a
117different and more direct manner than constraint equations. The previous
118set of equalities could also be written in this way as a set of constraint
119equations: ::
120
121  C1 * P1 - C2 * P2 = 0
122  C1 * P1 - C3 * P3 = 0
123  ...
124
125
126The first parameter (P1 above)
127is considered the independent variable
128and the remaining parameters are dependent variables. The dependent variables
129are set from the independent variable.
130An example of how this may be used woul be if, for example,
131a material has a number of O atoms, all in fairly similar bonding environments
132and the diffraction data are sparse, one my reduce the complexity of the model
133by defining Uiso for the first O atoms to be identical to the remaining atoms.
134The results of this refinement will be simpler to understand than if a set of
135constraint equations is used because the refined parameter will be the independent variable, which will be as ph::Uiso:n, corresponding to the first O atom.
136
137A parameter can be used in multiple
138equivalences as independent variable,
139but if parameter is used as both a dependent and independent variable or
140a parameter is used in equivalences and in "New Var" or "Const" constraints,
141this create conflicts that cannot be resolved within the equivalences implementation
142but can be handled as constraint equations.
143The equivalences that violate this are discovered in :func:`CheckEquivalences`
144and then :func:`MoveConfEquiv` is used to change these equivalences to "Const" equations.
145
146Equivalenced parameters ("EQUIV" constraints), when defined by users,
147are stored as described for type "e" in the
148:ref:`constraint definitions table <Constraint_definitions_table>`.
149Other equvalences are generated by symmetry prior
150to display or refinement in :func:`GSASIIstrIO.GetPhaseData`.
151These are not stored.
152
153Fixed parameters (Hold)
154^^^^^^^^^^^^^^^^^^^^^^^^
155
156When parameters are refined where a single refinement flag determines that several variables
157are refined at the same time (examples are: cell parameters, atom positions, anisotropic
158displacement parameters, magnetic moments,...) it can be useful to specify that a
159specific parameter should not be varied. These will most commonly be generated due to symmetry,
160but under specific conditions, there may be other good reasons to constrain a parameter.
161
162A "Hold" constraint is stored as described for type "h" in the
163:ref:`constraint definitions table <Constraint_definitions_table>`.
164
165.. _Constraints_processing:
166
167Constraint Processing
168---------------------
169
170When constraints will be used or edited, they are processed using a series of
171calls:
172
173* First all of the stored constraints are appended into a single list. They are
174  initially stored in separate lists only to improve their creation and display
175  in the GUI.
176
177* Then :func:`InitVars` is used to initialize the global variables in
178  this module (:mod:`GSASIImapvars`).
179
180* Then :func:`GSASIIstrIO.ProcessConstraints` is used to initially process the
181  constraints, as described below.
182
183* Symmetry-generated equivalences are then created in
184  :func:`GSASIIstrIO.GetPhaseData`, which also calls
185  :func:`GSASIIstrIO.cellVary` and for Pawley refinements
186  :func:`GSASIIstrIO.GetPawleyConstr`. These are entered directly into this
187  module's globals using :func:`StoreEquivalence`.
188* Constraints/equivalences are then checked for possible conflicts with
189  :func:`CheckConstraints`, this requires grouping the constraints,
190  as described below.
191* In refinements, :func:`GenerateConstraints` is then called to
192  create the constraints that will be used, see below for
193* For debugging constraints, :func:`VarRemapShow` can be called after
194  :func:`GenerateConstraints` to display the generated constraints.
195
196Constraint Reorganization (:func:`~GSASIIstrIO.ProcessConstraints`)
197^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
198
199:func:`GSASIIstrIO.ProcessConstraints` is used to initially process the
200constraints. This does these things:
201
2021. The "Hold", "Const" and "New Var" expressions are split between two paired lists,
203   :data:`constDictList` and :data:`fixedList` which are set:
204 
205   * For "Hold" entries a dict with a single entry is placed in constDictList where
206     the key is the parameter name (associated value is 0.0) and fixedList gets a
207     value of 0.0.
208   * For "Const" entries, a dict with multiple entries is placed in constDictList where
209     the key is the parameter name and the value is the multiplier for the parameter,
210     while fixedList gets a string value corresponding to the constant value for
211     the expression.
212   * For "New Var" entries, a dict with multiple entries is placed in constDictList
213     where the key is the parameter name and the value is the multiplier
214     for the parameter; an additional key "_vary" is given the value of True or False
215     depending on the refinement flag setting. The corresponding entry in
216     fixedList is None
217
218   The output from this will have this form where the first entry is a "Const", the
219   second is a "New Var" and the final is a "Hold".
220
221  .. code-block:: python
222
223    constrDict = [
224         {'0:12:Scale': 2.0, '0:14:Scale': 4.0, '0:13:Scale': 3.0, '0:0:Scale': 0.5},
225         {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0, '_vary':True},
226         {'0::A0': 0.0}]
227    fixedList = ['5.0', None, '0']
228
229
2302. Equivalences are stored using :func:`StoreEquivalence`
231   into this module's globals (:data:`~GSASIImapvars.arrayList`,
232   :data:`~GSASIImapvars.invarrayList`,    :data:`~GSASIImapvars.indParmList`,
233   :data:`~GSASIImapvars.dependentParmList` and  :data:`~GSASIImapvars.symGenList`)
234
235
236Parameter Grouping (:func:`GenerateConstraints`)
237^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
238
239Functions :func:`CheckConstraints` and :func:`GenerateConstraints` are used to
240process the parameter equivalences and constraint lists created in
241:func:`~GSASIIstrIO.ProcessConstraints`. The former is used to generate
242error messages and the latter to generate the internal information used to
243apply the constraints.
244
245Initially, in both a list of parameters that are fixed and those used in
246constraint relations are tabulated in :func:`CheckEquivalences`. The equivalence
247relations are the scanned for the following potential problems:
248
2491. a parameter is used as a dependent variable in more than one
250   equivalence relation
2512. a parameter is fixed and used in in an equivalence relation either
252   as a dependent or independent variable
2533. a parameter is used as a dependent variable in one equivalence relation and
254   as a independent variable in another
2554. a parameter is used in in an equivalence relation (either
256   as a dependent or independent variable) and is used in a constraint
257   expression
2585. a parameter is not defined in a particular refinement, but is used in an
259   equivalence relation
2606. a parameter uses a wildcard for the histogram number (sequential refinements)
261
262Cases 1 & 2 above cannot be corrected, and result in errors. Cases 3 & 4 are
263potentially corrected with :func:`MoveConfEquiv`, as described below.
264Case 5 causes the equivalence to
265be dropped. Case 6 causes the current histogram number to be substituted for
266the wildcard.
267
268For cases 3 & 4, :func:`MoveConfEquiv` is used to change these equivalences
269into "Const" equations. This can potentially mean that additional
270equivalences will be problematic, so if there are changes made by
271:func:`MoveConfEquiv`, :func:`CheckEquivalences` is repeated. If any problem
272cases are noted, the refinement cannot be performed.
273
274Constraint expressions ("Const" and "New Var") are sorted into
275groups so that each group contains the minimum number of entries that
276ensures each parameter is referenced in only one group in
277:func:`GroupConstraints`. This is done by scanning the
278list of dicts in :data:`constDictList` one by one and making a list
279of parameters used in that constraint expression. Any expression that contains
280a parameter in is in that list is added to the current group and those
281parameters are added to this list of parameters. The list of ungrouped
282expressions is then scanned again until no more expressions are added to the
283current group. This process is repeated until every expression has been
284placed in a group. Function :func:`GroupConstraints` returns two lists of lists.
285The first has, for each group, a list of the indices in :data:`constDictList`
286that comprise the group (there can be only one). The second list contains,
287for each group, the unique parameter names in that group.
288
289Each constraint group is then processed. First, wildcard parameters are
290renamed (in a sequential refinement). Any fixed parameters that are used
291in constraints are noted as errors. The number of refined parameters and
292the number of parameters that are not defined in the current refinement are
293also noted. It is fine if all parameters in a group are not defined or all are
294not varied, but if some are defined and others not or some are varied and
295others not, this constitutes an error.
296
297The contents of each group is then examined. Groups with a single
298parameter (holds) are ignored. Then for each group, the number
299of parameters in the group (Np) and the number of expressions in the
300group (Nc) are counted and for each expression. If Nc > Np, then the constraint
301is overdetermined, which also constitutes an error.
302
303The parameter multipliers for each expression are then assembled:
304
305::
306
307   M1a * P1 + M2a * P2 +... Mka * Pk
308   M1b * P1 + M2b * P2 +... Mkb * Pk
309   ...
310   M1j * P1 + M2j * P2 +... Mkj * Pk
311
312
313From this it becomes possible to create a Nc x Np matrix, which is
314called the constraint matrix:
315
316 .. math::
317
318   \\left( \\begin{matrix}
319   M_{1a}  & M_{2a} &... & M_{ka} \\\\
320   M_{1b}  & M_{2b} &... & M_{kb} \\\\
321   ...  \\\\
322   M_{1j}  & M_{2j}  &... & M_{kj}
323   \\end{matrix}\\right)
324
325When Nc<Np, then additional rows need to be added to the matrix and to
326the vector that contains the value for each row (:data:`fixedList`) where
327values are ``None`` for New Vars and a constant for fixed values.
328This then can describe a system of Np simultaneous equations:
329
330 .. math::
331
332   \\left( \\begin{matrix}
333   M_{1a}  & M_{2a} &... & M_{ka} \\\\
334   M_{1b}  & M_{2b} &... & M_{kb} \\\\
335   ...  \\\\
336   M_{1j}  & M_{2j}  &... & M_{kj}
337   \\end{matrix}\\right)
338   \\left( \\begin{matrix}
339   P_{1} \\\\
340   P_{2} \\\\
341   ...  \\\\
342   P_{k}
343   \\end{matrix}\\right)
344   =
345   \\left( \\begin{matrix}
346   C_{1} & C_{2} &  ... & C_{k}
347   \\end{matrix}\\right)
348
349This is done by creating a square matrix from the group using
350:func:`_FillArray` with parameter ``FillDiagonals=False`` (the default). Any
351unspecified rows are left as all zero. The first Nc rows in the
352array are then coverted to row-echelon form using :func:`_RowEchelon`. This
353will create an Exception if any two rows are linearly dependent (which means
354that no matter what values are used for the remaining rows, that the matrix
355will be singular). :func:`_FillArray` is then called with parameter
356``FillDiagonals=True``, which again creates a square matrix but where
357unspecified rows are zero except for the diagonal elements. The 
358`Gram-Schmidt process <http://en.wikipedia.org/wiki/Gram-Schmidt>`_,
359implemented  in :func:`GramSchmidtOrtho`, is used to find orthonormal unit
360vectors for the remaining Np-Nc rows of the matrix. This will fail with
361a ``ConstraintException`` if this is not possible (singular matrix) or
362the result is placed in :data:`constrArr` as a numpy array.
363
364Rows in the matrix corresponding to "New Var" constraints and those that
365were generated by the Gram-Schmidt process are provided with parameter names
366(this can be specified if a "New Var" entry by using a ``"_name"`` element
367in the constraint dict, but at present this is not implemented.) Names are
368generated using :data:`paramPrefix` which is set to ``"::constr"``, plus a
369number to make the new parameter name unique. Global dict :data:`genVarLookup`
370provides a lookup table, where the names of the parameters related to this new
371parameter can be looked up easily.
372
373Finally the parameters used as input to the constraint are placed in
374this module's globals
375:data:`~GSASIImapvars.dependentParmList` and the constraint matrix is
376placed in into  :data:`~GSASIImapvars.arrayList`. This can be used to compute
377the initial values for "New Var" parameters. The inverse of the
378constraint matrix is placed in :data:`~GSASIImapvars.invarrayList` and a list
379of the "New Var" parameters and a list of the fixed values (as str's)
380is placed in :data:`~GSASIImapvars.indParmList`. A lookup table for
381fixed values as floats is placed in :data:`~GSASIImapvars.fixedDict`.
382Finally the appropriate entry in :data:`~GSASIImapvars.symGenList` is set to
383False to indicate that this is not a symmetry generated constraint.
384
385
386*Externally-Accessible Routines*
387---------------------------------
388
389To define a set of constrained and unconstrained relations, one
390defines a list of dictionary defining constraint parameters and their
391values, a list of fixed values for each constraint and a list of
392parameters to be varied. In addition, one uses
393:func:`StoreEquivalence` to define parameters that are equivalent.
394Use :func:`EvaluateMultipliers` to convert formula-based constraint/equivalence
395multipliers to numbers and then
396use :func:`CheckConstraints` to check that the input is
397internally consistent and finally :func:`GroupConstraints` and
398:func:`GenerateConstraints` to generate the internally used
399tables. Routine :func:`Map2Dict` is used to initialize the parameter
400dictionary and routine :func:`Dict2Map`, :func:`Dict2Deriv`, and
401:func:`ComputeDepESD` are used to apply constraints. Routine
402:func:`VarRemapShow` is used to print out the constraint information,
403as stored by :func:`GenerateConstraints`. Further information on each routine
404is below:
405
406:func:`InitVars`
407  This is optionally used to clear out all defined previously defined constraint information
408 
409:func:`StoreEquivalence`
410  To implement parameter redefinition, one calls StoreEquivalence. This should be called for every set of
411  equivalence relationships. There is no harm in using StoreEquivalence with the same independent variable:
412
413  .. code-block:: python
414
415       StoreEquivalence('x',('y',))
416       StoreEquivalence('x',('z',))
417
418  or equivalently
419
420  .. code-block:: python
421
422       StoreEquivalence('x',('y','z'))
423
424  The latter will run more efficiently. Note that mixing independent
425  and dependent variables would require assignments, such as
426
427  .. code-block:: python
428
429        StoreEquivalence('x',('y',))
430        StoreEquivalence('y',('z',))
431
432  would require that equivalences be applied in a particular order and
433  thus is implemented as a constraint equation rather than an equivalence.
434       
435  Use StoreEquivalence before calling GenerateConstraints or CheckConstraints
436
437:func:`CheckConstraints`
438   check that input in internally consistent
439
440:func:`GenerateConstraints`
441   generate the internally used tables from constraints and equivalences
442
443:func:`EvaluateMultipliers`
444   Convert any string-specified (formula-based) multipliers to numbers. Call this before
445   using :func:`CheckConstraints` or :func:`GenerateConstraints`.
446   At present, the code may pass only the dict for phase (atom/cell) parameters,
447   but this could be expanded if needed.
448
449:func:`Map2Dict`
450   To determine values for the parameters created in this module, one
451   calls Map2Dict. This will not apply contraints.
452
453:func:`Dict2Map`
454   To take values from the new independent parameters and constraints,
455   one calls Dict2Map and set the parameter array, thus appling contraints.
456
457:func:`Dict2Deriv`
458   Use Dict2Deriv to determine derivatives on independent parameters
459   from those on dependent ones.
460
461:func:`ComputeDepESD`     
462   Use ComputeDepESD to compute uncertainties on dependent variables.
463
464:func:`VarRemapShow`
465   To show a summary of the parameter remapping, one calls VarRemapShow.
466
467*Global Variables*
468------------------
469
470dependentParmList:
471   a list containing group of lists of
472   parameters used in the group. Note that parameters listed in
473   dependentParmList should not be refined as they will not affect
474   the model
475
476indParmList:
477     a list containing groups of Independent parameters defined in
478     each group. This contains both parameters used in parameter
479     redefinitions as well as names of generated new parameters.
480
481arrayList:
482     a list containing group of relationship matrices to relate
483     parameters in dependentParmList to those in indParmList. Unlikely
484     to be used externally.
485
486invarrayList:
487     a list containing group of relationship matrices to relate
488     parameters in indParmList to those in dependentParmList. Unlikely
489     to be used externally.
490
491fixedVarList:
492     a list of parameters that have been 'fixed'
493     by defining them as equal to a constant (::var: = 0). Note that
494     the constant value is ignored at present. These parameters are
495     later removed from varyList which prevents them from being refined.
496     Unlikely to be used externally.
497
498fixedDict:
499     a dictionary containing the fixed values corresponding
500     to parameter equations.  The dict key is an ascii string, but the
501     dict value is a float.  Unlikely to be used externally.
502
503symGenList:
504     a list of boolean values that will be True to indicate that a constraint
505     (only equivalences) is generated by symmetry (or Pawley overlap)
506     
507problemVars:
508     a list containing parameters that show up in constraints producing errors
509
510
511
512*Routines/variables*
513---------------------
514
515Note that parameter names in GSAS-II are strings of form ``<ph#>:<hst#>:<nam>`` or ``<ph#>::<nam>:<at#>``.
516"""
517
518from __future__ import division, print_function
519import numpy as np
520import sys
521import GSASIIpath
522GSASIIpath.SetVersionNumber("$Revision: 4851 $")
523# data used for constraints;
524debug = False # turns on printing as constraint input is processed
525
526# note that constraints are stored listed by contraint groups,
527# where each constraint
528# group contains those parameters that must be handled together
529dependentParmList = []
530'''a list of lists where each item contains a list of parameters in each constraint group.
531note that parameters listed in dependentParmList should not be refined directly.'''
532indParmList = [] # a list of names for the new parameters
533'''a list of lists where each item contains a list for each constraint group with
534fixed values for constraint equations and names of generated (New Var) parameters.
535'''
536arrayList = []
537'''a list of of relationship matrices that map model parameters in each
538constraint group (in :data:`dependentParmList`) to
539generated (New Var) parameters.
540'''
541invarrayList = []
542'''a list of of inverse-relationship matrices that map constrained values and
543generated (New Var) parameters (in :data:`indParmList`) to model parameters
544(in :data:`dependentParmList`).
545'''
546fixedDict = {}
547'''A dict lookup-table containing the fixed values corresponding
548to defined parameter equations. Note the key is the original ascii string
549and the value in the dict is a float.
550'''
551fixedVarList = []
552'''List of parameters that should not be refined.
553'''
554symGenList = []
555'''A list of flags that if True indicates a constraint was generated by symmetry
556'''
557problemVars = []
558'''a list of parameters causing errors
559'''
560dependentVars = []
561'A list of dependent variables, taken from (:data:`dependentParmList`).'
562independentVars = []
563'A list of dependent variables, taken from (:data:`indParmList`).'
564genVarLookup = {}
565'provides a list of parameters that are related to each generated parameter'
566paramPrefix = "::constr"
567'A prefix for generated parameter names'
568consNum = 0
569'The number to be assigned to the next constraint to be created'
570
571class ConstraintException(Exception):
572    '''Defines an Exception that is used when an exception is raised processing constraints
573    '''
574    pass
575
576def InitVars():
577    '''Initializes all constraint information'''
578    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,consNum,symGenList
579    dependentParmList = [] # contains a list of parameters in each group
580    arrayList = [] # a list of of relationship matrices
581    invarrayList = [] # a list of inverse relationship matrices
582    indParmList = [] # a list of names for the new parameters
583    fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations
584    symGenList = [] # Flag if constraint is generated by symmetry
585    consNum = 0 # number of the next constraint to be created
586    global genVarLookup
587    genVarLookup = {}
588
589def VarKeys(constr):
590    """Finds the keys in a constraint that represent parameters
591    e.g. eliminates any that start with '_'
592
593    :param dict constr: a single constraint entry of form::
594
595        {'var1': mult1, 'var2': mult2,... '_notVar': val,...}
596
597        (see :func:`GroupConstraints`)
598    :returns: a list of keys where any keys beginning with '_' are
599      removed.
600    """
601    return [i for i in constr.keys() if not i.startswith('_')]
602
603
604def GroupConstraints(constrDict):
605    """divide the constraints into groups that share no parameters.
606
607    :param dict constrDict: a list of dicts defining relationships/constraints
608
609    ::
610   
611       constrDict = [{<constr1>}, {<constr2>}, ...]
612
613    where {<constr1>} is {'var1': mult1, 'var2': mult2,... }
614
615    :returns: two lists of lists:
616   
617      * a list of grouped contraints where each constraint grouped containts a list
618        of indices for constraint constrDict entries
619      * a list containing lists of parameter names contained in each group
620     
621      """
622    assignedlist = [] # relationships that have been used
623    groups = [] # contains a list of grouplists
624    ParmList = []
625    for i,consi in enumerate(constrDict):
626        if i in assignedlist: continue # already in a group, skip
627        # starting a new group
628        grouplist = [i,]
629        assignedlist.append(i)
630        groupset = set(VarKeys(consi))
631        changes = True # always loop at least once
632        while(changes): # loop until we can't find anything to add to the current group
633            changes = False # but don't loop again unless we find something
634            for j,consj in enumerate(constrDict):
635                if j in assignedlist: continue # already in a group, skip
636                if len(set(VarKeys(consj)) & groupset) > 0: # true if this needs to be added
637                    changes = True
638                    grouplist.append(j)
639                    assignedlist.append(j)
640                    groupset = groupset | set(VarKeys(consj))
641        group = sorted(grouplist)
642        varlist = sorted(list(groupset))
643        groups.append(group)
644        ParmList.append(varlist)
645    return groups,ParmList
646
647def CheckConstraints(varyList,constrDict,fixedList):
648    '''Takes a list of relationship entries comprising a group of
649    constraints and checks for inconsistencies such as conflicts in
650    parameter/variable definitions and or inconsistently varied parameters.
651
652    :param list varyList: a list of parameters names that will be varied
653
654    :param dict constrDict: a list of dicts defining relationships/constraints
655      (as created in :func:`GSASIIstrIO.ProcessConstraints` and
656      documented in :func:`GroupConstraints`)
657
658    :param list fixedList: a list of values specifying a fixed value for each
659      dict in constrDict. Values are either strings that can be converted to
660      floats or ``None`` if the constraint defines a new parameter rather
661      than a constant.
662
663    :returns: two strings:
664
665      * the first lists conflicts internal to the specified constraints
666      * the second lists conflicts where the varyList specifies some
667        parameters in a constraint, but not all
668       
669      If there are no errors, both strings will be empty
670    '''
671    import re
672    global dependentParmList,arrayList,invarrayList,indParmList,consNum
673    global problemVars
674    # Process the equivalences
675    #    If there are conflicting parameters, move them into constraints. This
676    #    may create new conflicts, requiring movement of other parameters, so
677    #    loop until there are no more changes to make.
678    parmsChanged = True
679    while parmsChanged:
680        parmsChanged = 0
681        errmsg,warnmsg,fixVlist,dropVarList,translateTable = CheckEquivalences(
682            constrDict,varyList)
683        #print('debug: using MoveConfEquiv to address =',errmsg)
684        if problemVars: parmsChanged = MoveConfEquiv(constrDict,fixedList)
685#    GSASIIpath.IPyBreak()
686
687    groups,parmlist = GroupConstraints(constrDict)
688    # scan through parameters in each relationship. Are all varied? If only some are
689    # varied, create a warning message.
690    for group,varlist in zip(groups,parmlist):
691        if len(varlist) == 1:   # process fixed (held) variables
692            var = varlist[0]
693            if var not in fixedVarList:
694                fixedVarList.append(var)
695            continue
696        for rel in group:
697            varied = 0
698            notvaried = ''
699            for var in constrDict[rel]:
700                if var.startswith('_'): continue
701                if not re.match('[0-9]*:[0-9\*]*:',var):
702                    warnmsg += "\nParameter "+str(var)+" does not begin with a ':'"
703                if var in varyList:
704                    varied += 1
705                else:
706                    if notvaried: notvaried += ', '
707                    notvaried += var
708                if var in fixVlist:
709                    errmsg += '\nParameter '+var+" is Fixed and used in a constraint:\n\t"
710                    if var not in problemVars: problemVars.append(var)
711                    errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
712            if varied > 0 and varied != len(VarKeys(constrDict[rel])):
713                warnmsg += "\nNot all parameters refined in constraint:\n\t"
714                warnmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
715                warnmsg += '\nNot refined: ' + notvaried + '\n'
716    if errmsg or warnmsg:
717        return errmsg,warnmsg
718
719    # now look for process each group and create the relations that are needed to form
720    # non-singular square matrix
721    for group,varlist in zip(groups,parmlist):
722        if len(varlist) == 1: continue # a constraint group with a single parameter can be ignored
723        if len(varlist) < len(group): # too many relationships -- no can do
724            errmsg += "\nOver-constrained input. "
725            errmsg += "There are more constraints " + str(len(group))
726            errmsg += "\n\tthan parameters " + str(len(varlist)) + "\n"
727            for rel in group:
728                errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
729                errmsg += "\n"
730                continue
731        try:
732            multarr = _FillArray(group,constrDict,varlist)
733            _RowEchelon(len(group),multarr,varlist)
734        except:
735            errmsg += "\nSingular input. "
736            errmsg += "There are internal inconsistencies in these constraints\n"
737            for rel in group:
738                errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
739                errmsg += "\n"
740            continue
741        try:
742            multarr = _FillArray(group,constrDict,varlist,FillDiagonals=True)
743            GramSchmidtOrtho(multarr,len(group))
744        except:
745            errmsg += "\nUnexpected singularity with constraints (in Gram-Schmidt)\n"
746            for rel in group:
747                errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
748                errmsg += "\n"
749            continue
750        mapvar = []
751        group = group[:]
752        # scan through all generated and input parameters
753        # Check again for inconsistent parameter use
754        # for new parameters -- where varied and unvaried parameters get grouped
755        # together. I don't think this can happen when not flagged before, but
756        # it does not hurt to check again.
757        for i in range(len(varlist)):
758            varied = 0
759            notvaried = ''
760            if len(group) > 0:
761                rel = group.pop(0)
762                fixedval = fixedList[rel]
763                for var in VarKeys(constrDict[rel]):
764                    if var in varyList:
765                        varied += 1
766                    else:
767                        if notvaried: notvaried += ', '
768                        notvaried += var
769            else:
770                fixedval = None
771            if fixedval is None:
772                varname = paramPrefix + str(consNum) # assign a name to a parameter
773                mapvar.append(varname)
774                consNum += 1
775            else:
776                mapvar.append(fixedval)
777            if varied > 0 and notvaried != '':
778                warnmsg += "\nNot all parameters refined in generated constraint"
779                warnmsg += '\nPlease report this unexpected error\n'
780                for rel in group:
781                    warnmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
782                    warnmsg += "\n"
783                warnmsg += '\n\tNot refined: ' + notvaried + '\n'
784        try:
785            np.linalg.inv(multarr)           
786        except:
787            errmsg += "\nSingular input. "
788            errmsg += "The following constraints are not "
789            errmsg += "linearly independent\n\tor do not "
790            errmsg += "allow for generation of a non-singular set\n"
791            errmsg += 'Please report this unexpected error\n'
792            for rel in group:
793                errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
794                errmsg += "\n"
795    _setVarLists([])
796    return errmsg,warnmsg
797
798def GenerateConstraints(varyList,constrDict,fixedList,parmDict=None,SeqHist=None):
799    '''Takes a list of relationship entries comprising a group of
800    constraints and builds the relationship lists and their inverse
801    and stores them in global parameters Also checks for internal
802    conflicts or inconsistencies in parameter/variable definitions.
803
804    :param list varyList: a list of parameters names (strings of form
805      ``<ph>:<hst>:<nam>``) that will be varied. Note that this is changed here.
806   
807    :param dict constrDict: a list of dicts defining relationships/constraints
808      (as defined in :func:`GroupConstraints`)
809
810    :param list fixedList: a list of values specifying a fixed value for each
811      dict in constrDict. Values are either strings that can be converted to
812      floats, float values or None if the constraint defines a new parameter.
813     
814    :param dict parmDict: a dict containing all parameters defined in current
815      refinement.
816
817    :param int SeqHist: number of current histogram, when used in a sequential
818      refinement. None (default) otherwise. Wildcard parameter names are
819      set to the current histogram, when found if not None.
820    '''
821    global dependentParmList,arrayList,invarrayList,indParmList,consNum
822    global genVarLookup
823    msg = ''
824    shortmsg = ''
825    changed = False
826   
827    # Process the equivalences
828    #    If there are conflicting parameters, move them into constraints. This
829    #    may create new conflicts, requiring movement of other parameters, so
830    #    loop until there are no more changes to make.
831    parmsChanged = True
832    while parmsChanged:
833        parmsChanged = 0
834        errmsg,warnmsg,fixVlist,dropVarList,translateTable = CheckEquivalences(
835            constrDict,varyList,parmDict,SeqHist)
836        if problemVars:
837            parmsChanged = MoveConfEquiv(constrDict,fixedList)
838            changed = True
839    if errmsg:
840        msg = errmsg
841    if warnmsg:
842        if msg: msg += '\n'
843        msg += warnmsg
844
845    # scan through parameters in each relationship. Are all varied? If only some are
846    # varied, create an error message.
847    groups,parmlist = GroupConstraints(constrDict)
848    for group,varlist in zip(groups,parmlist):
849        if len(varlist) == 1:   # process fixed (held) variables
850            var = varlist[0]
851            if var not in fixedVarList:
852                fixedVarList.append(var)
853            continue
854        for rel in group:
855            varied = 0
856            notvaried = ''
857            unused = 0
858            notused = ''
859            for var in constrDict[rel]:
860                if var.startswith('_'): continue
861                if var.split(':')[1] == '*' and SeqHist is not None:
862                    # convert wildcard var to reference current histogram; save translation in table
863                    sv = var.split(':')
864                    sv[1] = str(SeqHist)
865                    translateTable[var] = ':'.join(sv)
866                    var = translateTable[var]
867                if parmDict is not None and var not in parmDict:
868                    unused += 1
869                    if notvaried: notused += ', '
870                    notused += var
871                if var in varyList:
872                    varied += 1
873                else:
874                    if notvaried: notvaried += ', '
875                    notvaried += var
876                if var in fixedVarList:
877                    msg += '\nError: parameter '+var+" is Fixed and used in a constraint:\n\t"
878                    msg += _FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
879            #if unused > 0:# and unused != len(VarKeys(constrDict[rel])):
880            if unused > 0 and unused != len(VarKeys(constrDict[rel])):
881                #msg += "\nSome (but not all) parameters in constraint are not defined:\n\t"
882                #msg += _FormatConstraint(constrDict[rel],fixedList[rel])
883                #msg += '\nNot used: ' + notused + '\n'
884                shortmsg += notused+" not used in constraint "+_FormatConstraint(constrDict[rel],fixedList[rel])
885            elif varied > 0 and varied != len(VarKeys(constrDict[rel])):
886                #msg += "\nNot all parameters refined in constraint:\n\t"
887                #msg += _FormatConstraint(constrDict[rel],fixedList[rel])
888                #msg += '\nNot refined: ' + notvaried + '\n'
889                shortmsg += notvaried+" not varied in constraint "+_FormatConstraint(constrDict[rel],fixedList[rel])
890    # if there were errors found, go no farther
891    if shortmsg and SeqHist is not None:
892        if msg:
893            print (' *** ERROR in constraint definitions! ***')
894            print (msg)
895            raise ConstraintException
896        print ('*** Sequential refinement: ignoring constraint definition(s): ***')
897        print (shortmsg)
898        msg = ''
899    elif shortmsg:
900        msg += shortmsg
901    if msg:
902        print (' *** ERROR in constraint definitions! ***')
903        print (msg)
904        raise ConstraintException
905               
906    # now process each group and create the relations that are needed to form
907    # a non-singular square matrix
908    # If all are varied and this is a constraint equation, then set VaryFree flag
909    # so that the newly created relationships will be varied
910    for group,varlist in zip(groups,parmlist):
911        if len(varlist) == 1: continue
912        # for constraints, if all included parameters are refined,
913        # set the VaryFree flag, and remaining degrees of freedom will be
914        # varied (since consistency was checked, if any one parameter is
915        # refined, then assume that all are)
916        varsList = [] # make a list of all the referenced parameters as well
917        VaryFree = False
918        for rel in group:
919            varied = 0
920            unused = 0
921            for var in VarKeys(constrDict[rel]):
922                var = translateTable.get(var,var) # replace wildcards
923                if parmDict is not None and var not in parmDict:
924                    unused += 1                   
925                if var not in varsList: varsList.append(var)
926                if var in varyList: varied += 1
927            if fixedList[rel] is not None and varied > 0:
928                VaryFree = True
929        if len(varlist) < len(group): # too many relationships -- no can do
930            msg = 'too many relationships'
931            break
932        # Since we checked before, if any parameters are unused, then all must be.
933        # If so, this set of relationships can be ignored
934        if unused:
935            if debug: print('Constraint ignored (all parameters undefined)')
936            if debug: print ('    '+_FormatConstraint(constrDict[rel],fixedList[rel]))
937            continue
938        # fill in additional degrees of freedom
939        try:
940            arr = _FillArray(group,constrDict,varlist)
941            _RowEchelon(len(group),arr,varlist)
942            constrArr = _FillArray(group,constrDict,varlist,FillDiagonals=True)
943            GramSchmidtOrtho(constrArr,len(group))
944        except:
945            msg = 'Singular relationships'
946            break
947        mapvar = []
948        group = group[:]
949        # scan through all generated and input relationships, we need to add to the varied list
950        # all the new parameters where VaryFree has been set or where a New Var is varied.
951        #
952        # If a group does not contain any fixed values (constraint equations)
953        # and nothing in the group is varied, drop this group, so that the
954        # dependent parameters can be refined individually.
955        unused = True
956        for i in range(len(varlist)):
957            if len(group) > 0: # get the original equation reference
958                rel = group.pop(0)
959                fixedval = fixedList[rel]
960                varyflag = constrDict[rel].get('_vary',False)
961                varname = constrDict[rel].get('_name','')
962            else: # this relationship has been generated
963                varyflag = False
964                varname = ''
965                fixedval = None
966            if fixedval is None: # this is a new parameter, not a constraint
967                if not varname:
968                    varname = paramPrefix + str(consNum) # no assigned name, create one
969                    consNum += 1
970                mapvar.append(varname)
971                genVarLookup[varname] = varlist # save list of parameters related to this new var
972                # vary the new relationship if it is a degree of freedom in
973                # a set of contraint equations or if a New Var is flagged to be varied.
974                if VaryFree or varyflag: 
975                    unused = False
976                    varyList.append(varname)
977                    # fix (prevent varying) of all the parameters inside the constraint group
978                    # (dependent vars)
979                    for var in varsList:
980                        if var in varyList: varyList.remove(var)
981            else:
982                unused = False
983                mapvar.append(fixedval)
984        if unused: # skip over constraints that don't matter (w/o fixed value or any refined parameters)
985            if debug: print('Constraint ignored (all parameters unrefined)')
986            if debug: print ('   '+_FormatConstraint(constrDict[rel],fixedList[rel]))
987            continue 
988        dependentParmList.append([translateTable.get(var,var) for var in varlist])
989        arrayList.append(constrArr)
990        invarrayList.append(np.linalg.inv(constrArr))
991        indParmList.append(mapvar)
992        symGenList.append(False)
993    if msg:
994        print (' *** ERROR in constraint definitions! ***')
995        print (msg)
996        print (VarRemapShow(varyList))
997        raise ConstraintException
998    # setup dictionary containing the fixed values
999    global fixedDict
1000    # key is original ascii string, value is float
1001    for fixedval in fixedList:
1002        if fixedval:
1003            fixedDict[fixedval] = float(fixedval)
1004    _setVarLists(dropVarList)
1005    if changed:
1006        print(60*'=')
1007        print('Constraints were reclassified to avoid conflicts, as below:')
1008        print (VarRemapShow(varyList,True))
1009        print(60*'=')
1010    return groups,parmlist # saved for sequential fits
1011   
1012def _setVarLists(dropVarList):
1013    '''Make list of dependent and independent variables (after dropping unused vars in dropVarList)
1014    '''
1015    global dependentParmList,indParmList
1016    global dependentVars
1017    global independentVars
1018    dependentVars = []
1019    independentVars = []
1020    for varlist,mapvars in zip(dependentParmList,indParmList):  # process all constraints
1021        for mv in mapvars:
1022            if mv in dropVarList: continue
1023            if mv not in independentVars: independentVars.append(mv)
1024        for mv in varlist:
1025            if mv in dropVarList: continue
1026            if mv not in dependentVars: dependentVars.append(mv)
1027    if debug: # on debug, show what is parsed & generated, semi-readable
1028        print (50*'-')
1029        #print (VarRemapShow(varyList))
1030        #print ('Varied: ',varyList)
1031        print ('Not Varied: ',fixedVarList)
1032
1033# def CheckEquivalences(constrDict,varyList):
1034#     global dependentParmList,arrayList,invarrayList,indParmList,consNum
1035#     global problemVars
1036#     warnmsg = ''
1037#     errmsg = ''
1038#     problemVars = []
1039#     # process fixed variables (holds)
1040#     fixVlist = [] # list of Fixed vars
1041#     constrVars = [] # list of vars in constraint expressions
1042#     for cdict in constrDict:
1043#         # N.B. No "_" names in holds
1044#         if len(cdict) == 1:
1045#             fixVlist.append(list(cdict.keys())[0])
1046#         else:
1047#             constrVars += cdict.keys() # this will include _vary (not a problem)
1048#     # process equivalences: make a list of dependent and independent vars
1049#     #    and check for repeated uses (repetition of a parameter as an
1050#     #    independent var is OK)
1051#     indepVarList = []
1052#     depVarList = []
1053#     multdepVarList = []
1054#     for varlist,mapvars,multarr,invmultarr in zip(
1055#         dependentParmList,indParmList,arrayList,invarrayList):
1056#         if multarr is None: # an equivalence
1057#             zeromult = False
1058#             for mv in mapvars:
1059#                 varied = 0
1060#                 notvaried = ''
1061#                 if mv in varyList:
1062#                     varied += 1
1063#                 else:
1064#                     if notvaried: notvaried += ', '
1065#                     notvaried += mv
1066#                 if mv not in indepVarList: indepVarList.append(mv)
1067#                 for v,m in zip(varlist,invmultarr):
1068#                     if v in indepVarList:
1069#                         errmsg += '\nVariable '+v+' is used to set values in a constraint before its value is set in another constraint\n'
1070#                         if v not in problemVars: problemVars.append(v)
1071#                     if m == 0: zeromult = True
1072#                     if v in varyList:
1073#                         varied += 1
1074#                     else:
1075#                         if notvaried: notvaried += ', '
1076#                         notvaried += v
1077#                     if v in depVarList:
1078#                         multdepVarList.append(v)
1079#                     else:
1080#                         depVarList.append(v)
1081#             if varied > 0 and varied != len(varlist)+1:
1082#                 warnmsg += "\nNot all variables refined in equivalence:\n\t"
1083#                 s = ""
1084#                 for v in varlist:
1085#                     if s != "": s+= " & "
1086#                     s += str(v)           
1087#                 warnmsg += str(mv) + " => " + s
1088#                 warnmsg += '\nNot refined: ' + notvaried + '\n'
1089#             if zeromult:
1090#                 errmsg += "\nZero multiplier is invalid in equivalence:\n\t"
1091#                 s = ""
1092#                 for v in varlist:
1093#                     if s != "": s+= " & "
1094#                     s += str(v)           
1095#                 errmsg += str(mv) + " => " + s + '\n'
1096#     # check for errors:
1097#     if len(multdepVarList) > 0:
1098#         errmsg += "\nThe following parameters(s) are used in conflicting Equivalence relations as dependent variables:\n"
1099#         s = ''
1100#         for var in sorted(set(multdepVarList)):
1101#             if v not in problemVars: problemVars.append(v)
1102#             if s != "": s+= ", "
1103#             s += str(var)           
1104#         errmsg += '\t'+ s + '\n'
1105#     equivVarList = list(set(indepVarList).union(set(depVarList)))
1106#     if debug: print ('equivVarList',equivVarList)
1107#     # check for parameters that are both fixed and in an equivalence (not likely)
1108#     inboth = set(fixVlist).intersection(set(equivVarList))
1109#     if len(inboth) > 0:
1110#         errmsg += "\nThe following parameter(s) are used in both Equivalence and Fixed constraints:\n"
1111#         s = ''
1112#         for var in sorted(inboth):
1113#             if var not in problemVars: problemVars.append(var)
1114#             if s != "": s+= ", "
1115#             s += str(var)
1116#         errmsg += '\t'+ s + '\n'
1117#     # check for parameters that in both an equivalence and a constraint expression
1118#     inboth = set(constrVars).intersection(set(equivVarList))
1119#     if len(inboth) > 0:
1120#         errmsg += "\nThe following parameter(s) are used in both Equivalence and Equiv or new var constraints:\n"
1121#         s = ''
1122#         for var in sorted(inboth):
1123#             if var not in problemVars: problemVars.append(var)
1124#             if s != "": s+= ", "
1125#             s += str(var)
1126#         errmsg += '\t'+ s + '\n'
1127#     return errmsg,warnmsg,fixVlist
1128
1129def CheckEquivalences(constrDict,varyList,parmDict=None,SeqHist=None):
1130    '''Process equivalence constraints, looking for conflicts such as
1131    where a parameter is used in both an equivalence and a constraint expression
1132    or where chaining is done (A->B and B->C).
1133    When called during refinements, parmDict is defined, and for sequential refinement
1134    SeqHist ia also defined.
1135
1136      * parmDict is used to remove equivalences where a parameter is not present
1137        in a refinement
1138      * SeqHist is used to rename wild-card parameter names in sequential
1139        refinements to use the current histogram.
1140    '''
1141    global dependentParmList,arrayList,invarrayList,indParmList,consNum
1142    global problemVars
1143    warnmsg = ''
1144    errmsg = ''
1145    problemVars = []
1146    # process fixed parameters (holds)
1147    fixVlist = [] # list of Fixed vars
1148    constrVars = [] # list of vars in constraint expressions
1149    for cdict in constrDict:
1150        # N.B. No "_" names in holds
1151        if len(cdict) == 1:
1152            fixVlist.append(list(cdict.keys())[0])
1153        else:
1154            constrVars += cdict.keys() # this will include _vary (not a problem)
1155    # process equivalences: make a list of dependent and independent vars
1156    #    and check for repeated uses (repetition of a parameter as an
1157    #    independent var is OK)
1158    indepVarList = []
1159    depVarList = []
1160    multdepVarList = []
1161    dropVarList = []
1162    translateTable = {} # lookup table for wildcard referenced parameters
1163    for varlist,mapvars,multarr,invmultarr in zip(
1164        dependentParmList,indParmList,arrayList,invarrayList):
1165        if multarr is None: # an equivalence
1166            zeromult = False
1167            for i,mv in enumerate(mapvars):
1168                if mv.split(':')[1] == '*' and SeqHist is not None:
1169                    # convert wildcard var to reference current histogram; save translation in table
1170                    sv = mv.split(':')
1171                    sv[1] = str(SeqHist)
1172                    mv = translateTable[mv] = ':'.join(sv)
1173                    mapvars[i] = mv
1174                varied = 0
1175                notvaried = ''
1176                if mv in varyList:
1177                    varied += 1
1178                else:
1179                    if notvaried: notvaried += ', '
1180                    notvaried += mv
1181                if parmDict is not None and mv not in parmDict:
1182                    print ("Dropping equivalence for parameter "+str(mv)+". Not defined in this refinement")
1183                    if mv not in dropVarList: dropVarList.append(mv)
1184                if mv not in indepVarList: indepVarList.append(mv)
1185            for i,(v,m) in enumerate(zip(varlist,invmultarr)):
1186                if v.split(':')[1] == '*' and SeqHist is not None:
1187                    # convert wildcard var to reference current histogram; save translation in table
1188                    sv = v.split(':')
1189                    sv[1] = str(SeqHist)
1190                    varlist[i] = v = translateTable[v] = ':'.join(sv)
1191                if parmDict is not None and v not in parmDict:
1192                    print ("Dropping equivalence for dep. variable "+str(v)+". Not defined in this refinement")
1193                    if v not in dropVarList: dropVarList.append(v)
1194                    continue
1195                if m == 0: zeromult = True
1196                if v in varyList:
1197                    varied += 1
1198                else:
1199                    if notvaried: notvaried += ', '
1200                    notvaried += v
1201                if v in indepVarList:
1202                    errmsg += '\nParameter '+v+' is used to set values in a constraint before its value is set in another constraint\n'
1203                    if v not in problemVars: problemVars.append(v)
1204                if v in depVarList:
1205                    multdepVarList.append(v)
1206                else:
1207                    depVarList.append(v)
1208            if varied > 0 and varied != len(varlist)+1:
1209                warnmsg += "\nNot all parameters refined in equivalence:\n\t"
1210                s = ""
1211                for v in varlist:
1212                    if s != "": s+= " & "
1213                    s += str(v)           
1214                warnmsg += str(mv) + " => " + s
1215                warnmsg += '\nNot refined: ' + notvaried + '\n'
1216            if zeromult:
1217                errmsg += "\nZero multiplier is invalid in equivalence:\n\t"
1218                s = ""
1219                for v in varlist:
1220                    if s != "": s+= " & "
1221                    s += str(v)           
1222                errmsg += str(mv) + " => " + s + '\n'
1223    # check for errors:
1224    if len(multdepVarList) > 0:
1225        errmsg += "\nThe following parameters(s) are used in conflicting Equivalence relations as dependent variables:\n"
1226        s = ''
1227        for var in sorted(set(multdepVarList)):
1228            if v not in problemVars: problemVars.append(v)
1229            if s != "": s+= ", "
1230            s += str(var)           
1231        errmsg += '\t'+ s + '\n'
1232    equivVarList = list(set(indepVarList).union(set(depVarList)))
1233    if debug: print ('equivVarList',equivVarList)
1234    # check for parameters that are both fixed and in an equivalence (not likely)
1235    inboth = set(fixVlist).intersection(set(equivVarList))
1236    if len(inboth) > 0:
1237        errmsg += "\nThe following parameter(s) are used in both Equivalence and Fixed constraints:\n"
1238        s = ''
1239        for var in sorted(inboth):
1240            if var not in problemVars: problemVars.append(var)
1241            if s != "": s+= ", "
1242            s += str(var)
1243        errmsg += '\t'+ s + '\n'
1244    # check for parameters that in both an equivalence and a constraint expression
1245    inboth = set(constrVars).intersection(set(equivVarList))
1246    if len(inboth) > 0:
1247        errmsg += "\nThe following parameter(s) are used in both Equivalence and Equiv or new var constraints:\n"
1248        s = ''
1249        for var in sorted(inboth):
1250            if var not in problemVars: problemVars.append(var)
1251            if s != "": s+= ", "
1252            s += str(var)
1253        errmsg += '\t'+ s + '\n'
1254    return errmsg,warnmsg,fixVlist,dropVarList,translateTable
1255
1256def MoveConfEquiv(constrDict,fixedList):
1257    '''Address conflicts in Equivalence constraints by creating an constraint equation
1258    that has the same action as the equivalence and removing the Equivalence
1259    '''
1260    global dependentParmList,arrayList,invarrayList,indParmList,consNum
1261    global problemVars
1262    parmsChanged = 0
1263    for i,(varlist,mapvars) in enumerate(zip(dependentParmList,indParmList)):
1264        conf = False
1265        for mv in mapvars:
1266            if mv in problemVars:
1267                conf = True
1268                break
1269        for v in varlist:
1270            if v in problemVars:
1271                conf = True
1272                break
1273        if conf:
1274            parmsChanged += 1
1275            indvar = indParmList[i][0]
1276            for dep,mult in zip(dependentParmList[i],invarrayList[i]):
1277                #print('debug replacing equiv with constraint equation 0=',{indvar:-1.,dep:mult[0]})
1278                constrDict += [{indvar:-1.,dep:mult[0]}]
1279                fixedList += ['0.0']
1280            dependentParmList[i] = None
1281    if parmsChanged:
1282        for i in range(len(dependentParmList)-1,-1,-1):
1283            if dependentParmList[i] is None:
1284                del dependentParmList[i],indParmList[i],arrayList[i],invarrayList[i]
1285    return parmsChanged
1286
1287def StoreEquivalence(independentVar,dependentList,symGen=True):
1288    '''Takes a list of dependent parameter(s) and stores their
1289    relationship to a single independent parameter (independentVar).
1290
1291    Called with user-supplied constraints by :func:`GSASIIstrIO.ProcessConstraints,
1292    with Pawley constraints from :func:`GSASIIstrIO.GetPawleyConstr`,
1293    with Unit Cell constraints from :func:`GSASIIstrIO.cellVary`
1294    with symmetry-generated atom constraints from :func:`GSASIIstrIO.GetPhaseData`
1295
1296    :param str independentVar: name of master parameter that will be used to determine the value
1297      to set the dependent variables
1298
1299    :param list dependentList: a list of parameters that will set from
1300         independentVar. Each item in the list can be a string with the parameter
1301         name or a tuple containing a name and multiplier:
1302         ``['::parm1',('::parm2',.5),]``
1303
1304    '''
1305   
1306    global dependentParmList,arrayList,invarrayList,indParmList,symGenList
1307    mapList = []
1308    multlist = []
1309    allfloat = True
1310    for var in dependentList:
1311        if isinstance(var, str):
1312            mult = 1.0
1313        elif len(var) == 2:
1314            var,mult = var
1315        else:
1316            raise Exception("Cannot parse "+repr(var) + " as var or (var,multiplier)")
1317        mapList.append(var)
1318        try:           
1319            multlist.append(tuple((float(mult),)))
1320        except:
1321            allfloat = False
1322            multlist.append(tuple((mult,)))
1323    # added relationships to stored values
1324    arrayList.append(None)
1325    if allfloat:
1326        invarrayList.append(np.array(multlist))
1327    else:
1328        invarrayList.append(multlist)
1329    indParmList.append(list((independentVar,)))
1330    dependentParmList.append(mapList)
1331    symGenList.append(symGen)
1332    return
1333
1334def EvaluateMultipliers(constList,*dicts):
1335    '''Convert multipliers for constraints and equivalences that are specified
1336    as strings into values. The strings can specify values in the parameter dicts as
1337    well as normal Python functions, such as "2*np.cos(0::Ax:2/2.)"
1338   
1339    :param list constList: a list of dicts containing constraint expressions
1340    :param \*dicts: one or more dicts containing GSAS-II parameters and their values
1341       can be specified
1342    :returns: an empty string if there were no errors, or an error message listing
1343       the strings that could not be converted.
1344    '''
1345    def SubfromParmDict(s,prmDict):
1346        for key in prmDict:
1347            if key in s:
1348                s = s.replace(key,str(prmDict[key]))
1349        return eval(s)
1350    prmDict = {}
1351    for d in dicts: prmDict.update(d) # combine all passed parameter dicts
1352    problemList = ""
1353    # loop through multipliers in contraint expressions
1354    for const in constList:
1355        for key in const:
1356            if key.startswith('_'): continue
1357            try: # is this already a float, etc?
1358                1+const[key]
1359                continue
1360            except:
1361                pass
1362            try:
1363                newval = SubfromParmDict(const[key][:],prmDict)
1364                if GSASIIpath.GetConfigValue('debug'):
1365                    print('Changing ',const[key],'to',newval)
1366                const[key] = newval               
1367            except:
1368                if problemList: problemList += ", "
1369                problemList += const[key]
1370    # loop through multipliers in equivalences
1371    global arrayList,invarrayList
1372    for i,(a,valList) in enumerate(zip(arrayList,invarrayList)):
1373        if a is not None: continue # ignore if not equiv
1374        try:
1375            valList.shape
1376            continue # ignore if already a numpy array
1377        except:
1378            pass
1379        repList = []
1380        for v in valList:
1381            try:
1382                1+v[0]
1383                repList.append(tuple((v[0],)))
1384                continue
1385            except:
1386                pass
1387            try:
1388                newval = SubfromParmDict(v[0][:],prmDict)
1389                if GSASIIpath.GetConfigValue('debug'):
1390                    print('Changing ',v[0],'to',newval)
1391                repList.append(tuple((newval,)))
1392            except:
1393                if problemList: problemList += ", "
1394                problemList += v[0]
1395                repList.append(tuple(('error',)))
1396        invarrayList[i] = np.array(repList)
1397    return problemList
1398
1399def GetDependentVars():
1400    '''Return a list of dependent variables: e.g. parameters that are
1401    constrained in terms of other parameters
1402
1403    :returns: a list of parameter names
1404
1405    '''
1406    global dependentVars
1407    return dependentVars
1408
1409def GetIndependentVars():
1410    '''Return a list of independent variables: e.g. parameters that are
1411    slaved to other parameters by constraints
1412
1413    :returns: a list of parameter names
1414
1415    '''
1416    global independentVars
1417    return independentVars
1418
1419def PrintIndependentVars(parmDict,varyList,sigDict,PrintAll=False,pFile=None):
1420    '''Print the values and uncertainties on the independent parameters'''
1421    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
1422    printlist = []
1423    mapvars = GetIndependentVars()
1424    for i,name in enumerate(mapvars):
1425        if name in fixedDict: continue
1426        if PrintAll or name in varyList:
1427            sig = sigDict.get(name)
1428            printlist.append([name,parmDict[name],sig])
1429    if len(printlist) == 0: return
1430    s1 = ''
1431    s2 = ''
1432    s3 = ''
1433    pFile.write(130*'-'+'\n')
1434    pFile.write("Parameters generated by constraints\n")
1435    printlist.append(3*[None])
1436    for name,val,esd in printlist:
1437        if len(s1) > 120 or name is None:
1438            pFile.write(''+'\n')
1439            pFile.write(s1+'\n')
1440            pFile.write(s2+'\n')
1441            pFile.write(s3+'\n')
1442            s1 = ''
1443            if name is None: break
1444        if s1 == "":
1445            s1 = ' name  :'
1446            s2 = ' value :'
1447            s3 = ' sig   :'
1448        s1 += '%15s' % (name)
1449        s2 += '%15.5f' % (val)
1450        if esd is None:
1451            s3 += '%15s' % ('n/a')
1452        else:
1453            s3 += '%15.5f' % (esd)
1454
1455def ComputeDepESD(covMatrix,varyList,parmDict):
1456    '''Compute uncertainties for dependent parameters from independent ones
1457    returns a dictionary containing the esd values for dependent parameters
1458    '''
1459    sigmaDict = {}
1460    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
1461        #if invmultarr is None: continue # probably not needed
1462#        try:
1463#            valuelist = [parmDict[var] for var in mapvars]
1464#        except KeyError:
1465#            continue
1466        # get the v-covar matrix for independent parameters
1467        vcov = np.zeros((len(mapvars),len(mapvars)))
1468        for i1,name1 in enumerate(mapvars):
1469            if name1 not in varyList: continue
1470            iv1 = varyList.index(name1)
1471            for i2,name2 in enumerate(mapvars):
1472                if name2 not in varyList: continue
1473                iv2 = varyList.index(name2)
1474                vcov[i1][i2] = covMatrix[iv1][iv2]
1475        # vec is the vector that multiplies each of the independent values
1476        for v,vec in zip(varlist,invmultarr):
1477            sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec)))
1478    return sigmaDict
1479
1480def _FormatConstraint(RelDict,RelVal):
1481    '''Formats a Constraint or Function for use in a convenient way'''
1482    linelen = 45
1483    s = [""]
1484    for var,val in RelDict.items():
1485        if var.startswith('_'): continue
1486        if len(s[-1]) > linelen: s.append(' ')
1487        m = val
1488        if s[-1] != "" and m >= 0:
1489            s[-1] += ' + '
1490        elif s[-1] != "":
1491            s[-1] += ' - '
1492            m = abs(m)
1493        s[-1] += '%.3f*%s '%(m,var)
1494    if len(s[-1]) > linelen: s.append(' ')
1495    if RelVal is None:
1496        s[-1] += ' = New variable'
1497    else:
1498        s[-1] += ' = ' + RelVal
1499    s1 = ''
1500    for s2 in s:
1501        if s1 != '': s1 += '\n\t'
1502        s1 += s2
1503    return s1
1504
1505def VarRemapShow(varyList,inputOnly=False):
1506    '''List out the saved relationships. This should be done after the constraints have been
1507    defined using :func:`StoreEquivalence`, :func:`GroupConstraints` and :func:`GenerateConstraints`.
1508
1509    :returns: a string containing the details of the contraint relationships
1510    '''
1511    s = ''
1512    if len(fixedVarList) > 0:
1513        s += 'Fixed Parameters:\n'
1514        for v in fixedVarList:
1515            s += '    ' + v + '\n'
1516    if not inputOnly:
1517        s += 'User-supplied parameter mapping relations:\n'
1518    symout = ''
1519    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,symGenList
1520
1521    for varlist,mapvars,multarr,invmultarr,symFlag in zip(
1522        dependentParmList,indParmList,arrayList,invarrayList,symGenList):
1523        for i,mv in enumerate(mapvars):
1524            if multarr is None:
1525#                s1 = '  ' + str(mv) + ' is equivalent to parameter(s): '
1526                if len(varlist) == 1:
1527                    s1 = '   ' + str(mv) + ' is equivalent to '
1528                else:
1529                    s1 = '   ' + str(mv) + ' is equivalent to parameters: '
1530                j = 0
1531                for v,m in zip(varlist,invmultarr):
1532                    if debug: print ('v,m[0]: ',v,m[0])
1533                    if len(s1.split('\n')[-1]) > 75: s1 += '\n        '
1534                    if j > 0: s1 += ' & '
1535                    j += 1
1536                    s1 += str(v)
1537                    if m != 1:
1538                        s1 += " / " + str(m[0])
1539                if symFlag:
1540                    symout += s1 + '\n'
1541                else:
1542                    s += s1 + '\n'
1543                continue
1544            s += %s = ' % mv
1545            j = 0
1546            for m,v in zip(multarr[i,:],varlist):
1547                if m == 0: continue
1548                if j > 0: s += ' + '
1549                j += 1
1550                s += '(%s * %s)' % (m,v)
1551            if mv in varyList: s += ' VARY'
1552            s += '\n'
1553    if symout:
1554        s += 'Symmetry-generated relations:\n' + symout
1555    if inputOnly: return s
1556    s += 'Inverse parameter mapping relations:\n'
1557    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
1558        for i,mv in enumerate(varlist):
1559            s += %s = ' % mv
1560            j = 0
1561            for m,v in zip(invmultarr[i,:],mapvars):
1562                if m == 0: continue
1563                if j > 0: s += ' + '
1564                j += 1
1565                s += '(%s * %s)' % (m,v)
1566            s += '\n'
1567    return s
1568
1569def GetSymEquiv():
1570    '''Return the automatically generated (equivalence) relationships.
1571
1572    :returns: a list of strings containing the details of the contraint relationships
1573    '''
1574    symout = []
1575    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,symGenList
1576
1577    for varlist,mapvars,multarr,invmultarr,symFlag in zip(
1578        dependentParmList,indParmList,arrayList,invarrayList,symGenList):
1579        for i,mv in enumerate(mapvars):
1580            if not symFlag: continue
1581            if multarr is None:
1582                #s1 = str(mv) + ' = '
1583                s1 = ''
1584                s2 = ' = ' + str(mv)
1585                j = 0
1586                if len(varlist) == 1:
1587                    # format the way Bob prefers
1588                    if invmultarr[0][0] == 1: 
1589                        s1 = str(varlist[0]) + ' = ' + str(mv)
1590                    else:
1591                        s1 = str(varlist[0]) + ' = ' + str(
1592                            invmultarr[0][0]) + ' * '+ str(mv)
1593                    symout.append(s1)
1594                    continue
1595                for v,m in zip(varlist,invmultarr):
1596                    if debug: print ('v,m[0]: ',v,m[0])
1597                    if len(s1.split('\n')[-1]) > 75: s1 += '\n        '
1598                    if j > 0: s1 += ' =  '
1599                    j += 1
1600                    s1 += str(v)
1601                    if m != 1:
1602                        s1 += " / " + str(m[0])
1603                symout.append(s1+s2)
1604                continue
1605            else:
1606                s = %s = ' % mv
1607                j = 0
1608                for m,v in zip(multarr[i,:],varlist):
1609                    if m == 0: continue
1610                    if j > 0: s += ' + '
1611                    j += 1
1612                    s += '(%s * %s)' % (m,v)
1613                print ('unexpected sym op='+s)
1614    return symout
1615
1616def Dict2Deriv(varyList,derivDict,dMdv):
1617    '''Compute derivatives for Independent Parameters from the
1618    derivatives for the original parameters
1619
1620    :param list varyList: a list of parameters names that will be varied
1621
1622    :param dict derivDict: a dict containing derivatives for parameter values keyed by the
1623      parameter names.
1624
1625    :param list dMdv: a Jacobian, as a list of np.array containing derivatives for dependent
1626      parameter computed from derivDict
1627
1628    '''
1629    global dependentParmList,arrayList,invarrayList,indParmList,invarrayList
1630    for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
1631        for i,name in enumerate(mapvars):
1632            # grouped parameters: need to add in the derv. w/r
1633            # dependent variables to the independent ones
1634            if name not in varyList: continue # skip if independent var not varied
1635            if multarr is None:
1636                for v,m in zip(varlist,invmultarr):
1637                    if debug: print ('start dMdv',dMdv[varyList.index(name)])
1638                    if debug: print ('add derv',v,'/',m[0],'to derv',name,'add=',derivDict[v] / m[0])
1639                    if m == 0: continue
1640                    dMdv[varyList.index(name)] += derivDict[v] / m[0]
1641            else:
1642                for v,m in zip(varlist,multarr[i,:]):
1643                    if debug: print ('start dMdv',dMdv[varyList.index(name)])
1644                    if debug: print ('add derv',v,'*',m,'to derv',name,'add=',m * derivDict[v])
1645                    if m == 0: continue
1646                    dMdv[varyList.index(name)] += m * derivDict[v]
1647
1648def Map2Dict(parmDict,varyList):
1649    '''Create (or update) the Independent Parameters from the original
1650    set of Parameters
1651
1652    Removes dependent variables from the varyList
1653
1654    This should be done once, after the constraints have been
1655    defined using :func:`StoreEquivalence`,
1656    :func:`GroupConstraints` and :func:`GenerateConstraints` and
1657    before any parameter refinement is done. This completes the parameter
1658    dictionary by defining independent parameters and it satisfies the
1659    constraint equations in the initial parameters
1660
1661    :param dict parmDict: a dict containing parameter values keyed by the
1662      parameter names.
1663      This will contain updated values for both dependent and independent
1664      parameters after Dict2Map is called. It will also contain some
1665      unexpected entries of every constant value {'0':0.0} & {'1.0':1.0},
1666      which do not cause any problems.
1667
1668    :param list varyList: a list of parameters names that will be varied
1669   
1670
1671    '''
1672    # process the Independent vars: remove dependent ones from varylist
1673    # and then compute values for the Independent ones from their dependents
1674    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
1675    for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
1676        for item in varlist:
1677            try:
1678                varyList.remove(item)
1679            except ValueError:
1680                pass
1681        if multarr is None: continue
1682        valuelist = [parmDict[var] for var in varlist]
1683        parmDict.update(zip(mapvars,
1684                            np.dot(multarr,np.array(valuelist)))
1685                        )
1686    # now remove fixed parameters from the varyList
1687    global fixedVarList
1688    for item in fixedVarList:
1689        try:
1690            varyList.remove(item)
1691        except ValueError:
1692            pass
1693    # Set constrained parameters to their fixed values
1694    parmDict.update(fixedDict)
1695
1696def Dict2Map(parmDict,varyList):
1697    '''Applies the constraints defined using :func:`StoreEquivalence`,
1698    :func:`GroupConstraints` and :func:`GenerateConstraints` by changing
1699    values in a dict containing the parameters. This should be
1700    done before the parameters are used for any computations
1701
1702    :param dict parmDict: a dict containing parameter values keyed by the
1703      parameter names.
1704      This will contain updated values for both dependent and independent
1705      parameters after Dict2Map is called. It will also contain some
1706      unexpected entries of every constant value {'0':0.0} & {'1.0':1.0},
1707      which do not cause any problems.
1708
1709    :param list varyList: a list of parameters names that will be varied
1710   
1711    '''
1712    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
1713    # reset fixed values (should not be needed, but very quick)
1714    # - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended
1715    # not needed, but also not a problem - BHT
1716    parmDict.update(fixedDict)
1717    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
1718        #if invmultarr is None: continue
1719        try: 
1720            valuelist = [parmDict[var] for var in mapvars]
1721        except KeyError:
1722            continue
1723        parmDict.update(zip(varlist,np.dot(invmultarr,np.array(valuelist))))
1724
1725#======================================================================
1726# internal routines follow (these routines are unlikely to be called
1727# from outside the module)
1728
1729def GramSchmidtOrtho(a,nkeep=0):
1730    '''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to
1731    find orthonormal unit vectors relative to first row.
1732
1733    If nkeep is non-zero, the first nkeep rows in the array are not changed
1734   
1735    input:
1736       arrayin: a 2-D non-singular square array
1737    returns:
1738       a orthonormal set of unit vectors as a square array
1739    '''
1740    def proj(a,b):
1741        'Projection operator'
1742        return a*(np.dot(a,b)/np.dot(a,a))
1743    for j in range(nkeep,len(a)):
1744        for i in range(j):
1745            a[j] -= proj(a[i],a[j])
1746        if np.allclose(np.linalg.norm(a[j]),0.0):
1747            raise ConstraintException("Singular input to GramSchmidtOrtho")
1748        a[j] /= np.linalg.norm(a[j])
1749    return a
1750
1751def _FillArray(sel,dict,collist,FillDiagonals=False):
1752    '''Construct a n by n matrix (n = len(collist)
1753    filling in the rows using the relationships defined
1754    in the dictionaries selected by sel
1755
1756    If FillDiagonals is True, diagonal elements in the
1757    array are set to 1.0
1758    '''
1759    n = len(collist)
1760    if FillDiagonals:
1761        arr = np.eye(n)
1762    else:
1763        arr = np.zeros(2*[n,])
1764    # fill the top rows
1765    for i,cnum in enumerate(sel):
1766        for j,var in enumerate(collist):
1767            arr[i,j] = dict[cnum].get(var,0)
1768    return arr
1769
1770def _SwapColumns(i,m,v):
1771    '''Swap columns in matrix m as well as the labels in v
1772    so that element (i,i) is replaced by the first non-zero element in row i after that element
1773
1774    Throws an exception if there are no non-zero elements in that row
1775    '''
1776    for j in range(i+1,len(v)):
1777        if not np.allclose(m[i,j],0):
1778            m[:,(i,j)] = m[:,(j,i)]
1779            v[i],v[j] = v[j],v[i]
1780            return
1781    else:
1782        raise ConstraintException('Singular input')
1783
1784def _RowEchelon(m,arr,collist):
1785    '''Convert the first m rows in Matrix arr to row-echelon form
1786    exchanging columns in the matrix and collist as needed.
1787
1788    throws an exception if the matrix is singular because
1789    the first m rows are not linearly independent
1790    '''
1791    for i in range(m):
1792        if np.allclose(arr[i,i],0):
1793            _SwapColumns(i,arr,collist)
1794        arr[i,:] /= arr[i,i] # normalize row
1795        # subtract current row from subsequent rows to set values to left of diagonal to 0
1796        for j in range(i+1,m):
1797            arr[j,:] -= arr[i,:] * arr[j,i]
1798
1799if __name__ == "__main__":
1800    parmdict = {}
1801    constrDict = [
1802        {'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},
1803        {'0:0:eA': 0.0},
1804        {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
1805        {'1::C(10,0,1)': 1.0, '2::C(10,0,1)': 1.0},
1806        {'1::AUiso:0': 1.0, '0::AUiso:0': 1.0},
1807        {'0::A0': 0.0}
1808        ]
1809    fixedList = ['5.0', '0', None, None, '1.0', '0']
1810    StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), ))
1811    #StoreEquivalence('1::atomx:3',('2::atomx:3', ('2::atomz:3',2,), )) # error: dependent & independent vars mixed
1812    #StoreEquivalence('1::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) # error: dependent vars repeated
1813    #StoreEquivalence('0:1:eA',('0:0:eA',)) # error: equiv & fixed
1814    #StoreEquivalence('0:99:Scale',('0:12:Scale',)) # error: equiv & constrained
1815    #StoreEquivalence('0:12:Scale',('0:99:Scale',)) # error: equiv & constrained
1816    varylist = ['2::atomx:3',
1817                '2::C(10,6,1)', '1::C(10,6,1)',
1818                '2::atomy:3', '2::atomz:3',
1819                '0:12:Scale', '0:11:Scale', '0:14:Scale', '0:13:Scale', '0:0:Scale']
1820#    e,w = CheckConstraints([,
1821#                     [{'2:0:Scale': 1.0, '5:0:Scale': 1.0, '10:0:Scale': 1.0, '6:0:Scale': 1.0, '9:0:Scale': 1.0, '8:0:Scale': 1.0,# '3:0:Scale': 1.0, '4:0:Scale': 1.0, '7:0:Scale': 1.0, '1:0:Scale': 1.0, '0:0:Scale': 1.0}],
1822#                     ['1.0'])
1823#    if e: print 'error=',e
1824#    if w: print 'error=',w
1825#    varyList = ['0::A0', '0::AUiso:0', '0::Afrac:1', '0::Afrac:2', '0::Afrac:3', '0::Afrac:4',
1826#       '0::dAx:5', '0::dAy:5', '0::dAz:5', '0::AUiso:5', ':0:Back;0', ':0:Back;1', ':0:Back;2', ':0:Back;3',
1827#       ':0:Back;4', ':0:Back;5', ':0:Back;6', ':0:Back;7', ':0:Back;8', ':0:Back;9', ':0:Back;10', ':0:Back;11'
1828#       :0:U', ':0:V', ':0:W', ':0:X', ':0:Y', ':0:Scale', ':0:DisplaceX', ':0:DisplaceY']
1829#    constrDict = [
1830#        {'0::Afrac:4': 24.0, '0::Afrac:1': 16.0, '0::Afrac:3': 24.0, '0::Afrac:2': 16.0},
1831#        {'0::Afrac:1': 1.0, '0::Afrac:2': 1.0},
1832#        {'0::Afrac:4': 1.0, '0::Afrac:3': 1.0}]
1833#    fixedList = ['40.0', '1.0', '1.0']
1834
1835    errmsg, warnmsg = CheckConstraints(varylist,constrDict,fixedList)
1836    if errmsg:
1837        print ("*** Error ********************")
1838        print (errmsg)
1839    if warnmsg:
1840        print ("*** Warning ********************")
1841        print (warnmsg)
1842    if errmsg or warnmsg:
1843        sys.exit()
1844    groups,parmlist = GroupConstraints(constrDict)
1845    GenerateConstraints(groups,parmlist,varylist,constrDict,fixedList)
1846    print (VarRemapShow(varylist))
1847    parmdict.update( {
1848        '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,
1849        '0:0:eA': 0.0,
1850        '2::C(10,6,1)': 0.2, '1::C(10,6,1)': 0.3,
1851        '1::C(10,0,1)': 0.2, '2::C(10,0,1)': 0.3,
1852        '1::AUiso:0': 0.02, '0::AUiso:0': 0.03,
1853        '0::A0': 0.0,
1854        '2::atomx:3':0.23,'2::atomy:3':-.23, '2::atomz:3':-0.11,
1855        })
1856    print ('parmdict start',parmdict)
1857    print ('varylist start',varylist)
1858    before = parmdict.copy()
1859    Map2Dict(parmdict,varylist)
1860    print ('parmdict before and after Map2Dict')
1861    print ('  key / before / after')
1862    for key in sorted(list(parmdict.keys())):
1863        print ('  '+key,'\t',before.get(key),'\t',parmdict[key])
1864    print ('varylist after',varylist)
1865    before = parmdict.copy()
1866    Dict2Map(parmdict,varylist)
1867    print ('after Dict2Map')
1868    print ('  key / before / after')
1869    for key in sorted(list(parmdict.keys())):
1870        print ('  '+key,'\t',before.get(key),'\t',parmdict[key])
1871#    dMdv = len(varylist)*[0]
1872#    deriv = {}
1873#    for i,v in enumerate(parmdict.keys()): deriv[v]=i
1874#    Dict2Deriv(varylist,deriv,dMdv)
Note: See TracBrowser for help on using the repository browser.