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