source: trunk/GSASIImapvars.py @ 3136

Last change on this file since 3136 was 3136, checked in by vondreele, 4 years ago

make GSAS-II python 3.6 compliant & preserve python 2.7 use;changes:
do from future import division, print_function for all GSAS-II py sources
all menu items revised to be py 2.7/3.6 compliant
all wx.OPEN --> wx.FD_OPEN in file dialogs
all integer divides (typically for image pixel math) made explicit with ; ambiguous ones made floats as appropriate
all print "stuff" --> print (stuff)
all print >> pFile,'stuff' --> pFile.writeCIFtemplate('stuff')
all read file opens made explicit 'r' or 'rb'
all cPickle imports made for py2.7 or 3.6 as cPickle or _pickle; test for '2' platform.version_tuple[0] for py 2.7
define cPickleload to select load(fp) or load(fp,encoding='latin-1') for loading gpx files; provides cross compatibility between py 2.7/3.6 gpx files
make dict.keys() as explicit list(dict.keys()) as needed (NB: possible source of remaining py3.6 bugs)
make zip(a,b) as explicit list(zip(a,b)) as needed (NB: possible source of remaining py3.6 bugs)
select unichr/chr according test for '2' platform.version_tuple[0] for py 2.7 (G2pwdGUI * G2plot) for special characters
select wg.EVT_GRID_CELL_CHANGE (classic) or wg.EVT_GRID_CELL_CHANGED (phoenix) in grid Bind
maxint --> maxsize; used in random number stuff
raise Exception,"stuff" --> raise Exception("stuff")
wx 'classic' sizer.DeleteWindows?() or 'phoenix' sizer.Clear(True)
wx 'classic' SetToolTipString?(text) or 'phoenix' SetToolTip?(wx.ToolTip?(text)); define SetToolTipString?(self,text) to handle the choice in plots
status.SetFields? --> status.SetStatusText?
'classic' AddSimpleTool? or 'phoenix' self.AddTool? for plot toolbar; Bind different as well
define GetItemPydata? as it doesn't exist in wx 'phoenix'
allow python versions 2.7 & 3.6 to run GSAS-II
Bind override commented out - no logging capability (NB: remove all logging code?)
all import ContentsValidator? open filename & test if valid then close; filepointer removed from Reader
binary importers (mostly images) test for 'byte' type & convert as needed to satisfy py 3.6 str/byte rules

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 53.2 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2017-10-23 16:39:16 +0000 (Mon, 23 Oct 2017) $
4# $Author: vondreele $
5# $Revision: 3136 $
6# $URL: trunk/GSASIImapvars.py $
7# $Id: GSASIImapvars.py 3136 2017-10-23 16:39:16Z vondreele $
8########### SVN repository information ###################
9"""
10*GSASIImapvars: Parameter constraints*
11======================================
12
13Module to implements algebraic contraints, parameter redefinition
14and parameter simplification contraints.
15
16Parameter redefinition (new vars) is done by creating one or more relationships
17between a set of parameters
18
19::
20
21   Mx1 * Px + My1 * Py +...
22   Mx2 * Px + Mz2 * Pz + ...
23
24where Pj is a parameter name and Mjk is a constant.
25
26Constant constraint Relations can also be supplied in the form of an equation:
27
28::
29
30  nx1 * Px + ny1 * Py +... = C1
31
32where Cn is a constant. These equations define an algebraic
33constrant.
34
35Parameters can also be "fixed" (held), which prevents them from being refined.
36
37All of the above three cases are input using routines
38GroupConstraints and GenerateConstraints. The input consists of a list of
39relationship dictionaries:
40
41.. code-block:: python
42
43    constrDict = [
44         {'0:12:Scale': 2.0, '0:14:Scale': 4.0, '0:13:Scale': 3.0, '0:0:Scale': 0.5},
45         {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
46         {'0::A0': 0.0}]
47    fixedList = ['5.0', None, '0']
48
49Where the dictionary defines the first part of an expression and the corresponding fixedList
50item is either None (for parameter redefinition) or the constant value, for a constant
51constraint equation. A dictionary that contains a single term defines a variable that
52will be fixed (held). The multiplier and the fixedList value in this case are ignored.
53
54Parameters can also be equivalenced or "slaved" to another parameter, such that one
55(independent) parameter is equated to several (now dependent) parameters. In
56algebraic form this is:
57
58::
59
60   P0 = M1 * P1 = M2 * P2 = ...
61
62Thus parameters P0, P1 and P2,... are linearly equivalent. Routine StoreEquivalence is
63used to specify these equivalences.
64
65Parameter redefinition (new vars) describes a new, independent, parameter, which
66is defined in terms of dependent parameters that are defined in the
67Model, while fixed constrained relations effectively reduce the complexity
68of the Model by removing a degree of freedom. It is possible for a parameter to appear
69in both a parameter redefinition expression and a fixed constraint equation, but a
70parameter cannot be used a parameter equivalance cannot be used elsewhere (not fixed,
71constrained or redefined). Likewise a fixed parameter cannot be used elsewhere (not
72equivalanced, constrained or redefined).
73
74Relationships are grouped so that a set of dependent parameters appear
75in only one group (done in routine GroupConstraints.) Note that if a
76group contains relationships/equations that involve N dependent
77parameters, there must exist N-C new parameters, where C is the number
78of contraint equations in the group. Routine GenerateConstraints takes
79the output from GroupConstraints and generates the
80"missing" relationships and saves that information in the module's
81global variables. Each generated parameter is named sequentially using paramPrefix.
82
83A list of parameters that will be varied is specified as input to GenerateConstraints
84(varyList). A fixed parameter will simply be removed from this list preventing that
85parameter from being varied. Note that all parameters in a constraint relationship
86must specified as varied (appear in varyList) or none can be varied. This is
87checked in GenerateConstraints. Likewise, if all parameters in a constraint are
88not referenced in a refinement, the constraint is ignored, but if some parameters
89in a constraint group are not referenced in a refinement, but others are this
90constitutes and error.
91
92* When a new variable is created, the variable is assigned the name associated
93  in the constraint definition or it is assigned a default name of form
94  ``::constr<n>`` (see paramPrefix). The vary setting for variables used in the
95  constraint are ignored.
96  Note that any generated "missing" relations are not varied. Only
97  the input relations can be are varied.
98 
99* If all parameters in a fixed constraint equation are varied, the generated "missing"
100  relations in the group are all varied. This provides the N-C degrees of freedom.
101
102*External Routines*
103-------------------
104
105To define a set of constrained and unconstrained relations, one
106defines a list of dictionary defining constraint parameters and their
107values, a list of fixed values for each constraint and a list of
108parameters to be varied. In addition, one uses
109:func:`StoreEquivalence` to define parameters that are equivalent. One
110can then use :func:`CheckConstraints` to check that the input is
111internally consistent and finally :func:`GroupConstraints` and
112:func:`GenerateConstraints` to generate the internally used
113tables. Routines :func:`Map2Dict` is used to initialize the parameter
114dictionary and :func:`Dict2Map`, :func:`Dict2Deriv`, and
115:func:`ComputeDepESD` are used to apply constraints. Routine
116:func:`VarRemapShow` is used to print out the constraint information,
117as stored by :func:`GenerateConstraints`.
118
119:func:`InitVars`
120  This is optionally used to clear out all defined previously defined constraint information
121 
122:func:`StoreEquivalence`
123  To implement parameter redefinition, one calls StoreEquivalence. This should be called for every set of
124  equivalence relationships. There is no harm in using StoreEquivalence with the same independent variable:
125
126  .. code-block:: python
127
128       StoreEquivalence('x',('y',))
129       StoreEquivalence('x',('z',))
130
131  or equivalently
132
133  .. code-block:: python
134
135       StoreEquivalence('x',('y','z'))
136
137  The latter will run more efficiently. Note that mixing independent and dependent variables would
138  require assignments to be done in a particular order and thus is not is not allowed. This will cause an
139  error:
140
141  .. code-block:: python
142
143        StoreEquivalence('x',('y',))
144        StoreEquivalence('y',('z',))
145       
146  Use StoreEquivalence before calling GenerateConstraints or CheckConstraints
147
148:func:`CheckConstraints`
149   To check that input in internally consistent, use CheckConstraints
150
151:func:`Map2Dict`
152   To determine values for the parameters created in this module, one
153   calls Map2Dict. This will not apply contraints.
154
155:func:`Dict2Map`
156   To take values from the new independent parameters and constraints,
157   one calls Dict2Map. This will apply contraints.
158
159:func:`Dict2Deriv`
160   Use Dict2Deriv to determine derivatives on independent parameters
161   from those on dependent ones
162
163:func:`ComputeDepESD`     
164   Use ComputeDepESD to compute uncertainties on dependent variables
165
166:func:`VarRemapShow`
167   To show a summary of the parameter remapping, one calls VarRemapShow
168
169*Global Variables*
170------------------
171
172dependentParmList:
173   contains a list by group of lists of
174   parameters used in the group. Note that parameters listed in
175   dependentParmList should not be refined as they will not affect
176   the model
177
178indParmList:
179     a list of groups of Independent parameters defined in
180     each group. This contains both parameters used in parameter
181     redefinitions as well as names of generated new parameters.
182
183fixedVarList:
184     a list of variables that have been 'fixed'
185     by defining them as equal to a constant (::var: = 0). Note that
186     the constant value is ignored at present. These variables are
187     later removed from varyList which prevents them from being refined.
188     Unlikely to be used externally.
189
190arrayList:
191     a list by group of relationship matrices to relate
192     parameters in dependentParmList to those in indParmList. Unlikely
193     to be used externally.
194
195invarrayList:
196     a list by group of relationship matrices to relate
197     parameters in indParmList to those in dependentParmList. Unlikely
198     to be used externally.
199
200fixedDict:
201     a dictionary containing the fixed values corresponding
202     to parameter equations.  The dict key is an ascii string, but the
203     dict value is a float.  Unlikely to be used externally.
204
205symGenList:
206     a list of boolean values that will be True to indicate that a constraint
207     (only equivalences) is generated by symmetry (or Pawley overlap)
208     
209problemVars:
210     a list containing variables that show up in constraints producing errors
211
212*Routines*
213----------
214
215Note that parameter names in GSAS-II are strings of form ``<ph>:<hst>:<nam>``
216
217"""
218
219from __future__ import division, print_function
220import numpy as np
221import sys
222import GSASIIpath
223GSASIIpath.SetVersionNumber("$Revision: 3136 $")
224# data used for constraints;
225debug = False # turns on printing as constraint input is processed
226# note that constraints are stored listed by contraint groups, where each constraint
227# group contains those parameters that must be handled together
228dependentParmList = [] # contains a list of parameters in each group
229# note that parameters listed in dependentParmList should not be refined
230arrayList = [] # a list of of relationship matrices
231invarrayList = [] # a list of inverse relationship matrices
232indParmList = [] # a list of names for the new parameters
233fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations
234               # key is original ascii string, value is float
235fixedVarList = [] # List of variables that should not be refined
236symGenList = [] # Flag if constraint is generated by symmetry
237problemVars = [] # variables causing errors
238dependentVars = []
239independentVars = []
240
241# prefix for parameter names
242paramPrefix = "::constr"
243consNum = 0 # number of the next constraint to be created
244
245def InitVars():
246    '''Initializes all constraint information'''
247    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,consNum,symGenList
248    dependentParmList = [] # contains a list of parameters in each group
249    arrayList = [] # a list of of relationship matrices
250    invarrayList = [] # a list of inverse relationship matrices
251    indParmList = [] # a list of names for the new parameters
252    fixedDict = {} # a dictionary containing the fixed values corresponding to defined parameter equations
253    symGenList = [] # Flag if constraint is generated by symmetry
254    consNum = 0 # number of the next constraint to be created
255
256def VarKeys(constr):
257    """Finds the keys in a constraint that represent variables
258    e.g. eliminates any that start with '_'
259
260    :param dict constr: a single constraint entry of form::
261
262        {'var1': mult1, 'var2': mult2,... '_notVar': val,...}
263
264        (see :func:`GroupConstraints`)
265    :returns: a list of keys where any keys beginning with '_' are
266      removed.
267    """
268    return [i for i in constr.keys() if not i.startswith('_')]
269
270
271def GroupConstraints(constrDict):
272    """divide the constraints into groups that share no parameters.
273
274    :param dict constrDict: a list of dicts defining relationships/constraints
275
276    ::
277   
278       constrDict = [{<constr1>}, {<constr2>}, ...]
279
280    where {<constr1>} is {'var1': mult1, 'var2': mult2,... }
281
282    :returns: two lists of lists:
283   
284      * a list of grouped contraints where each constraint grouped containts a list
285        of indices for constraint constrDict entries
286      * a list containing lists of parameter names contained in each group
287     
288      """
289    assignedlist = [] # relationships that have been used
290    groups = [] # contains a list of grouplists
291    ParmList = []
292    for i,consi in enumerate(constrDict):
293        if i in assignedlist: continue # already in a group, skip
294        # starting a new group
295        grouplist = [i,]
296        assignedlist.append(i)
297        groupset = set(VarKeys(consi))
298        changes = True # always loop at least once
299        while(changes): # loop until we can't find anything to add to the current group
300            changes = False # but don't loop again unless we find something
301            for j,consj in enumerate(constrDict):
302                if j in assignedlist: continue # already in a group, skip
303                if len(set(VarKeys(consj)) & groupset) > 0: # true if this needs to be added
304                    changes = True
305                    grouplist.append(j)
306                    assignedlist.append(j)
307                    groupset = groupset | set(VarKeys(consj))
308        group = sorted(grouplist)
309        varlist = sorted(list(groupset))
310        groups.append(group)
311        ParmList.append(varlist)
312    return groups,ParmList
313
314def CheckConstraints(varyList,constrDict,fixedList):
315    '''Takes a list of relationship entries comprising a group of
316    constraints and checks for inconsistencies such as conflicts in
317    parameter/variable definitions and or inconsistently varied parameters.
318
319    :param list varyList: a list of parameters names that will be varied
320
321    :param dict constrDict: a list of dicts defining relationships/constraints
322      (as created in :func:`GSASIIstrIO.ProcessConstraints` and
323      documented in :func:`GroupConstraints`)
324
325    :param list fixedList: a list of values specifying a fixed value for each
326      dict in constrDict. Values are either strings that can be converted to
327      floats or ``None`` if the constraint defines a new parameter rather
328      than a constant.
329
330    :returns: two strings:
331
332      * the first lists conflicts internal to the specified constraints
333      * the second lists conflicts where the varyList specifies some
334        parameters in a constraint, but not all
335       
336      If there are no errors, both strings will be empty
337    '''
338    import re
339    global dependentParmList,arrayList,invarrayList,indParmList,consNum
340    errmsg = ''
341    warnmsg = ''
342    global problemVars
343    problemVars = []
344    fixVlist = []
345    # process fixed variables (holds)
346    for cdict in constrDict:
347        # N.B. No "_" names in holds
348        if len(cdict) == 1:
349            fixVlist.append(list(cdict.keys())[0])
350   
351    # process equivalences: make a list of dependent and independent vars
352    #    and check for repeated uses (repetition of a parameter as an
353    #    independent var is OK)
354    indepVarList = []
355    depVarList = []
356    multdepVarList = []
357    for varlist,mapvars,multarr,invmultarr in zip(
358        dependentParmList,indParmList,arrayList,invarrayList):
359        if multarr is None: # an equivalence
360            zeromult = False
361            for mv in mapvars:
362                varied = 0
363                notvaried = ''
364                if mv in varyList:
365                    varied += 1
366                else:
367                    if notvaried: notvaried += ', '
368                    notvaried += mv
369                if mv not in indepVarList: indepVarList.append(mv)
370                for v,m in zip(varlist,invmultarr):
371                    if v in indepVarList:
372                        errmsg += '\nVariable '+v+' is used to set values in a constraint before its value is set in another constraint\n'
373                        if v not in problemVars: problemVars.append(v)
374                    if m == 0: zeromult = True
375                    if v in varyList:
376                        varied += 1
377                    else:
378                        if notvaried: notvaried += ', '
379                        notvaried += v
380                    if v in depVarList:
381                        multdepVarList.append(v)
382                    else:
383                        depVarList.append(v)
384            if varied > 0 and varied != len(varlist)+1:
385                warnmsg += "\nNot all variables refined in equivalence:\n\t"
386                s = ""
387                for v in varlist:
388                    if s != "": s+= " & "
389                    s += str(v)           
390                warnmsg += str(mv) + " => " + s
391                warnmsg += '\nNot refined: ' + notvaried + '\n'
392            if zeromult:
393                errmsg += "\nZero multiplier is invalid in equivalence:\n\t"
394                s = ""
395                for v in varlist:
396                    if s != "": s+= " & "
397                    s += str(v)           
398                errmsg += str(mv) + " => " + s + '\n'
399
400    # check for errors:
401    if len(multdepVarList) > 0:
402        errmsg += "\nThe following parameters(s) are used in conflicting Equivalence relations as dependent variables:\n"
403        s = ''
404        for var in sorted(set(multdepVarList)):
405            if v not in problemVars: problemVars.append(v)
406            if s != "": s+= ", "
407            s += str(var)           
408        errmsg += '\t'+ s + '\n'
409    equivVarList = list(set(indepVarList).union(set(depVarList)))
410    if debug: print ('equivVarList',equivVarList)
411    inboth = set(fixVlist).intersection(set(equivVarList))
412    if len(inboth) > 0:
413        errmsg += "\nThe following parameter(s) are used in both Equivalence and Fixed constraints:\n"
414        s = ''
415        for var in sorted(inboth):
416            if var not in problemVars: problemVars.append(var)
417            if s != "": s+= ", "
418            s += str(var)
419        errmsg += '\t'+ s + '\n'
420
421    groups,parmlist = GroupConstraints(constrDict)
422    # scan through parameters in each relationship. Are all varied? If only some are
423    # varied, create a warning message.
424    for group,varlist in zip(groups,parmlist):
425        if len(varlist) == 1: continue
426        for rel in group:
427            varied = 0
428            notvaried = ''
429            for var in constrDict[rel]:
430                if var.startswith('_'): continue
431                if not re.match('[0-9]*:[0-9\*]*:',var):
432                    warnmsg += "\nVariable "+str(var)+" does not begin with a ':'"
433                if var in varyList:
434                    varied += 1
435                else:
436                    if notvaried: notvaried += ', '
437                    notvaried += var
438                if var in fixVlist:
439                    errmsg += '\nParameter '+var+" is Fixed and used in a constraint:\n\t"
440                    if var not in problemVars: problemVars.append(var)
441                    errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
442            if varied > 0 and varied != len(VarKeys(constrDict[rel])):
443                warnmsg += "\nNot all variables refined in constraint:\n\t"
444                warnmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
445                warnmsg += '\nNot refined: ' + notvaried + '\n'
446    if errmsg or warnmsg:
447        return errmsg,warnmsg
448
449    # now look for process each group and create the relations that are needed to form
450    # non-singular square matrix
451    for group,varlist in zip(groups,parmlist):
452        if len(varlist) == 1: continue # a constraint group with a single variable can be ignored
453        if len(varlist) < len(group): # too many relationships -- no can do
454            errmsg += "\nOver-constrained input. "
455            errmsg += "There are more constraints " + str(len(group))
456            errmsg += "\n\tthan variables " + str(len(varlist)) + "\n"
457            for rel in group:
458                errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
459                errmsg += "\n"
460                continue
461        try:
462            multarr = _FillArray(group,constrDict,varlist)
463            _RowEchelon(len(group),multarr,varlist)
464        except:
465            errmsg += "\nSingular input. "
466            errmsg += "There are internal inconsistencies in these constraints\n"
467            for rel in group:
468                errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
469                errmsg += "\n"
470            continue
471        try:
472            multarr = _FillArray(group,constrDict,varlist,FillDiagonals=True)
473            GramSchmidtOrtho(multarr,len(group))
474        except:
475            errmsg += "\nUnexpected singularity with constraints (in Gram-Schmidt)\n"
476            for rel in group:
477                errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
478                errmsg += "\n"
479            continue
480        mapvar = []
481        group = group[:]
482        # scan through all generated and input variables
483        # Check again for inconsistent variable use
484        # for new variables -- where varied and unvaried parameters get grouped
485        # together. I don't think this can happen when not flagged before, but
486        # it does not hurt to check again.
487        for i in range(len(varlist)):
488            varied = 0
489            notvaried = ''
490            if len(group) > 0:
491                rel = group.pop(0)
492                fixedval = fixedList[rel]
493                for var in VarKeys(constrDict[rel]):
494                    if var in varyList:
495                        varied += 1
496                    else:
497                        if notvaried: notvaried += ', '
498                        notvaried += var
499            else:
500                fixedval = None
501            if fixedval is None:
502                varname = paramPrefix + str(consNum) # assign a name to a variable
503                mapvar.append(varname)
504                consNum += 1
505            else:
506                mapvar.append(fixedval)
507            if varied > 0 and notvaried != '':
508                warnmsg += "\nNot all variables refined in generated constraint"
509                warnmsg += '\nPlease report this unexpected error\n'
510                for rel in group:
511                    warnmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
512                    warnmsg += "\n"
513                warnmsg += '\n\tNot refined: ' + notvaried + '\n'
514        try:
515            np.linalg.inv(multarr)           
516        except:
517            errmsg += "\nSingular input. "
518            errmsg += "The following constraints are not "
519            errmsg += "linearly independent\n\tor do not "
520            errmsg += "allow for generation of a non-singular set\n"
521            errmsg += 'Please report this unexpected error\n'
522            for rel in group:
523                errmsg += _FormatConstraint(constrDict[rel],fixedList[rel])
524                errmsg += "\n"
525    _setVarLists([])
526    return errmsg,warnmsg
527
528def GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict=None,SeqHist=None):
529    '''Takes a list of relationship entries comprising a group of
530    constraints and builds the relationship lists and their inverse
531    and stores them in global variables Also checks for internal
532    conflicts or inconsistencies in parameter/variable definitions.
533
534    :param list groups: a list of grouped contraints where each constraint
535      grouped containts a list of indices for constraint constrDict entries,
536      created in :func:`GroupConstraints` (returned as 1st value)
537
538    :param list parmlist: a list containing lists of parameter names
539      contained in each group, created in :func:`GroupConstraints`
540      (returned as 2nd value)
541
542    :param list varyList: a list of parameters names (strings of form
543      ``<ph>:<hst>:<nam>``) that will be varied. Note that this is changed here.
544   
545    :param dict constrDict: a list of dicts defining relationships/constraints
546      (as defined in :func:`GroupConstraints`)
547
548    :param list fixedList: a list of values specifying a fixed value for each
549      dict in constrDict. Values are either strings that can be converted to
550      floats, float values or None if the constraint defines a new parameter.
551     
552    :param dict parmDict: a dict containing all parameters defined in current
553      refinement.
554
555    :param int SeqHist: number of current histogram, when used in a sequential
556      refinement. None (default) otherwise. Wildcard variable names are
557      set to the current histogram, when found if not None.
558    '''
559    global dependentParmList,arrayList,invarrayList,indParmList,consNum
560    msg = ''
561    shortmsg = ''
562
563    # process fixed (held) variables
564    for cdict in constrDict:
565        if len(cdict) == 1:
566            fixedVarList.append(list(cdict.keys())[0])
567   
568    # process equivalences: make a list of dependent and independent vars
569    #    and check for repeated uses (repetition of a parameter as an
570    #    independent var is OK [A=B; A=C], but chaining: [A=B; B=C] is not good)
571    dropVarList = []
572    translateTable = {} # lookup table for wildcard referenced variables
573    for varlist,mapvars,multarr,invmultarr in zip(       # process equivalences
574        dependentParmList,indParmList,arrayList,invarrayList):
575        if multarr is None: # true only if an equivalence
576            zeromult = False
577            for mv in mapvars:
578                #s = ''
579                varied = 0
580                notvaried = ''
581                if mv in varyList:
582                    varied += 1
583                else:
584                    if notvaried: notvaried += ', '
585                    notvaried += mv
586                if parmDict is not None and mv not in parmDict:
587                    print ("Dropping equivalence for variable "+str(mv)+". Not defined in this refinement")
588                    if mv not in dropVarList: dropVarList.append(mv)
589                    #msg += "\nCannot equivalence to variable "+str(mv)+". Not defined in this refinement"
590                    #continue
591            for v,m in zip(varlist,invmultarr):
592                if parmDict is not None and v not in parmDict:
593                    print ("Dropping equivalence for dep. variable "+str(v)+". Not defined in this refinement")
594                    if v not in dropVarList: dropVarList.append(v)
595                    continue
596                if m == 0: zeromult = True
597                if v in varyList:
598                    varied += 1
599                else:
600                    if notvaried: notvaried += ', '
601                    notvaried += v
602            if varied > 0 and varied != len(varlist)+1:
603                msg += "\nNot all variables refined in equivalence:\n\t"
604                s = ""
605                for v in varlist:
606                    if s != "": s+= " & "
607                    s += str(v)           
608                msg += str(mv) + " => " + s
609                msg += '\nNot refined: ' + notvaried + '\n'
610            if zeromult:
611                msg += "\nZero multiplier is invalid in equivalence:\n\t"
612                s = ""
613                for v in varlist:
614                    if s != "": s+= " & "
615                    s += str(v)           
616                msg += str(mv) + " => " + s + '\n'
617
618    # scan through parameters in each relationship. Are all varied? If only some are
619    # varied, create an error message.
620    for group,varlist in zip(groups,parmlist):
621        if len(varlist) == 1: continue
622        for rel in group:
623            varied = 0
624            notvaried = ''
625            unused = 0
626            notused = ''
627            for var in constrDict[rel]:
628                if var.startswith('_'): continue
629                if var.split(':')[1] == '*' and SeqHist is not None:
630                    # convert wildcard var to reference current histogram; save translation in table
631                    sv = var.split(':')
632                    sv[1] = str(SeqHist)
633                    translateTable[var] = ':'.join(sv)
634                    var = translateTable[var]
635                if parmDict is not None and var not in parmDict:
636                    unused += 1
637                    if notvaried: notused += ', '
638                    notused += var
639                if var in varyList:
640                    varied += 1
641                else:
642                    if notvaried: notvaried += ', '
643                    notvaried += var
644                if var in fixedVarList:
645                    msg += '\nError: parameter '+var+" is Fixed and used in a constraint:\n\t"
646                    msg += _FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
647            #if unused > 0:# and unused != len(VarKeys(constrDict[rel])):
648            if unused > 0 and unused != len(VarKeys(constrDict[rel])):
649                #msg += "\nSome (but not all) variables in constraint are not defined:\n\t"
650                #msg += _FormatConstraint(constrDict[rel],fixedList[rel])
651                #msg += '\nNot used: ' + notused + '\n'
652                shortmsg += notused+" not used in constraint "+_FormatConstraint(constrDict[rel],fixedList[rel])
653            elif varied > 0 and varied != len(VarKeys(constrDict[rel])):
654                #msg += "\nNot all variables refined in constraint:\n\t"
655                #msg += _FormatConstraint(constrDict[rel],fixedList[rel])
656                #msg += '\nNot refined: ' + notvaried + '\n'
657                shortmsg += notvaried+" not varied in constraint "+_FormatConstraint(constrDict[rel],fixedList[rel])
658    # if there were errors found, go no farther
659    if shortmsg and SeqHist is not None:
660        if msg:
661            print (' *** ERROR in constraint definitions! ***')
662            print (msg)
663            raise Exception
664        print ('*** Sequential refinement: ignoring constraint definition(s): ***')
665        print (shortmsg)
666        msg = ''
667    elif shortmsg:
668        msg += shortmsg
669    if msg:
670        print (' *** ERROR in constraint definitions! ***')
671        print (msg)
672        raise Exception
673               
674    # now process each group and create the relations that are needed to form
675    # a non-singular square matrix
676    # If all are varied and this is a constraint equation, then set VaryFree flag
677    # so that the newly created relationships will be varied
678    for group,varlist in zip(groups,parmlist):
679        if len(varlist) == 1: continue
680        # for constraints, if all included variables are refined,
681        # set the VaryFree flag, and remaining degrees of freedom will be
682        # varied (since consistency was checked, if any one variable is
683        # refined, then assume that all are)
684        varsList = [] # make a list of all the referenced variables as well
685        VaryFree = False
686        for rel in group:
687            varied = 0
688            unused = 0
689            for var in VarKeys(constrDict[rel]):
690                var = translateTable.get(var,var) # replace wildcards
691                if parmDict is not None and var not in parmDict:
692                    unused += 1                   
693                if var not in varsList: varsList.append(var)
694                if var in varyList: varied += 1
695            if fixedList[rel] is not None and varied > 0:
696                VaryFree = True
697        if len(varlist) < len(group): # too many relationships -- no can do
698            msg = 'too many relationships'
699            break
700        # Since we checked before, if any variables are unused, then all must be.
701        # If so, this set of relationships can be ignored
702        if unused:
703            if debug: print('Constraint ignored (all variables undefined)')
704            if debug: print ('    '+_FormatConstraint(constrDict[rel],fixedList[rel]))
705            continue
706        # fill in additional degrees of freedom
707        try:
708            arr = _FillArray(group,constrDict,varlist)
709            _RowEchelon(len(group),arr,varlist)
710            constrArr = _FillArray(group,constrDict,varlist,FillDiagonals=True)
711            GramSchmidtOrtho(constrArr,len(group))
712        except:
713            msg = 'Singular relationships'
714            break
715        mapvar = []
716        group = group[:]
717        # scan through all generated and input relationships, we need to add to the varied list
718        # all the new parameters where VaryFree has been set or where a New Var is varied.
719        #
720        # If a group does not contain any fixed values (constraint equations)
721        # and nothing in the group is varied, drop this group, so that the
722        # dependent parameters can be refined individually.
723        unused = True
724        for i in range(len(varlist)):
725            if len(group) > 0: # get the original equation reference
726                rel = group.pop(0)
727                fixedval = fixedList[rel]
728                varyflag = constrDict[rel].get('_vary',False)
729                varname = constrDict[rel].get('_name','')
730            else: # this relationship has been generated
731                varyflag = False
732                varname = ''
733                fixedval = None
734            if fixedval is None: # this is a new variable, not a constraint
735                if not varname:
736                    varname = paramPrefix + str(consNum) # no assigned name, create one
737                    consNum += 1
738                mapvar.append(varname)
739                # vary the new relationship if it is a degree of freedom in
740                # a set of contraint equations or if a New Var is flagged to be varied.
741                if VaryFree or varyflag: 
742                    unused = False
743                    varyList.append(varname)
744                    # fix (prevent varying) of all the variables inside the constraint group
745                    # (dependent vars)
746                    for var in varsList:
747                        if var in varyList: varyList.remove(var)
748            else:
749                unused = False
750                mapvar.append(fixedval)
751        if unused: # skip over constraints that don't matter (w/o fixed value or any refined variables)
752            if debug: print('Constraint ignored (all variables unrefined)')
753            if debug: print ('   '+_FormatConstraint(constrDict[rel],fixedList[rel]))
754            continue 
755        dependentParmList.append([translateTable.get(var,var) for var in varlist])
756        arrayList.append(constrArr)
757        invarrayList.append(np.linalg.inv(constrArr))
758        indParmList.append(mapvar)
759        symGenList.append(False)
760    if msg:
761        print (' *** ERROR in constraint definitions! ***')
762        print (msg)
763        print (VarRemapShow(varyList))
764        raise Exception
765    # setup dictionary containing the fixed values
766    global fixedDict
767    # key is original ascii string, value is float
768    for fixedval in fixedList:
769        if fixedval:
770            fixedDict[fixedval] = float(fixedval)
771    _setVarLists(dropVarList)
772   
773def _setVarLists(dropVarList):
774    '''Make list of dependent and independent variables (after dropping unused vars in dropVarList)
775    '''
776    global dependentParmList,indParmList
777    global dependentVars
778    global independentVars
779    dependentVars = []
780    independentVars = []
781    for varlist,mapvars in zip(dependentParmList,indParmList):  # process all constraints
782        for mv in mapvars:
783            if mv in dropVarList: continue
784            if mv not in independentVars: independentVars.append(mv)
785        for mv in varlist:
786            if mv in dropVarList: continue
787            if mv not in dependentVars: dependentVars.append(mv)
788    if debug: # on debug, show what is parsed & generated, semi-readable
789        print (50*'-')
790        print (VarRemapShow(varyList))
791        print ('Varied: ',varyList)
792        print ('Not Varied: ',fixedVarList)
793
794def StoreEquivalence(independentVar,dependentList,symGen=True):
795    '''Takes a list of dependent parameter(s) and stores their
796    relationship to a single independent parameter (independentVar)
797
798    :param str independentVar: name of master parameter that will be used to determine the value
799      to set the dependent variables
800
801    :param list dependentList: a list of parameters that will set from
802         independentVar. Each item in the list can be a string with the parameter
803         name or a tuple containing a name and multiplier:
804         ``['::parm1',('::parm2',.5),]``
805
806    '''
807   
808    global dependentParmList,arrayList,invarrayList,indParmList,symGenList
809    mapList = []
810    multlist = []
811    for var in dependentList:
812        if isinstance(var, str):
813            mult = 1.0
814        elif len(var) == 2:
815            var,mult = var
816        else:
817            raise Exception("Cannot parse "+repr(var) + " as var or (var,multiplier)")
818        mapList.append(var)
819        multlist.append(tuple((mult,)))
820    # added relationships to stored values
821    arrayList.append(None)
822    invarrayList.append(np.array(multlist))
823    indParmList.append(tuple((independentVar,)))
824    dependentParmList.append(mapList)
825    symGenList.append(symGen)
826    return
827
828def GetDependentVars():
829    '''Return a list of dependent variables: e.g. variables that are
830    constrained in terms of other variables
831
832    :returns: a list of variable names
833
834    '''
835    global dependentVars
836    return dependentVars
837
838def GetIndependentVars():
839    '''Return a list of independent variables: e.g. variables that are
840    created by constraints of other variables
841
842    :returns: a list of variable names
843
844    '''
845    global independentVars
846    return independentVars
847
848def PrintIndependentVars(parmDict,varyList,sigDict,PrintAll=False,pFile=None):
849    '''Print the values and uncertainties on the independent variables'''
850    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
851    printlist = []
852    mapvars = GetIndependentVars()
853    for i,name in enumerate(mapvars):
854        if name in fixedDict: continue
855        if PrintAll or name in varyList:
856            sig = sigDict.get(name)
857            printlist.append([name,parmDict[name],sig])
858    if len(printlist) == 0: return
859    s1 = ''
860    s2 = ''
861    s3 = ''
862    pFile.write(130*'-'+'\n')
863    pFile.write("Variables generated by constraints\n")
864    printlist.append(3*[None])
865    for name,val,esd in printlist:
866        if len(s1) > 120 or name is None:
867            pFile.write(''+'\n')
868            pFile.write(s1+'\n')
869            pFile.write(s2+'\n')
870            pFile.write(s3+'\n')
871            s1 = ''
872            if name is None: break
873        if s1 == "":
874            s1 = ' name  :'
875            s2 = ' value :'
876            s3 = ' sig   :'
877        s1 += '%15s' % (name)
878        s2 += '%15.4f' % (val)
879        if esd is None:
880            s3 += '%15s' % ('n/a')
881        else:   
882            s3 += '%15.4f' % (esd)
883
884def ComputeDepESD(covMatrix,varyList,parmDict):
885    '''Compute uncertainties for dependent parameters from independent ones
886    returns a dictionary containing the esd values for dependent parameters
887    '''
888    sigmaDict = {}
889    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
890        #if invmultarr is None: continue # probably not needed
891#        try:
892#            valuelist = [parmDict[var] for var in mapvars]
893#        except KeyError:
894#            continue
895        # get the v-covar matrix for independent parameters
896        vcov = np.zeros((len(mapvars),len(mapvars)))
897        for i1,name1 in enumerate(mapvars):
898            if name1 not in varyList: continue
899            iv1 = varyList.index(name1)
900            for i2,name2 in enumerate(mapvars):
901                if name2 not in varyList: continue
902                iv2 = varyList.index(name2)
903                vcov[i1][i2] = covMatrix[iv1][iv2]
904        # vec is the vector that multiplies each of the independent values
905        for v,vec in zip(varlist,invmultarr):
906            sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec)))
907    return sigmaDict
908
909def _FormatConstraint(RelDict,RelVal):
910    '''Formats a Constraint or Function for use in a convenient way'''
911    linelen = 45
912    s = [""]
913    for var,val in RelDict.items():
914        if var.startswith('_'): continue
915        if len(s[-1]) > linelen: s.append(' ')
916        m = val
917        if s[-1] != "" and m >= 0:
918            s[-1] += ' + '
919        elif s[-1] != "":
920            s[-1] += ' - '
921            m = abs(m)
922        s[-1] += '%.3f*%s '%(m,var)
923    if len(s[-1]) > linelen: s.append(' ')
924    if RelVal is None:
925        s[-1] += ' = New variable'
926    else:
927        s[-1] += ' = ' + RelVal
928    s1 = ''
929    for s2 in s:
930        if s1 != '': s1 += '\n\t'
931        s1 += s2
932    return s1
933
934def VarRemapShow(varyList,inputOnly=False):
935    '''List out the saved relationships. This should be done after the constraints have been
936    defined using :func:`StoreEquivalence`, :func:`GroupConstraints` and :func:`GenerateConstraints`.
937
938    :returns: a string containing the details of the contraint relationships
939    '''
940    s = ''
941    if len(fixedVarList) > 0:
942        s += 'Fixed Variables:\n'
943        for v in fixedVarList:
944            s += '    ' + v + '\n'
945    s += 'User-supplied variable mapping relations:\n'
946    symout = ''
947    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,symGenList
948
949    for varlist,mapvars,multarr,invmultarr,symFlag in zip(
950        dependentParmList,indParmList,arrayList,invarrayList,symGenList):
951        for i,mv in enumerate(mapvars):
952            if multarr is None:
953#                s1 = '  ' + str(mv) + ' is equivalent to parameter(s): '
954                if len(varlist) == 1:
955                    s1 = '   ' + str(mv) + ' is equivalent to '
956                else:
957                    s1 = '   ' + str(mv) + ' is equivalent to parameters: '
958                j = 0
959                for v,m in zip(varlist,invmultarr):
960                    if debug: print ('v,m[0]: ',v,m[0])
961                    if len(s1.split('\n')[-1]) > 75: s1 += '\n        '
962                    if j > 0: s1 += ' & '
963                    j += 1
964                    s1 += str(v)
965                    if m != 1:
966                        s1 += " / " + str(m[0])
967                if symFlag:
968                    symout += s1 + '\n'
969                else:
970                    s += s1 + '\n'
971                continue
972            s += %s = ' % mv
973            j = 0
974            for m,v in zip(multarr[i,:],varlist):
975                if m == 0: continue
976                if j > 0: s += ' + '
977                j += 1
978                s += '(%s * %s)' % (m,v)
979            if mv in varyList: s += ' VARY'
980            s += '\n'
981    if symout:
982        s += 'Symmetry-generated relations:\n' + symout
983    if inputOnly: return s
984    s += 'Inverse variable mapping relations:\n'
985    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
986        for i,mv in enumerate(varlist):
987            s += %s = ' % mv
988            j = 0
989            for m,v in zip(invmultarr[i,:],mapvars):
990                if m == 0: continue
991                if j > 0: s += ' + '
992                j += 1
993                s += '(%s * %s)' % (m,v)
994            s += '\n'
995    return s
996
997def GetSymEquiv():
998    '''Return the automatically generated (equivalence) relationships.
999
1000    :returns: a list of strings containing the details of the contraint relationships
1001    '''
1002    symout = []
1003    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict,symGenList
1004
1005    for varlist,mapvars,multarr,invmultarr,symFlag in zip(
1006        dependentParmList,indParmList,arrayList,invarrayList,symGenList):
1007        for i,mv in enumerate(mapvars):
1008            if not symFlag: continue
1009            if multarr is None:
1010                s1 = str(mv) + ' = '
1011                j = 0
1012                for v,m in zip(varlist,invmultarr):
1013                    if debug: print ('v,m[0]: ',v,m[0])
1014                    if len(s1.split('\n')[-1]) > 75: s1 += '\n        '
1015                    if j > 0: s1 += ' =  '
1016                    j += 1
1017                    s1 += str(v)
1018                    if m != 1:
1019                        s1 += " / " + str(m[0])
1020                    symout.append(s1)
1021                continue
1022            else:
1023                s = %s = ' % mv
1024                j = 0
1025                for m,v in zip(multarr[i,:],varlist):
1026                    if m == 0: continue
1027                    if j > 0: s += ' + '
1028                    j += 1
1029                    s += '(%s * %s)' % (m,v)
1030                print ('unexpected sym op='+s)
1031    return symout
1032
1033def Dict2Deriv(varyList,derivDict,dMdv):
1034    '''Compute derivatives for Independent Parameters from the
1035    derivatives for the original parameters
1036
1037    :param list varyList: a list of parameters names that will be varied
1038
1039    :param dict derivDict: a dict containing derivatives for parameter values keyed by the
1040      parameter names.
1041
1042    :param list dMdv: a Jacobian, as a list of np.array containing derivatives for dependent
1043      parameter computed from derivDict
1044
1045    '''
1046    global dependentParmList,arrayList,invarrayList,indParmList,invarrayList
1047    for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
1048        for i,name in enumerate(mapvars):
1049            # grouped variables: need to add in the derv. w/r
1050            # dependent variables to the independent ones
1051            if name not in varyList: continue # skip if independent var not varied
1052            if multarr is None:
1053                for v,m in zip(varlist,invmultarr):
1054                    if debug: print ('start dMdv',dMdv[varyList.index(name)])
1055                    if debug: print ('add derv',v,'/',m[0],'to derv',name,'add=',derivDict[v] / m[0])
1056                    if m == 0: continue
1057                    dMdv[varyList.index(name)] += derivDict[v] / m[0]
1058            else:
1059                for v,m in zip(varlist,multarr[i,:]):
1060                    if debug: print ('start dMdv',dMdv[varyList.index(name)])
1061                    if debug: print ('add derv',v,'*',m,'to derv',name,'add=',m * derivDict[v])
1062                    if m == 0: continue
1063                    dMdv[varyList.index(name)] += m * derivDict[v]
1064
1065def Map2Dict(parmDict,varyList):
1066    '''Create (or update) the Independent Parameters from the original
1067    set of Parameters
1068
1069    Removes dependent variables from the varyList
1070
1071    This should be done once, after the constraints have been
1072    defined using :func:`StoreEquivalence`,
1073    :func:`GroupConstraints` and :func:`GenerateConstraints` and
1074    before any variable refinement is done. This completes the parameter
1075    dictionary by defining independent parameters and it satisfies the
1076    constraint equations in the initial parameters
1077
1078    :param dict parmDict: a dict containing parameter values keyed by the
1079      parameter names.
1080      This will contain updated values for both dependent and independent
1081      parameters after Dict2Map is called. It will also contain some
1082      unexpected entries of every constant value {'0':0.0} & {'1.0':1.0},
1083      which do not cause any problems.
1084
1085    :param list varyList: a list of parameters names that will be varied
1086   
1087
1088    '''
1089    # process the Independent vars: remove dependent ones from varylist
1090    # and then compute values for the Independent ones from their dependents
1091    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
1092    for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
1093        for item in varlist:
1094            try:
1095                varyList.remove(item)
1096            except ValueError:
1097                pass
1098        if multarr is None: continue
1099        valuelist = [parmDict[var] for var in varlist]
1100        parmDict.update(zip(mapvars,
1101                            np.dot(multarr,np.array(valuelist)))
1102                        )
1103    # now remove fixed variables from the varyList
1104    global fixedVarList
1105    for item in fixedVarList:
1106        try:
1107            varyList.remove(item)
1108        except ValueError:
1109            pass
1110    # Set constrained parameters to their fixed values
1111    parmDict.update(fixedDict)
1112
1113def Dict2Map(parmDict,varyList):
1114    '''Applies the constraints defined using :func:`StoreEquivalence`,
1115    :func:`GroupConstraints` and :func:`GenerateConstraints` by changing
1116    values in a dict containing the parameters. This should be
1117    done before the parameters are used for any computations
1118
1119    :param dict parmDict: a dict containing parameter values keyed by the
1120      parameter names.
1121      This will contain updated values for both dependent and independent
1122      parameters after Dict2Map is called. It will also contain some
1123      unexpected entries of every constant value {'0':0.0} & {'1.0':1.0},
1124      which do not cause any problems.
1125
1126    :param list varyList: a list of parameters names that will be varied
1127   
1128    '''
1129    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
1130    # reset fixed values (should not be needed, but very quick)
1131    # - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended
1132    # not needed, but also not a problem - BHT
1133    parmDict.update(fixedDict)
1134    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
1135        #if invmultarr is None: continue
1136        try: 
1137            valuelist = [parmDict[var] for var in mapvars]
1138        except KeyError:
1139            continue
1140        parmDict.update(zip(varlist,np.dot(invmultarr,np.array(valuelist))))
1141
1142#======================================================================
1143# internal routines follow (these routines are unlikely to be called
1144# from outside the module)
1145
1146def GramSchmidtOrtho(a,nkeep=0):
1147    '''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to
1148    find orthonormal unit vectors relative to first row.
1149
1150    If nkeep is non-zero, the first nkeep rows in the array are not changed
1151   
1152    input:
1153       arrayin: a 2-D non-singular square array
1154    returns:
1155       a orthonormal set of unit vectors as a square array
1156    '''
1157    def proj(a,b):
1158        'Projection operator'
1159        return a*(np.dot(a,b)/np.dot(a,a))
1160    for j in range(nkeep,len(a)):
1161        for i in range(j):
1162            a[j] -= proj(a[i],a[j])
1163        if np.allclose(np.linalg.norm(a[j]),0.0):
1164            raise Exception("Singular input to GramSchmidtOrtho")
1165        a[j] /= np.linalg.norm(a[j])
1166    return a
1167
1168def _FillArray(sel,dict,collist,FillDiagonals=False):
1169    '''Construct a n by n matrix (n = len(collist)
1170    filling in the rows using the relationships defined
1171    in the dictionaries selected by sel
1172
1173    If FillDiagonals is True, diagonal elements in the
1174    array are set to 1.0
1175    '''
1176    n = len(collist)
1177    if FillDiagonals:
1178        arr = np.eye(n)
1179    else:
1180        arr = np.zeros(2*[n,])
1181    # fill the top rows
1182    for i,cnum in enumerate(sel):
1183        for j,var in enumerate(collist):
1184            arr[i,j] = dict[cnum].get(var,0)
1185    return arr
1186
1187def _SwapColumns(i,m,v):
1188    '''Swap columns in matrix m as well as the labels in v
1189    so that element (i,i) is replaced by the first non-zero element in row i after that element
1190
1191    Throws an exception if there are no non-zero elements in that row
1192    '''
1193    for j in range(i+1,len(v)):
1194        if not np.allclose(m[i,j],0):
1195            m[:,(i,j)] = m[:,(j,i)]
1196            v[i],v[j] = v[j],v[i]
1197            return
1198    else:
1199        raise Exception('Singular input')
1200
1201def _RowEchelon(m,arr,collist):
1202    '''Convert the first m rows in Matrix arr to row-echelon form
1203    exchanging columns in the matrix and collist as needed.
1204
1205    throws an exception if the matrix is singular because
1206    the first m rows are not linearly independent
1207    '''
1208    for i in range(m):
1209        if np.allclose(arr[i,i],0):
1210            _SwapColumns(i,arr,collist)
1211        arr[i,:] /= arr[i,i] # normalize row
1212        # subtract current row from subsequent rows to set values to left of diagonal to 0
1213        for j in range(i+1,m):
1214            arr[j,:] -= arr[i,:] * arr[j,i]
1215
1216if __name__ == "__main__":
1217    parmdict = {}
1218    constrDict = [
1219        {'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},
1220        {'0:0:eA': 0.0},
1221        {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
1222        {'1::C(10,0,1)': 1.0, '2::C(10,0,1)': 1.0},
1223        {'1::AUiso:0': 1.0, '0::AUiso:0': 1.0},
1224        {'0::A0': 0.0}
1225        ]
1226    fixedList = ['5.0', '0', None, None, '1.0', '0']
1227    StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), ))
1228    #StoreEquivalence('1::atomx:3',('2::atomx:3', ('2::atomz:3',2,), )) # error: dependent & independent vars mixed
1229    #StoreEquivalence('1::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) # error: dependent vars repeated
1230    #StoreEquivalence('0:1:eA',('0:0:eA',)) # error: equiv & fixed
1231    #StoreEquivalence('0:99:Scale',('0:12:Scale',)) # error: equiv & constrained
1232    #StoreEquivalence('0:12:Scale',('0:99:Scale',)) # error: equiv & constrained
1233    varylist = ['2::atomx:3',
1234                '2::C(10,6,1)', '1::C(10,6,1)',
1235                '2::atomy:3', '2::atomz:3',
1236                '0:12:Scale', '0:11:Scale', '0:14:Scale', '0:13:Scale', '0:0:Scale']
1237#    e,w = CheckConstraints([,
1238#                     [{'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}],
1239#                     ['1.0'])
1240#    if e: print 'error=',e
1241#    if w: print 'error=',w
1242#    varyList = ['0::A0', '0::AUiso:0', '0::Afrac:1', '0::Afrac:2', '0::Afrac:3', '0::Afrac:4', '0::dAx:5', '0::dAy:5', '0::dAz:5', '0::AUiso:5', ':0:Back;0', ':0:Back;1', ':0:Back;2', ':0:Back;3', ':0:Back;4', ':0:Back;5', ':0:Back;6', ':0:Back;7', ':0:Back;8', ':0:Back;9', ':0:Back;10', ':0:Back;11', ':0:U', ':0:V', ':0:W', ':0:X', ':0:Y', ':0:Scale', ':0:DisplaceX', ':0:DisplaceY']
1243#    constrDict = [
1244#        {'0::Afrac:4': 24.0, '0::Afrac:1': 16.0, '0::Afrac:3': 24.0, '0::Afrac:2': 16.0},
1245#        {'0::Afrac:1': 1.0, '0::Afrac:2': 1.0},
1246#        {'0::Afrac:4': 1.0, '0::Afrac:3': 1.0}]
1247#    fixedList = ['40.0', '1.0', '1.0']
1248
1249    errmsg, warnmsg = CheckConstraints(varylist,constrDict,fixedList)
1250    if errmsg:
1251        print ("*** Error ********************")
1252        print (errmsg)
1253    if warnmsg:
1254        print ("*** Warning ********************")
1255        print (warnmsg)
1256    if errmsg or warnmsg:
1257        sys.exit()
1258    groups,parmlist = GroupConstraints(constrDict)
1259    GenerateConstraints(groups,parmlist,varylist,constrDict,fixedList)
1260    print (VarRemapShow(varylist))
1261    parmdict.update( {
1262        '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,
1263        '0:0:eA': 0.0,
1264        '2::C(10,6,1)': 0.2, '1::C(10,6,1)': 0.3,
1265        '1::C(10,0,1)': 0.2, '2::C(10,0,1)': 0.3,
1266        '1::AUiso:0': 0.02, '0::AUiso:0': 0.03,
1267        '0::A0': 0.0,
1268        '2::atomx:3':0.23,'2::atomy:3':-.23, '2::atomz:3':-0.11,
1269        })
1270    print ('parmdict start',parmdict)
1271    print ('varylist start',varylist)
1272    before = parmdict.copy()
1273    Map2Dict(parmdict,varylist)
1274    print ('parmdict before and after Map2Dict')
1275    print ('  key / before / after')
1276    for key in sorted(list(parmdict.keys())):
1277        print ('  '+key,'\t',before.get(key),'\t',parmdict[key])
1278    print ('varylist after',varylist)
1279    before = parmdict.copy()
1280    Dict2Map(parmdict,varylist)
1281    print ('after Dict2Map')
1282    print ('  key / before / after')
1283    for key in sorted(list(parmdict.keys())):
1284        print ('  '+key,'\t',before.get(key),'\t',parmdict[key])
1285#    dMdv = len(varylist)*[0]
1286#    deriv = {}
1287#    for i,v in enumerate(parmdict.keys()): deriv[v]=i
1288#    Dict2Deriv(varylist,deriv,dMdv)
Note: See TracBrowser for help on using the repository browser.