source: trunk/GSASIImapvars.py @ 1160

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

finish ISODISPLACE fixes; improve show var window; improve help window; add refine checkbox for newvars in constraints display

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