source: trunk/GSASIImapvars.py @ 1046

Last change on this file since 1046 was 1046, checked in by vondreele, 8 years ago

put in missing GSASIIpath & version finder

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