source: trunk/GSASIImapvars.py @ 927

Last change on this file since 927 was 903, checked in by toby, 12 years ago

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