source: trunk/GSASIImapvars.py @ 762

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

update py files to use native line ends, add update from svn to help menu when svn or svn.exe is in path or ./svn/bin

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