source: trunk/GSASIImapvars.py @ 1138

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

major constraints revision

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