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