source: branch/2frame/GSASIImapvars.py @ 2997

Last change on this file since 2997 was 2895, checked in by toby, 8 years ago

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