source: trunk/GSASIImapvars.py @ 1143

Last change on this file since 1143 was 1143, checked in by toby, 9 years ago

rework constraints to handle names and refine flag for new var (input linear constraints); redo powder 'Sample Parameters'
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 47.1 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2013-11-11 02:32:26 +0000 (Mon, 11 Nov 2013) $
4# $Author: toby $
5# $Revision: 1143 $
6# $URL: trunk/GSASIImapvars.py $
7# $Id: GSASIImapvars.py 1143 2013-11-11 02:32:26Z 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: 1143 $")
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:
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    # 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 all the
674        # dependent parameters are varied. Check again for inconsistent variable use
675        # for new variables -- where varied and unvaried parameters get grouped
676        # together. I don't think this can happen when not flagged before, but
677        # it does not hurt to check again.
678        for i in range(len(varlist)):
679            if len(group) > 0: # get the original equation reference
680                rel = group.pop(0)
681                fixedval = fixedList[rel]
682                varyflag = constrDict[rel].get('_vary',False)
683                varname = constrDict[rel].get('_name','')
684            else: # relationship has been generated
685                varyflag = False
686                varname = ''
687                fixedval = None
688            if fixedval is None: # this is a new variable, not a constraint
689                if not varname:
690                    varname = paramPrefix + str(consNum) # no assigned name, create one
691                    consNum += 1
692                mapvar.append(varname)
693                # vary the new relationship if it is a degree of freedom in
694                # a set of contraint equations or if a new variable is flagged to be varied.
695                if VaryFree or varyflag: 
696                    varyList.append(varname)
697                    # fix (prevent varying) of all the variables inside the constraint group
698                    for var in varsList:
699                        if var in varyList: varyList.remove(var)
700            else:
701                mapvar.append(fixedval)
702        dependentParmList.append(varlist)
703        arrayList.append(constrArr)
704        invarrayList.append(np.linalg.inv(constrArr))
705        indParmList.append(mapvar)
706    if msg:
707        print ' *** ERROR in constraint definitions! ***'
708        print msg
709        print VarRemapShow(varyList)
710        raise Exception
711    # setup dictionary containing the fixed values
712    global fixedDict
713    # key is original ascii string, value is float
714    for fixedval in fixedList:
715        if fixedval:
716            fixedDict[fixedval] = float(fixedval)
717
718    if debug: # on debug, show what is parsed & generated, semi-readable
719        print 50*'-'
720        print VarRemapShow(varyList)
721        print 'Varied: ',varyList
722        print 'Not Varied: ',fixedVarList
723
724def StoreEquivalence(independentVar,dependentList):
725    '''Takes a list of dependent parameter(s) and stores their
726    relationship to a single independent parameter (independentVar)
727
728    :param str independentVar: name of master parameter that will be used to determine the value
729      to set the dependent variables
730
731    :param list dependentList: a list of parameters that will set from
732         independentVar. Each item in the list can be a string with the parameter
733         name or a tuple containing a name and multiplier:
734         ``['parm1',('parm2',.5),]``
735
736    '''
737   
738    global dependentParmList,arrayList,invarrayList,indParmList
739    mapList = []
740    multlist = []
741    for var in dependentList:
742        if isinstance(var, basestring):
743            mult = 1.0
744        elif len(var) == 2:
745            var,mult = var
746        else:
747            raise Exception, "Cannot parse "+repr(var) + " as var or (var,multiplier)"
748        mapList.append(var)
749        multlist.append(tuple((mult,)))
750    # added relationships to stored values
751    arrayList.append(None)
752    invarrayList.append(np.array(multlist))
753    indParmList.append(tuple((independentVar,)))
754    dependentParmList.append(mapList)
755    return
756
757def GetDependentVars():
758    '''Return a list of dependent variables: e.g. variables that are
759    constrained in terms of other variables
760
761    :returns: a list of variable names
762
763    '''
764    dependentVars = []
765    global dependentParmList
766    for lst in dependentParmList:
767        for itm in lst: dependentVars.append(itm)
768    return dependentVars
769
770def GetIndependentVars():
771    '''Return a list of independent variables: e.g. variables that are
772    created by constraints of other variables
773
774    :returns: a list of variable names
775
776    '''
777    independentVars = []
778    global indParmList,fixedDict
779    for lst in indParmList:
780        for name in lst:
781            if name in fixedDict: continue
782            independentVars.append(name)
783    return independentVars
784
785def PrintIndependentVars(parmDict,varyList,sigDict,PrintAll=False,pFile=None):
786    '''Print the values and uncertainties on the independent variables'''
787    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
788    printlist = []
789    mapvars = GetIndependentVars()
790    for i,name in enumerate(mapvars):
791        if name in fixedDict: continue
792        if PrintAll or name in varyList:
793            sig = sigDict.get(name)
794            printlist.append([name,parmDict[name],sig])
795    if len(printlist) == 0: return
796    s1 = ''
797    print >>pFile,130*'-'
798    print >>pFile,"Variables generated by constraints"
799    printlist.append(3*[None])
800    for name,val,esd in printlist:
801        if len(s1) > 120 or name is None:
802            print >>pFile,''
803            print >>pFile,s1
804            print >>pFile,s2
805            print >>pFile,s3
806            s1 = ''
807            if name is None: break
808        if s1 == "":
809            s1 = ' name  :'
810            s2 = ' value :'
811            s3 = ' sig   :'
812        s1 += '%15s' % (name)
813        s2 += '%15.4f' % (val)
814        if esd is None:
815            s3 += '%15s' % ('n/a')
816        else:   
817            s3 += '%15.4f' % (esd)
818
819def ComputeDepESD(covMatrix,varyList,parmDict):
820    '''Compute uncertainties for dependent parameters from independent ones
821    returns a dictionary containing the esd values for dependent parameters
822    '''
823    sigmaDict = {}
824    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
825        #if invmultarr is None: continue # probably not needed
826        valuelist = [parmDict[var] for var in mapvars]
827        # get the v-covar matrix for independent parameters
828        vcov = np.zeros((len(mapvars),len(mapvars)))
829        for i1,name1 in enumerate(mapvars):
830            if name1 not in varyList: continue
831            iv1 = varyList.index(name1)
832            for i2,name2 in enumerate(mapvars):
833                if name2 not in varyList: continue
834                iv2 = varyList.index(name2)
835                vcov[i1][i2] = covMatrix[iv1][iv2]
836        # vec is the vector that multiplies each of the independent values
837        for v,vec in zip(varlist,invmultarr):
838            sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec)))
839    return sigmaDict
840
841def _FormatConstraint(RelDict,RelVal):
842    '''Formats a Constraint or Function for use in a convenient way'''
843    linelen = 45
844    s = [""]
845    for var,val in RelDict.items():
846        if var.startswith('_'): continue
847        if len(s[-1]) > linelen: s.append(' ')
848        m = val
849        if s[-1] != "" and m >= 0:
850            s[-1] += ' + '
851        elif s[-1] != "":
852            s[-1] += ' - '
853            m = abs(m)
854        s[-1] += '%.3f*%s '%(m,var)
855    if len(s[-1]) > linelen: s.append(' ')
856    if RelVal is None:
857        s[-1] += ' = New variable'
858    else:
859        s[-1] += ' = ' + RelVal
860    s1 = ''
861    for s2 in s:
862        if s1 != '': s1 += '\n\t'
863        s1 += s2
864    return s1
865
866def VarRemapShow(varyList):
867    '''List out the saved relationships. This should be done after the constraints have been
868    defined using :func:`StoreEquivalence`, :func:`GroupConstraints` and :func:`GenerateConstraints`.
869
870    :returns: a string containing the details of the contraint relationships
871    '''
872    s = ''
873    if len(fixedVarList) > 0:
874        s += 'Fixed Variables:\n'
875        for v in fixedVarList:
876            s += '    ' + v + '\n'
877    s += 'Variable mapping relations:\n'
878    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
879    for varlist,mapvars,multarr,invmultarr in zip(
880        dependentParmList,indParmList,arrayList,invarrayList):
881        for i,mv in enumerate(mapvars):
882            if multarr is None:
883                s += '  ' + str(mv) + ' is equivalent to parameter(s): '
884                j = 0
885                for v,m in zip(varlist,invmultarr):
886                    #print v,m[0]
887                    if j > 0: s += '  & '
888                    j += 1
889                    s += str(v)
890                    if m != 1:
891                        s += " / " + str(m[0])                       
892                s += '\n'
893                continue
894            s += '  %s = ' % mv
895            j = 0
896            for m,v in zip(multarr[i,:],varlist):
897                if m == 0: continue
898                if j > 0: s += ' + '
899                j += 1
900                s += '(%s * %s)' % (m,v)
901            if mv in varyList: s += ' VARY'
902            s += '\n'
903    s += 'Inverse variable mapping relations:\n'
904    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
905        for i,mv in enumerate(varlist):
906            s += '  %s = ' % mv
907            j = 0
908            for m,v in zip(invmultarr[i,:],mapvars):
909                if m == 0: continue
910                if j > 0: s += ' + '
911                j += 1
912                s += '(%s * %s)' % (m,v)
913            s += '\n'
914    return s
915
916def Dict2Deriv(varyList,derivDict,dMdv):
917    '''Compute derivatives for Independent Parameters from the
918    derivatives for the original parameters
919
920    :param list varyList: a list of parameters names that will be varied
921
922    :param dict derivDict: a dict containing derivatives for parameter values keyed by the
923      parameter names.
924
925    :param list dMdv: a Jacobian, as a list of np.array containing derivatives for dependent
926      parameter computed from derivDict
927
928    '''
929    global dependentParmList,arrayList,invarrayList,indParmList,invarrayList
930    for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
931        for i,name in enumerate(mapvars):
932            # grouped variables: need to add in the derv. w/r
933            # dependent variables to the independent ones
934            if name not in varyList: continue # skip if independent var not varied
935            if multarr is None:
936                for v,m in zip(varlist,invmultarr):
937                    #print 'add derv',v,'/',m[0],'to derv',name
938                    if m == 0: continue
939                    dMdv[varyList.index(name)] += derivDict[v] / m[0]
940            else:
941                for v,m in zip(varlist,multarr[i,:]):
942                    #print 'add derv',v,'*',m,'to derv',name
943                    if m == 0: continue
944                    dMdv[varyList.index(name)] += m * derivDict[v]
945
946def Map2Dict(parmDict,varyList):
947    '''Create (or update) the Independent Parameters from the original
948    set of Parameters
949
950    Removes dependent variables from the varyList
951
952    This should be done once, after the constraints have been
953    defined using :func:`StoreEquivalence`,
954    :func:`GroupConstraints` and :func:`GenerateConstraints` and
955    before any variable refinement is done
956    to complete the parameter dictionary by defining independent
957    parameters and satisfying the constraint equations.
958
959    :param dict parmDict: a dict containing parameter values keyed by the
960      parameter names.
961      This will contain updated values for both dependent and independent
962      parameters after Dict2Map is called. It will also contain some
963      unexpected entries of every constant value {'0':0.0} & {'1.0':1.0},
964      which do not cause any problems.
965
966    :param list varyList: a list of parameters names that will be varied
967   
968
969    '''
970    # process the Independent vars: remove dependent ones from varylist
971    # and then compute values for the Independent ones from their dependents
972    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
973    for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
974        for item in varlist:
975            try:
976                varyList.remove(item)
977            except ValueError:
978                pass
979        if multarr is None: continue
980        valuelist = [parmDict[var] for var in varlist]
981        parmDict.update(zip(mapvars,
982                            np.dot(multarr,np.array(valuelist)))
983                        )
984    # now remove fixed variables from the varyList
985    global fixedVarList
986    for item in fixedVarList:
987        try:
988            varyList.remove(item)
989        except ValueError:
990            pass
991    # Set constrained parameters to their fixed values
992    parmDict.update(fixedDict)
993
994def Dict2Map(parmDict,varyList):
995    '''Applies the constraints defined using :func:`StoreEquivalence`,
996    :func:`GroupConstraints` and :func:`GenerateConstraints` by changing
997    values in a dict containing the parameters. This should be
998    done before the parameters are used for any computations
999
1000    :param dict parmDict: a dict containing parameter values keyed by the
1001      parameter names.
1002      This will contain updated values for both dependent and independent
1003      parameters after Dict2Map is called. It will also contain some
1004      unexpected entries of every constant value {'0':0.0} & {'1.0':1.0},
1005      which do not cause any problems.
1006
1007    :param list varyList: a list of parameters names that will be varied
1008   
1009    '''
1010    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
1011    # reset fixed values (should not be needed, but very quick)
1012    # - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended
1013    # not needed, but also not a problem - BHT
1014    parmDict.update(fixedDict)
1015    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
1016        #if invmultarr is None: continue
1017        valuelist = [parmDict[var] for var in mapvars]
1018        parmDict.update(zip(varlist,
1019                            np.dot(invmultarr,np.array(valuelist)))
1020                        )
1021
1022#======================================================================
1023# internal routines follow (these routines are unlikely to be called
1024# from outside the module)
1025
1026def GramSchmidtOrtho(a,nkeep=0):
1027    '''Use the Gram-Schmidt process (http://en.wikipedia.org/wiki/Gram-Schmidt) to
1028    find orthonormal unit vectors relative to first row.
1029
1030    If nkeep is non-zero, the first nkeep rows in the array are not changed
1031   
1032    input:
1033       arrayin: a 2-D non-singular square array
1034    returns:
1035       a orthonormal set of unit vectors as a square array
1036    '''
1037    def proj(a,b):
1038        'Projection operator'
1039        return a*(np.dot(a,b)/np.dot(a,a))
1040    for j in range(nkeep,len(a)):
1041        for i in range(j):
1042            a[j] -= proj(a[i],a[j])
1043        if np.allclose(np.linalg.norm(a[j]),0.0):
1044            raise Exception,"Singular input to GramSchmidtOrtho"
1045        a[j] /= np.linalg.norm(a[j])
1046    return a
1047
1048def _FillArray(sel,dict,collist,FillDiagonals=False):
1049    '''Construct a n by n matrix (n = len(collist)
1050    filling in the rows using the relationships defined
1051    in the dictionaries selected by sel
1052
1053    If FillDiagonals is True, diagonal elements in the
1054    array are set to 1.0
1055    '''
1056    n = len(collist)
1057    if FillDiagonals:
1058        arr = np.eye(n)
1059    else:
1060        arr = np.zeros(2*[n,])
1061    # fill the top rows
1062    for i,cnum in enumerate(sel):
1063        for j,var in enumerate(collist):
1064            arr[i,j] = dict[cnum].get(var,0)
1065    return arr
1066
1067def _SwapColumns(i,m,v):
1068    '''Swap columns in matrix m as well as the labels in v
1069    so that element (i,i) is replaced by the first non-zero element in row i after that element
1070
1071    Throws an exception if there are no non-zero elements in that row
1072    '''
1073    for j in range(i+1,len(v)):
1074        if not np.allclose(m[i,j],0):
1075            m[:,(i,j)] = m[:,(j,i)]
1076            v[i],v[j] = v[j],v[i]
1077            return
1078    else:
1079        raise Exception,'Singular input'
1080
1081def _RowEchelon(m,arr,collist):
1082    '''Convert the first m rows in Matrix arr to row-echelon form
1083    exchanging columns in the matrix and collist as needed.
1084
1085    throws an exception if the matrix is singular because
1086    the first m rows are not linearly independent
1087    '''
1088    n = len(collist)
1089    for i in range(m):
1090        if np.allclose(arr[i,i],0):
1091            _SwapColumns(i,arr,collist)
1092        arr[i,:] /= arr[i,i] # normalize row
1093        # subtract current row from subsequent rows to set values to left of diagonal to 0
1094        for j in range(i+1,m):
1095            arr[j,:] -= arr[i,:] * arr[j,i]
1096
1097if __name__ == "__main__":
1098    parmdict = {}
1099    constrDict = [
1100        {'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},
1101        {'0:0:eA': 0.0},
1102        {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
1103        {'1::C(10,0,1)': 1.0, '2::C(10,0,1)': 1.0},
1104        {'1::AUiso:0': 1.0, '0::AUiso:0': 1.0},
1105        {'0::A0': 0.0}
1106        ]
1107    fixedList = ['5.0', '0', None, None, '1.0', '0']
1108    StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), ))
1109    #StoreEquivalence('1::atomx:3',('2::atomx:3', ('2::atomz:3',2,), )) # error: dependent & independent vars mixed
1110    #StoreEquivalence('1::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) # error: dependent vars repeated
1111    #StoreEquivalence('0:1:eA',('0:0:eA',)) # error: equiv & fixed
1112    #StoreEquivalence('0:99:Scale',('0:12:Scale',)) # error: equiv & constrained
1113    #StoreEquivalence('0:12:Scale',('0:99:Scale',)) # error: equiv & constrained
1114    varylist = ['2::atomx:3',
1115                '2::C(10,6,1)', '1::C(10,6,1)',
1116                '2::atomy:3', '2::atomz:3',
1117                '0:12:Scale', '0:11:Scale', '0:14:Scale', '0:13:Scale', '0:0:Scale']
1118#    e,w = CheckConstraints([,
1119#                     [{'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}],
1120#                     ['1.0'])
1121#    if e: print 'error=',e
1122#    if w: print 'error=',w
1123#    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']
1124#    constrDict = [
1125#        {'0::Afrac:4': 24.0, '0::Afrac:1': 16.0, '0::Afrac:3': 24.0, '0::Afrac:2': 16.0},
1126#        {'0::Afrac:1': 1.0, '0::Afrac:2': 1.0},
1127#        {'0::Afrac:4': 1.0, '0::Afrac:3': 1.0}]
1128#    fixedList = ['40.0', '1.0', '1.0']
1129
1130    errmsg, warnmsg = CheckConstraints(varylist,constrDict,fixedList)
1131    if errmsg:
1132        print "*** Error ********************"
1133        print errmsg
1134    if warnmsg:
1135        print "*** Warning ********************"
1136        print warnmsg
1137    if errmsg or warnmsg:
1138        sys.exit()
1139    groups,parmlist = GroupConstraints(constrDict)
1140    GenerateConstraints(groups,parmlist,varylist,constrDict,fixedList)
1141    print VarRemapShow(varylist)
1142    parmdict.update( {
1143        '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,
1144        '0:0:eA': 0.0,
1145        '2::C(10,6,1)': 0.2, '1::C(10,6,1)': 0.3,
1146        '1::C(10,0,1)': 0.2, '2::C(10,0,1)': 0.3,
1147        '1::AUiso:0': 0.02, '0::AUiso:0': 0.03,
1148        '0::A0': 0.0,
1149        '2::atomx:3':0.23,'2::atomy:3':-.23, '2::atomz:3':-0.11,
1150        })
1151    print 'parmdict start',parmdict
1152    print 'varylist start',varylist
1153    before = parmdict.copy()
1154    Map2Dict(parmdict,varylist)
1155    print 'parmdict before and after Map2Dict'
1156    print '  key / before / after'
1157    for key in sorted(parmdict.keys()):
1158        print '  '+key,'\t',before.get(key),'\t',parmdict[key]
1159    print 'varylist after',varylist
1160    before = parmdict.copy()
1161    Dict2Map(parmdict,varylist)
1162    print 'after Dict2Map'
1163    print '  key / before / after'
1164    for key in sorted(parmdict.keys()):
1165        print '  '+key,'\t',before.get(key),'\t',parmdict[key]
1166#    dMdv = len(varylist)*[0]
1167#    deriv = {}
1168#    for i,v in enumerate(parmdict.keys()): deriv[v]=i
1169#    Dict2Deriv(varylist,deriv,dMdv)
Note: See TracBrowser for help on using the repository browser.