source: trunk/GSASIImapvars.py @ 5000

Last change on this file since 5000 was 4998, checked in by toby, 3 months ago

partial fix of magnetic constraints (still need to deal w/unvaried items; fullrmc minor stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 80.5 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2021-07-21 20:04:47 +0000 (Wed, 21 Jul 2021) $
4# $Author: toby $
5# $Revision: 4998 $
6# $URL: trunk/GSASIImapvars.py $
7# $Id: GSASIImapvars.py 4998 2021-07-21 20:04:47Z 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: 4998 $")
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,mvMsg = 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,mvMsg = 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\n\t"+_FormatConstraint(constrDict[rel],fixedList[rel])+'\n'
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\n\t"+_FormatConstraint(constrDict[rel],fixedList[rel])+'\n'
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 found while processing constraints group:'
946            for rel in group:
947                msg += '\n  ' + _FormatConstraint(constrDict[rel],fixedList[rel])
948            break
949        mapvar = []
950        group = group[:]
951        # scan through all generated and input relationships, we need to add to the varied list
952        # all the new parameters where VaryFree has been set or where a New Var is varied.
953        #
954        # If a group does not contain any fixed values (constraint equations)
955        # and nothing in the group is varied, drop this group, so that the
956        # dependent parameters can be refined individually.
957        unused = True
958        for i in range(len(varlist)):
959            if len(group) > 0: # get the original equation reference
960                rel = group.pop(0)
961                fixedval = fixedList[rel]
962                varyflag = constrDict[rel].get('_vary',False)
963                varname = constrDict[rel].get('_name','')
964            else: # this relationship has been generated
965                varyflag = False
966                varname = ''
967                fixedval = None
968            if fixedval is None: # this is a new parameter, not a constraint
969                if not varname:
970                    varname = paramPrefix + str(consNum) # no assigned name, create one
971                    consNum += 1
972                mapvar.append(varname)
973                genVarLookup[varname] = varlist # save list of parameters related to this new var
974                # vary the new relationship if it is a degree of freedom in
975                # a set of contraint equations or if a New Var is flagged to be varied.
976                if VaryFree or varyflag: 
977                    unused = False
978                    varyList.append(varname)
979                    # fix (prevent varying) of all the parameters inside the constraint group
980                    # (dependent vars)
981                    for var in varsList:
982                        if var in varyList: varyList.remove(var)
983            else:
984                unused = False
985                mapvar.append(fixedval)
986        if unused: # skip over constraints that don't matter (w/o fixed value or any refined parameters)
987            if debug: print('Constraint ignored (all parameters unrefined)')
988            if debug: print ('   '+_FormatConstraint(constrDict[rel],fixedList[rel]))
989            continue 
990        dependentParmList.append([translateTable.get(var,var) for var in varlist])
991        arrayList.append(constrArr)
992        invarrayList.append(np.linalg.inv(constrArr))
993        indParmList.append(mapvar)
994        symGenList.append(False)
995    if msg:
996        print (' *** ERROR in constraint definitions! ***')
997        print (msg)
998        print (VarRemapShow(varyList))
999        raise ConstraintException
1000    # setup dictionary containing the fixed values
1001    global fixedDict
1002    # key is original ascii string, value is float
1003    for fixedval in fixedList:
1004        if fixedval:
1005            fixedDict[fixedval] = float(fixedval)
1006    _setVarLists(dropVarList)
1007    if changed:
1008        print(60*'=')
1009        print('Constraints were reclassified to avoid conflicts, as below:')
1010        print(mvMsg)
1011        print('New constraints are:')
1012        print (VarRemapShow(varyList,True))
1013        print(60*'=')
1014    return groups,parmlist # saved for sequential fits
1015   
1016def _setVarLists(dropVarList):
1017    '''Make list of dependent and independent variables (after dropping unused vars in dropVarList)
1018    '''
1019    global dependentParmList,indParmList
1020    global dependentVars
1021    global independentVars
1022    dependentVars = []
1023    independentVars = []
1024    for varlist,mapvars in zip(dependentParmList,indParmList):  # process all constraints
1025        for mv in mapvars:
1026            if mv in dropVarList: continue
1027            if mv not in independentVars: independentVars.append(mv)
1028        for mv in varlist:
1029            if mv in dropVarList: continue
1030            if mv not in dependentVars: dependentVars.append(mv)
1031    if debug: # on debug, show what is parsed & generated, semi-readable
1032        print (50*'-')
1033        #print (VarRemapShow(varyList))
1034        #print ('Varied: ',varyList)
1035        print ('Not Varied: ',fixedVarList)
1036
1037# def CheckEquivalences(constrDict,varyList):
1038#     global dependentParmList,arrayList,invarrayList,indParmList,consNum
1039#     global problemVars
1040#     warnmsg = ''
1041#     errmsg = ''
1042#     problemVars = []
1043#     # process fixed variables (holds)
1044#     fixVlist = [] # list of Fixed vars
1045#     constrVars = [] # list of vars in constraint expressions
1046#     for cdict in constrDict:
1047#         # N.B. No "_" names in holds
1048#         if len(cdict) == 1:
1049#             fixVlist.append(list(cdict.keys())[0])
1050#         else:
1051#             constrVars += cdict.keys() # this will include _vary (not a problem)
1052#     # process equivalences: make a list of dependent and independent vars
1053#     #    and check for repeated uses (repetition of a parameter as an
1054#     #    independent var is OK)
1055#     indepVarList = []
1056#     depVarList = []
1057#     multdepVarList = []
1058#     for varlist,mapvars,multarr,invmultarr in zip(
1059#         dependentParmList,indParmList,arrayList,invarrayList):
1060#         if multarr is None: # an equivalence
1061#             zeromult = False
1062#             for mv in mapvars:
1063#                 varied = 0
1064#                 notvaried = ''
1065#                 if mv in varyList:
1066#                     varied += 1
1067#                 else:
1068#                     if notvaried: notvaried += ', '
1069#                     notvaried += mv
1070#                 if mv not in indepVarList: indepVarList.append(mv)
1071#                 for v,m in zip(varlist,invmultarr):
1072#                     if v in indepVarList:
1073#                         errmsg += '\nVariable '+v+' is used to set values in a constraint before its value is set in another constraint\n'
1074#                         if v not in problemVars: problemVars.append(v)
1075#                     if m == 0: zeromult = True
1076#                     if v in varyList:
1077#                         varied += 1
1078#                     else:
1079#                         if notvaried: notvaried += ', '
1080#                         notvaried += v
1081#                     if v in depVarList:
1082#                         multdepVarList.append(v)
1083#                     else:
1084#                         depVarList.append(v)
1085#             if varied > 0 and varied != len(varlist)+1:
1086#                 warnmsg += "\nNot all variables refined in equivalence:\n\t"
1087#                 s = ""
1088#                 for v in varlist:
1089#                     if s != "": s+= " & "
1090#                     s += str(v)           
1091#                 warnmsg += str(mv) + " => " + s
1092#                 warnmsg += '\nNot refined: ' + notvaried + '\n'
1093#             if zeromult:
1094#                 errmsg += "\nZero multiplier is invalid in equivalence:\n\t"
1095#                 s = ""
1096#                 for v in varlist:
1097#                     if s != "": s+= " & "
1098#                     s += str(v)           
1099#                 errmsg += str(mv) + " => " + s + '\n'
1100#     # check for errors:
1101#     if len(multdepVarList) > 0:
1102#         errmsg += "\nThe following parameters(s) are used in conflicting Equivalence relations as dependent variables:\n"
1103#         s = ''
1104#         for var in sorted(set(multdepVarList)):
1105#             if v not in problemVars: problemVars.append(v)
1106#             if s != "": s+= ", "
1107#             s += str(var)           
1108#         errmsg += '\t'+ s + '\n'
1109#     equivVarList = list(set(indepVarList).union(set(depVarList)))
1110#     if debug: print ('equivVarList',equivVarList)
1111#     # check for parameters that are both fixed and in an equivalence (not likely)
1112#     inboth = set(fixVlist).intersection(set(equivVarList))
1113#     if len(inboth) > 0:
1114#         errmsg += "\nThe following parameter(s) are used in both Equivalence and Fixed constraints:\n"
1115#         s = ''
1116#         for var in sorted(inboth):
1117#             if var not in problemVars: problemVars.append(var)
1118#             if s != "": s+= ", "
1119#             s += str(var)
1120#         errmsg += '\t'+ s + '\n'
1121#     # check for parameters that in both an equivalence and a constraint expression
1122#     inboth = set(constrVars).intersection(set(equivVarList))
1123#     if len(inboth) > 0:
1124#         errmsg += "\nThe following parameter(s) are used in both Equivalence and Equiv or new var constraints:\n"
1125#         s = ''
1126#         for var in sorted(inboth):
1127#             if var not in problemVars: problemVars.append(var)
1128#             if s != "": s+= ", "
1129#             s += str(var)
1130#         errmsg += '\t'+ s + '\n'
1131#     return errmsg,warnmsg,fixVlist
1132
1133def CheckEquivalences(constrDict,varyList,parmDict=None,SeqHist=None):
1134    '''Process equivalence constraints, looking for conflicts such as
1135    where a parameter is used in both an equivalence and a constraint expression
1136    or where chaining is done (A->B and B->C).
1137    When called during refinements, parmDict is defined, and for sequential refinement
1138    SeqHist ia also defined.
1139
1140      * parmDict is used to remove equivalences where a parameter is not present
1141        in a refinement
1142      * SeqHist is used to rename wild-card parameter names in sequential
1143        refinements to use the current histogram.
1144    '''
1145    global dependentParmList,arrayList,invarrayList,indParmList,consNum
1146    global problemVars
1147    warnmsg = ''
1148    errmsg = ''
1149    problemVars = []
1150    # process fixed parameters (holds)
1151    fixVlist = [] # list of Fixed vars
1152    constrVars = [] # list of vars in constraint expressions
1153    for cdict in constrDict:
1154        # N.B. No "_" names in holds
1155        if len(cdict) == 1:
1156            fixVlist.append(list(cdict.keys())[0])
1157        else:
1158            constrVars += cdict.keys() # this will include _vary (not a problem)
1159    # process equivalences: make a list of dependent and independent vars
1160    #    and check for repeated uses (repetition of a parameter as an
1161    #    independent var is OK)
1162    indepVarList = []
1163    depVarList = []
1164    multdepVarList = []
1165    dropVarList = []
1166    translateTable = {} # lookup table for wildcard referenced parameters
1167    for varlist,mapvars,multarr,invmultarr in zip(
1168        dependentParmList,indParmList,arrayList,invarrayList):
1169        if multarr is None: # an equivalence
1170            zeromult = False
1171            for i,mv in enumerate(mapvars):
1172                if mv.split(':')[1] == '*' and SeqHist is not None:
1173                    # convert wildcard var to reference current histogram; save translation in table
1174                    sv = mv.split(':')
1175                    sv[1] = str(SeqHist)
1176                    mv = translateTable[mv] = ':'.join(sv)
1177                    mapvars[i] = mv
1178                varied = 0
1179                notvaried = ''
1180                if mv in varyList:
1181                    varied += 1
1182                else:
1183                    if notvaried: notvaried += ', '
1184                    notvaried += mv
1185                if parmDict is not None and mv not in parmDict:
1186                    print ("Dropping equivalence for parameter "+str(mv)+". Not defined in this refinement")
1187                    if mv not in dropVarList: dropVarList.append(mv)
1188                if mv not in indepVarList: indepVarList.append(mv)
1189            for i,(v,m) in enumerate(zip(varlist,invmultarr)):
1190                if v.split(':')[1] == '*' and SeqHist is not None:
1191                    # convert wildcard var to reference current histogram; save translation in table
1192                    sv = v.split(':')
1193                    sv[1] = str(SeqHist)
1194                    varlist[i] = v = translateTable[v] = ':'.join(sv)
1195                if parmDict is not None and v not in parmDict:
1196                    print ("Dropping equivalence for dep. variable "+str(v)+". Not defined in this refinement")
1197                    if v not in dropVarList: dropVarList.append(v)
1198                    continue
1199                if m == 0: zeromult = True
1200                if v in varyList:
1201                    varied += 1
1202                else:
1203                    if notvaried: notvaried += ', '
1204                    notvaried += v
1205                if v in indepVarList:
1206                    errmsg += '\nParameter '+v+' is used to set values in a constraint before its value is set in another constraint\n'
1207                    if v not in problemVars: problemVars.append(v)
1208                if v in depVarList:
1209                    multdepVarList.append(v)
1210                else:
1211                    depVarList.append(v)
1212            if varied > 0 and varied != len(varlist)+1:
1213                warnmsg += "\nNot all parameters refined in equivalence:\n\t"
1214                s = ""
1215                for v in varlist:
1216                    if s != "": s+= " & "
1217                    s += str(v)           
1218                warnmsg += str(mv) + " => " + s
1219                warnmsg += '\nNot refined: ' + notvaried + '\n'
1220            if zeromult:
1221                errmsg += "\nZero multiplier is invalid in equivalence:\n\t"
1222                s = ""
1223                for v in varlist:
1224                    if s != "": s+= " & "
1225                    s += str(v)           
1226                errmsg += str(mv) + " => " + s + '\n'
1227    # check for errors:
1228    if len(multdepVarList) > 0:
1229        errmsg += "\nThe following parameters(s) are used in conflicting Equivalence relations as dependent variables:\n"
1230        s = ''
1231        for var in sorted(set(multdepVarList)):
1232            if v not in problemVars: problemVars.append(v)
1233            if s != "": s+= ", "
1234            s += str(var)           
1235        errmsg += '\t'+ s + '\n'
1236    equivVarList = list(set(indepVarList).union(set(depVarList)))
1237    if debug: print ('equivVarList',equivVarList)
1238    # check for parameters that are both fixed and in an equivalence (not likely)
1239    inboth = set(fixVlist).intersection(set(equivVarList))
1240    if len(inboth) > 0:
1241        errmsg += "\nThe following parameter(s) are used in both Equivalence and Fixed constraints:\n"
1242        s = ''
1243        for var in sorted(inboth):
1244            if var not in problemVars: problemVars.append(var)
1245            if s != "": s+= ", "
1246            s += str(var)
1247        errmsg += '\t'+ s + '\n'
1248    # check for parameters that in both an equivalence and a constraint expression
1249    inboth = set(constrVars).intersection(set(equivVarList))
1250    if len(inboth) > 0:
1251        errmsg += "\nThe following parameter(s) are used in both Equivalence and Equiv or new var constraints:\n"
1252        s = ''
1253        for var in sorted(inboth):
1254            if var not in problemVars: problemVars.append(var)
1255            if s != "": s+= ", "
1256            s += str(var)
1257        errmsg += '\t'+ s + '\n'
1258    return errmsg,warnmsg,fixVlist,dropVarList,translateTable
1259
1260def MoveConfEquiv(constrDict,fixedList):
1261    '''Address conflicts in Equivalence constraints by creating an constraint equation
1262    that has the same action as the equivalence and removing the Equivalence
1263    '''
1264    global dependentParmList,arrayList,invarrayList,indParmList,consNum
1265    global problemVars
1266    parmsChanged = 0
1267    msg = ''
1268    if problemVars:
1269        msg = 'Conflict: variable(s) used in both equivalences and constraints: '
1270        for i1,v1 in enumerate(problemVars):
1271            if i1 > 0: msg += ', '
1272            msg += v1
1273    for i,(varlist,mapvars) in enumerate(zip(dependentParmList,indParmList)):
1274        conf = False
1275        for mv in mapvars:
1276            if mv in problemVars:
1277                conf = True
1278                break
1279        for v in varlist:
1280            if v in problemVars:
1281                conf = True
1282                break
1283        if conf:
1284            parmsChanged += 1
1285            indvar = indParmList[i][0]
1286            msg += '\n  Removing equivalence:\n    ' + _showEquiv(
1287                dependentParmList[i],indParmList[i],invarrayList[i])
1288            msg += '\n  Creating new constraint(s):'
1289            for dep,mult in zip(dependentParmList[i],invarrayList[i]):
1290                constrDict += [{indvar:-1.,dep:mult[0]}]
1291                fixedList += ['0.0']
1292                msg += '\n    ' + _FormatConstraint(constrDict[-1],fixedList[-1])
1293            dependentParmList[i] = None
1294    if parmsChanged:
1295        for i in range(len(dependentParmList)-1,-1,-1):
1296            if dependentParmList[i] is None:
1297                del dependentParmList[i],indParmList[i],arrayList[i],invarrayList[i]
1298    return parmsChanged,msg
1299
1300def StoreEquivalence(independentVar,dependentList,symGen=True):
1301    '''Takes a list of dependent parameter(s) and stores their
1302    relationship to a single independent parameter (independentVar).
1303
1304    Called with user-supplied constraints by :func:`GSASIIstrIO.ProcessConstraints,
1305    with Pawley constraints from :func:`GSASIIstrIO.GetPawleyConstr`,
1306    with Unit Cell constraints from :func:`GSASIIstrIO.cellVary`
1307    with symmetry-generated atom constraints from :func:`GSASIIstrIO.GetPhaseData`
1308
1309    :param str independentVar: name of master parameter that will be used to determine the value
1310      to set the dependent variables
1311
1312    :param list dependentList: a list of parameters that will set from
1313         independentVar. Each item in the list can be a string with the parameter
1314         name or a tuple containing a name and multiplier:
1315         ``['::parm1',('::parm2',.5),]``
1316
1317    '''
1318   
1319    global dependentParmList,arrayList,invarrayList,indParmList,symGenList
1320    mapList = []
1321    multlist = []
1322    allfloat = True
1323    for var in dependentList:
1324        if isinstance(var, str):
1325            mult = 1.0
1326        elif len(var) == 2:
1327            var,mult = var
1328        else:
1329            raise Exception("Cannot parse "+repr(var) + " as var or (var,multiplier)")
1330        mapList.append(var)
1331        try:           
1332            multlist.append(tuple((float(mult),)))
1333        except:
1334            allfloat = False
1335            multlist.append(tuple((mult,)))
1336    # added relationships to stored values
1337    arrayList.append(None)
1338    if allfloat:
1339        invarrayList.append(np.array(multlist))
1340    else:
1341        invarrayList.append(multlist)
1342    indParmList.append(list((independentVar,)))
1343    dependentParmList.append(mapList)
1344    symGenList.append(symGen)
1345    return
1346
1347def EvaluateMultipliers(constList,*dicts):
1348    '''Convert multipliers for constraints and equivalences that are specified
1349    as strings into values. The strings can specify values in the parameter dicts as
1350    well as normal Python functions, such as "2*np.cos(0::Ax:2/2.)"
1351   
1352    :param list constList: a list of dicts containing constraint expressions
1353    :param \*dicts: one or more dicts containing GSAS-II parameters and their values
1354       can be specified
1355    :returns: an empty string if there were no errors, or an error message listing
1356       the strings that could not be converted.
1357    '''
1358    def SubfromParmDict(s,prmDict):
1359        for key in prmDict:
1360            if key in s:
1361                s = s.replace(key,str(prmDict[key]))
1362        return eval(s)
1363    prmDict = {}
1364    for d in dicts: prmDict.update(d) # combine all passed parameter dicts
1365    problemList = ""
1366    # loop through multipliers in contraint expressions
1367    for const in constList:
1368        for key in const:
1369            if key.startswith('_'): continue
1370            try: # is this already a float, etc?
1371                1+const[key]
1372                continue
1373            except:
1374                pass
1375            try:
1376                newval = SubfromParmDict(const[key][:],prmDict)
1377                if GSASIIpath.GetConfigValue('debug'):
1378                    print('Changing ',const[key],'to',newval)
1379                const[key] = newval               
1380            except:
1381                if problemList: problemList += ", "
1382                problemList += const[key]
1383    # loop through multipliers in equivalences
1384    global arrayList,invarrayList
1385    for i,(a,valList) in enumerate(zip(arrayList,invarrayList)):
1386        if a is not None: continue # ignore if not equiv
1387        try:
1388            valList.shape
1389            continue # ignore if already a numpy array
1390        except:
1391            pass
1392        repList = []
1393        for v in valList:
1394            try:
1395                1+v[0]
1396                repList.append(tuple((v[0],)))
1397                continue
1398            except:
1399                pass
1400            try:
1401                newval = SubfromParmDict(v[0][:],prmDict)
1402                if GSASIIpath.GetConfigValue('debug'):
1403                    print('Changing ',v[0],'to',newval)
1404                repList.append(tuple((newval,)))
1405            except:
1406                if problemList: problemList += ", "
1407                problemList += v[0]
1408                repList.append(tuple(('error',)))
1409        invarrayList[i] = np.array(repList)
1410    return problemList
1411
1412def GetDependentVars():
1413    '''Return a list of dependent variables: e.g. parameters that are
1414    constrained in terms of other parameters
1415
1416    :returns: a list of parameter names
1417
1418    '''
1419    global dependentVars
1420    return dependentVars
1421
1422def GetIndependentVars():
1423    '''Return a list of independent variables: e.g. parameters that are
1424    slaved to other parameters by constraints
1425
1426    :returns: a list of parameter names
1427
1428    '''
1429    global independentVars
1430    return independentVars
1431
1432def PrintIndependentVars(parmDict,varyList,sigDict,PrintAll=False,pFile=None):
1433    '''Print the values and uncertainties on the independent parameters'''
1434    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
1435    printlist = []
1436    mapvars = GetIndependentVars()
1437    for i,name in enumerate(mapvars):
1438        if name in fixedDict: continue
1439        if PrintAll or name in varyList:
1440            sig = sigDict.get(name)
1441            printlist.append([name,parmDict[name],sig])
1442    if len(printlist) == 0: return
1443    s1 = ''
1444    s2 = ''
1445    s3 = ''
1446    pFile.write(130*'-'+'\n')
1447    pFile.write("Parameters generated by constraints\n")
1448    printlist.append(3*[None])
1449    for name,val,esd in printlist:
1450        if len(s1) > 120 or name is None:
1451            pFile.write(''+'\n')
1452            pFile.write(s1+'\n')
1453            pFile.write(s2+'\n')
1454            pFile.write(s3+'\n')
1455            s1 = ''
1456            if name is None: break
1457        if s1 == "":
1458            s1 = ' name  :'
1459            s2 = ' value :'
1460            s3 = ' sig   :'
1461        s1 += '%15s' % (name)
1462        s2 += '%15.5f' % (val)
1463        if esd is None:
1464            s3 += '%15s' % ('n/a')
1465        else:
1466            s3 += '%15.5f' % (esd)
1467
1468def ComputeDepESD(covMatrix,varyList,parmDict):
1469    '''Compute uncertainties for dependent parameters from independent ones
1470    returns a dictionary containing the esd values for dependent parameters
1471    '''
1472    sigmaDict = {}
1473    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
1474        #if invmultarr is None: continue # probably not needed
1475#        try:
1476#            valuelist = [parmDict[var] for var in mapvars]
1477#        except KeyError:
1478#            continue
1479        # get the v-covar matrix for independent parameters
1480        vcov = np.zeros((len(mapvars),len(mapvars)))
1481        for i1,name1 in enumerate(mapvars):
1482            if name1 not in varyList: continue
1483            iv1 = varyList.index(name1)
1484            for i2,name2 in enumerate(mapvars):
1485                if name2 not in varyList: continue
1486                iv2 = varyList.index(name2)
1487                vcov[i1][i2] = covMatrix[iv1][iv2]
1488        # vec is the vector that multiplies each of the independent values
1489        for v,vec in zip(varlist,invmultarr):
1490            sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec)))
1491    return sigmaDict
1492
1493def _FormatConstraint(RelDict,RelVal):
1494    '''Formats a Constraint or Function for use in a convenient way'''
1495    linelen = 45
1496    s = [""]
1497    for var,val in RelDict.items():
1498        if var.startswith('_'): continue
1499        if len(s[-1]) > linelen: s.append(' ')
1500        m = val
1501        if s[-1] != "" and m >= 0:
1502            s[-1] += ' + '
1503        elif s[-1] != "":
1504            s[-1] += ' - '
1505            m = abs(m)
1506        s[-1] += '%.3f*%s '%(m,var)
1507    if len(s[-1]) > linelen: s.append(' ')
1508    if RelVal is None:
1509        s[-1] += ' = New variable'
1510    else:
1511        s[-1] += ' = ' + RelVal
1512    s1 = ''
1513    for s2 in s:
1514        if s1 != '': s1 += '\n\t'
1515        s1 += s2
1516    return s1
1517
1518def _showEquiv(varlist,mapvars,invmultarr):
1519    '''Format an equivalence relationship
1520    note that
1521    varlist,           mapvars,     invmultarr
1522    are elements of
1523    dependentParmList, indParmList, invarrayList
1524    '''
1525    for i,mv in enumerate(mapvars):
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    return s1
1540
1541def VarRemapShow(varyList,inputOnly=False):
1542    '''List out the saved relationships. This should be done after the constraints have been
1543    defined using :func:`StoreEquivalence`, :func:`GroupConstraints` and :func:`GenerateConstraints`.
1544
1545    :returns: a string containing the details of the contraint relationships
1546    '''
1547    s = ''
1548    if len(fixedVarList) > 0:
1549        s += 'Fixed Parameters:\n'
1550        for v in fixedVarList:
1551            s += '    ' + v + '\n'
1552    if not inputOnly:
1553        s += 'User-supplied parameter mapping relations:\n'
1554    symout = ''
1555    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,symGenList
1556
1557    for varlist,mapvars,multarr,invmultarr,symFlag in zip(
1558        dependentParmList,indParmList,arrayList,invarrayList,symGenList):
1559        for i,mv in enumerate(mapvars):
1560            if multarr is None:
1561#                s1 = '  ' + str(mv) + ' is equivalent to parameter(s): '
1562                if len(varlist) == 1:
1563                    s1 = '   ' + str(mv) + ' is equivalent to '
1564                else:
1565                    s1 = '   ' + str(mv) + ' is equivalent to parameters: '
1566                j = 0
1567                for v,m in zip(varlist,invmultarr):
1568                    if debug: print ('v,m[0]: ',v,m[0])
1569                    if len(s1.split('\n')[-1]) > 75: s1 += '\n        '
1570                    if j > 0: s1 += ' & '
1571                    j += 1
1572                    s1 += str(v)
1573                    if m != 1:
1574                        s1 += " / " + str(m[0])
1575                if symFlag:
1576                    symout += s1 + '\n'
1577                else:
1578                    s += s1 + '\n'
1579                continue
1580            s += %s = ' % mv
1581            j = 0
1582            for m,v in zip(multarr[i,:],varlist):
1583                if m == 0: continue
1584                if j > 0: s += ' + '
1585                j += 1
1586                s += '(%s * %s)' % (m,v)
1587            if mv in varyList: s += ' VARY'
1588            s += '\n'
1589    if symout:
1590        s += 'Symmetry-generated relations:\n' + symout
1591    if inputOnly: return s
1592    s += 'Inverse parameter mapping relations:\n'
1593    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
1594        for i,mv in enumerate(varlist):
1595            s += %s = ' % mv
1596            j = 0
1597            for m,v in zip(invmultarr[i,:],mapvars):
1598                if m == 0: continue
1599                if j > 0: s += ' + '
1600                j += 1
1601                s += '(%s * %s)' % (m,v)
1602            s += '\n'
1603    return s
1604
1605def GetSymEquiv():
1606    '''Return the automatically generated (equivalence) relationships.
1607
1608    :returns: a list of strings containing the details of the contraint relationships
1609    '''
1610    symout = []
1611    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,symGenList
1612
1613    for varlist,mapvars,multarr,invmultarr,symFlag in zip(
1614        dependentParmList,indParmList,arrayList,invarrayList,symGenList):
1615        for i,mv in enumerate(mapvars):
1616            if not symFlag: continue
1617            if multarr is None:
1618                s1 = ''
1619                s2 = ' = ' + str(mv)
1620                j = 0
1621                if len(varlist) == 1:
1622                    # format the way Bob prefers
1623                    if invmultarr[0][0] == 1: 
1624                        s1 = str(varlist[0]) + ' = ' + str(mv)
1625                    else:
1626                        s1 = str(varlist[0]) + ' = ' + str(
1627                            invmultarr[0][0]) + ' * '+ str(mv)
1628                    symout.append(s1)
1629                    continue
1630                for v,m in zip(varlist,invmultarr):
1631                    if debug: print ('v,m[0]: ',v,m[0])
1632                    if len(s1.split('\n')[-1]) > 75: s1 += '\n        '
1633                    if j > 0: s1 += ' =  '
1634                    j += 1
1635                    s1 += str(v)
1636                    if m != 1:
1637                        s1 += " / " + str(m[0])
1638                symout.append(s1+s2)
1639                continue
1640            else:
1641                s = %s = ' % mv
1642                j = 0
1643                for m,v in zip(multarr[i,:],varlist):
1644                    if m == 0: continue
1645                    if j > 0: s += ' + '
1646                    j += 1
1647                    s += '(%s * %s)' % (m,v)
1648                print ('unexpected sym op='+s)
1649    return symout
1650
1651def Dict2Deriv(varyList,derivDict,dMdv):
1652    '''Compute derivatives for Independent Parameters from the
1653    derivatives for the original parameters
1654
1655    :param list varyList: a list of parameters names that will be varied
1656
1657    :param dict derivDict: a dict containing derivatives for parameter values keyed by the
1658      parameter names.
1659
1660    :param list dMdv: a Jacobian, as a list of np.array containing derivatives for dependent
1661      parameter computed from derivDict
1662
1663    '''
1664    global dependentParmList,arrayList,invarrayList,indParmList,invarrayList
1665    for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
1666        for i,name in enumerate(mapvars):
1667            # grouped parameters: need to add in the derv. w/r
1668            # dependent variables to the independent ones
1669            if name not in varyList: continue # skip if independent var not varied
1670            if multarr is None:
1671                for v,m in zip(varlist,invmultarr):
1672                    if debug: print ('start dMdv',dMdv[varyList.index(name)])
1673                    if debug: print ('add derv',v,'/',m[0],'to derv',name,'add=',derivDict[v] / m[0])
1674                    if m == 0: continue
1675                    dMdv[varyList.index(name)] += derivDict[v] / m[0]
1676            else:
1677                for v,m in zip(varlist,multarr[i,:]):
1678                    if debug: print ('start dMdv',dMdv[varyList.index(name)])
1679                    if debug: print ('add derv',v,'*',m,'to derv',name,'add=',m * derivDict[v])
1680                    if m == 0: continue
1681                    dMdv[varyList.index(name)] += m * derivDict[v]
1682
1683def Map2Dict(parmDict,varyList):
1684    '''Create (or update) the Independent Parameters from the original
1685    set of Parameters
1686
1687    Removes dependent variables from the varyList
1688
1689    This should be done once, after the constraints have been
1690    defined using :func:`StoreEquivalence`,
1691    :func:`GroupConstraints` and :func:`GenerateConstraints` and
1692    before any parameter refinement is done. This completes the parameter
1693    dictionary by defining independent parameters and it satisfies the
1694    constraint equations in the initial parameters
1695
1696    :param dict parmDict: a dict containing parameter values keyed by the
1697      parameter names.
1698      This will contain updated values for both dependent and independent
1699      parameters after Dict2Map is called. It will also contain some
1700      unexpected entries of every constant value {'0':0.0} & {'1.0':1.0},
1701      which do not cause any problems.
1702
1703    :param list varyList: a list of parameters names that will be varied
1704   
1705
1706    '''
1707    # process the Independent vars: remove dependent ones from varylist
1708    # and then compute values for the Independent ones from their dependents
1709    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
1710    for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
1711        for item in varlist:
1712            try:
1713                varyList.remove(item)
1714            except ValueError:
1715                pass
1716        if multarr is None: continue
1717        valuelist = [parmDict[var] for var in varlist]
1718        parmDict.update(zip(mapvars,
1719                            np.dot(multarr,np.array(valuelist)))
1720                        )
1721    # now remove fixed parameters from the varyList
1722    global fixedVarList
1723    for item in fixedVarList:
1724        try:
1725            varyList.remove(item)
1726        except ValueError:
1727            pass
1728    # Set constrained parameters to their fixed values
1729    parmDict.update(fixedDict)
1730
1731def Dict2Map(parmDict,varyList):
1732    '''Applies the constraints defined using :func:`StoreEquivalence`,
1733    :func:`GroupConstraints` and :func:`GenerateConstraints` by changing
1734    values in a dict containing the parameters. This should be
1735    done before the parameters are used for any computations
1736
1737    :param dict parmDict: a dict containing parameter values keyed by the
1738      parameter names.
1739      This will contain updated values for both dependent and independent
1740      parameters after Dict2Map is called. It will also contain some
1741      unexpected entries of every constant value {'0':0.0} & {'1.0':1.0},
1742      which do not cause any problems.
1743
1744    :param list varyList: a list of parameters names that will be varied
1745   
1746    '''
1747    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
1748    # reset fixed values (should not be needed, but very quick)
1749    # - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended
1750    # not needed, but also not a problem - BHT
1751    parmDict.update(fixedDict)
1752    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
1753        #if invmultarr is None: continue
1754        try: 
1755            valuelist = [parmDict[var] for var in mapvars]
1756        except KeyError:
1757            continue
1758        parmDict.update(zip(varlist,np.dot(invmultarr,np.array(valuelist))))
1759
1760#======================================================================
1761# internal routines follow (these routines are unlikely to be called
1762# from outside the module)
1763
1764def GramSchmidtOrtho(a,nkeep=0):
1765    '''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to
1766    find orthonormal unit vectors relative to first row.
1767
1768    If nkeep is non-zero, the first nkeep rows in the array are not changed
1769   
1770    input:
1771       arrayin: a 2-D non-singular square array
1772    returns:
1773       a orthonormal set of unit vectors as a square array
1774    '''
1775    def proj(a,b):
1776        'Projection operator'
1777        return a*(np.dot(a,b)/np.dot(a,a))
1778    for j in range(nkeep,len(a)):
1779        for i in range(j):
1780            a[j] -= proj(a[i],a[j])
1781        if np.allclose(np.linalg.norm(a[j]),0.0):
1782            raise ConstraintException("Singular input to GramSchmidtOrtho")
1783        a[j] /= np.linalg.norm(a[j])
1784    return a
1785
1786def _FillArray(sel,dict,collist,FillDiagonals=False):
1787    '''Construct a n by n matrix (n = len(collist)
1788    filling in the rows using the relationships defined
1789    in the dictionaries selected by sel
1790
1791    If FillDiagonals is True, diagonal elements in the
1792    array are set to 1.0
1793    '''
1794    n = len(collist)
1795    if FillDiagonals:
1796        arr = np.eye(n)
1797    else:
1798        arr = np.zeros(2*[n,])
1799    # fill the top rows
1800    for i,cnum in enumerate(sel):
1801        for j,var in enumerate(collist):
1802            arr[i,j] = dict[cnum].get(var,0)
1803    return arr
1804
1805def _SwapColumns(i,m,v):
1806    '''Swap columns in matrix m as well as the labels in v
1807    so that element (i,i) is replaced by the first non-zero element in row i after that element
1808
1809    Throws an exception if there are no non-zero elements in that row
1810    '''
1811    for j in range(i+1,len(v)):
1812        if not np.allclose(m[i,j],0):
1813            m[:,(i,j)] = m[:,(j,i)]
1814            v[i],v[j] = v[j],v[i]
1815            return
1816    else:
1817        raise ConstraintException('Singular input')
1818
1819def _RowEchelon(m,arr,collist):
1820    '''Convert the first m rows in Matrix arr to row-echelon form
1821    exchanging columns in the matrix and collist as needed.
1822
1823    throws an exception if the matrix is singular because
1824    the first m rows are not linearly independent
1825    '''
1826    for i in range(m):
1827        if np.allclose(arr[i,i],0):
1828            _SwapColumns(i,arr,collist)
1829        arr[i,:] /= arr[i,i] # normalize row
1830        # subtract current row from subsequent rows to set values to left of diagonal to 0
1831        for j in range(i+1,m):
1832            arr[j,:] -= arr[i,:] * arr[j,i]
1833
1834if __name__ == "__main__":
1835    parmdict = {}
1836    constrDict = [
1837        {'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},
1838        {'0:0:eA': 0.0},
1839        {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
1840        {'1::C(10,0,1)': 1.0, '2::C(10,0,1)': 1.0},
1841        {'1::AUiso:0': 1.0, '0::AUiso:0': 1.0},
1842        {'0::A0': 0.0}
1843        ]
1844    fixedList = ['5.0', '0', None, None, '1.0', '0']
1845    StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), ))
1846    #StoreEquivalence('1::atomx:3',('2::atomx:3', ('2::atomz:3',2,), )) # error: dependent & independent vars mixed
1847    #StoreEquivalence('1::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) # error: dependent vars repeated
1848    #StoreEquivalence('0:1:eA',('0:0:eA',)) # error: equiv & fixed
1849    #StoreEquivalence('0:99:Scale',('0:12:Scale',)) # error: equiv & constrained
1850    #StoreEquivalence('0:12:Scale',('0:99:Scale',)) # error: equiv & constrained
1851    varylist = ['2::atomx:3',
1852                '2::C(10,6,1)', '1::C(10,6,1)',
1853                '2::atomy:3', '2::atomz:3',
1854                '0:12:Scale', '0:11:Scale', '0:14:Scale', '0:13:Scale', '0:0:Scale']
1855#    e,w = CheckConstraints([,
1856#                     [{'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}],
1857#                     ['1.0'])
1858#    if e: print 'error=',e
1859#    if w: print 'error=',w
1860#    varyList = ['0::A0', '0::AUiso:0', '0::Afrac:1', '0::Afrac:2', '0::Afrac:3', '0::Afrac:4',
1861#       '0::dAx:5', '0::dAy:5', '0::dAz:5', '0::AUiso:5', ':0:Back;0', ':0:Back;1', ':0:Back;2', ':0:Back;3',
1862#       ':0:Back;4', ':0:Back;5', ':0:Back;6', ':0:Back;7', ':0:Back;8', ':0:Back;9', ':0:Back;10', ':0:Back;11'
1863#       :0:U', ':0:V', ':0:W', ':0:X', ':0:Y', ':0:Scale', ':0:DisplaceX', ':0:DisplaceY']
1864#    constrDict = [
1865#        {'0::Afrac:4': 24.0, '0::Afrac:1': 16.0, '0::Afrac:3': 24.0, '0::Afrac:2': 16.0},
1866#        {'0::Afrac:1': 1.0, '0::Afrac:2': 1.0},
1867#        {'0::Afrac:4': 1.0, '0::Afrac:3': 1.0}]
1868#    fixedList = ['40.0', '1.0', '1.0']
1869
1870    errmsg, warnmsg = CheckConstraints(varylist,constrDict,fixedList)
1871    if errmsg:
1872        print ("*** Error ********************")
1873        print (errmsg)
1874    if warnmsg:
1875        print ("*** Warning ********************")
1876        print (warnmsg)
1877    if errmsg or warnmsg:
1878        sys.exit()
1879    groups,parmlist = GroupConstraints(constrDict)
1880    GenerateConstraints(groups,parmlist,varylist,constrDict,fixedList)
1881    print (VarRemapShow(varylist))
1882    parmdict.update( {
1883        '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,
1884        '0:0:eA': 0.0,
1885        '2::C(10,6,1)': 0.2, '1::C(10,6,1)': 0.3,
1886        '1::C(10,0,1)': 0.2, '2::C(10,0,1)': 0.3,
1887        '1::AUiso:0': 0.02, '0::AUiso:0': 0.03,
1888        '0::A0': 0.0,
1889        '2::atomx:3':0.23,'2::atomy:3':-.23, '2::atomz:3':-0.11,
1890        })
1891    print ('parmdict start',parmdict)
1892    print ('varylist start',varylist)
1893    before = parmdict.copy()
1894    Map2Dict(parmdict,varylist)
1895    print ('parmdict before and after Map2Dict')
1896    print ('  key / before / after')
1897    for key in sorted(list(parmdict.keys())):
1898        print ('  '+key,'\t',before.get(key),'\t',parmdict[key])
1899    print ('varylist after',varylist)
1900    before = parmdict.copy()
1901    Dict2Map(parmdict,varylist)
1902    print ('after Dict2Map')
1903    print ('  key / before / after')
1904    for key in sorted(list(parmdict.keys())):
1905        print ('  '+key,'\t',before.get(key),'\t',parmdict[key])
1906#    dMdv = len(varylist)*[0]
1907#    deriv = {}
1908#    for i,v in enumerate(parmdict.keys()): deriv[v]=i
1909#    Dict2Deriv(varylist,deriv,dMdv)
Note: See TracBrowser for help on using the repository browser.