source: trunk/GSASIImapvars.py @ 3046

Last change on this file since 3046 was 3046, checked in by toby, 4 years ago

label user-supplied & code-generated equivalences (constraints); improve error conflicting constraint error message

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