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