Changeset 567


Ignore:
Timestamp:
Apr 20, 2012 3:34:58 PM (11 years ago)
Author:
toby
Message:

more work on constraints

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASII.py

    r564 r567  
    859859        '''Window to edit Constraint values
    860860        '''
    861         def __init__(self,parent,title,text,data):
     861        def __init__(self,parent,title,text,data,separator='*'):
    862862            wx.Dialog.__init__(self,parent,-1,title,
    863863                pos=wx.DefaultPosition,style=wx.DEFAULT_DIALOG_STYLE)
    … …  
    871871            dataGridSizer = wx.FlexGridSizer(rows=len(data),cols=2,hgap=2,vgap=2)
    872872            for id,item in enumerate(self.data[:-1]):
    873                 #name = wx.TextCtrl(panel,-1,item[1],size=wx.Size(200,20))
    874                 name = wx.StaticText(panel,-1,item[1],size=wx.Size(200,20))
    875                 #name.SetEditable(False)
     873                lbl = item[1]
     874                if lbl[-1] != '=': lbl += ' ' + separator + ' '
     875                name = wx.StaticText(panel,-1,lbl,size=wx.Size(200,20),
     876                                     style=wx.ALIGN_RIGHT)
    876877                scale = wx.TextCtrl(panel,id,'%.3f'%(item[0]),style=wx.TE_PROCESS_ENTER)
    877878                scale.Bind(wx.EVT_TEXT_ENTER,self.OnScaleChange)
    878879                scale.Bind(wx.EVT_KILL_FOCUS,self.OnScaleChange)
    879                 dataGridSizer.Add(scale,0,wx.LEFT,10)
    880                 dataGridSizer.Add(name,0,wx.RIGHT,10)
     880                dataGridSizer.Add(name,0,wx.LEFT,10)
     881                dataGridSizer.Add(scale,0,wx.RIGHT,10)
    881882            mainSizer.Add(dataGridSizer,0,wx.EXPAND)
    882883            OkBtn = wx.Button(panel,-1,"Ok")
    … …  
    895896            panel.Fit()
    896897            self.Fit()
     898            self.CenterOnParent()
    897899           
    898900        def OnNameChange(self,event):

  • trunk/GSASIIgrid.py

    r549 r567  
    12281228        '''This creates a sizer displaying all of the constraints entered
    12291229        '''
    1230         constSizer = wx.FlexGridSizer(1,5,0,0)
     1230        constSizer = wx.FlexGridSizer(1,4,0,0)
     1231        maxlen = 70 # characters before wrapping a constraint
    12311232        for Id,item in enumerate(data[name]):
    12321233            eqString = ['',]
    … …  
    12421243                if item[-1] == 'f':
    12431244                    for term in item[:-3]:
    1244                         if len(eqString[-1]) > 60:
     1245                        if len(eqString[-1]) > maxlen:
    12451246                            eqString.append(' ')
     1247                        m = term[0]
    12461248                        if eqString[-1] != '':
    1247                             if term[0] > 0:
     1249                            if m >= 0:
    12481250                                eqString[-1] += ' + '
    12491251                            else:
    12501252                                eqString[-1] += ' - '
    1251                         eqString[-1] += '%.3f*%s '%(abs(term[0]),term[1])
     1253                                m = abs(m)
     1254                        eqString[-1] += '%.3f*%s '%(m,term[1])
    12521255                    typeString = ' NEWVAR  '
    12531256                    eqString[-1] += ' = New Variable   '
    12541257                elif item[-1] == 'c':
    12551258                    for term in item[:-3]:
    1256                         if len(eqString[-1]) > 60:
     1259                        if len(eqString[-1]) > maxlen:
    12571260                            eqString.append(' ')
    12581261                        if eqString[-1] != '':
    … …  
    12651268                    eqString[-1] += ' = %.3f'%(item[-3])+'  '
    12661269                elif item[-1] == 'e':
    1267                     first = 1.0
    12681270                    for term in item[:-3]:
    1269                         if len(eqString[-1]) > 60:
     1271                        if term[0] == 0: term[0] = 1.0
     1272                        if len(eqString[-1]) > maxlen:
    12701273                            eqString.append(' ')
    12711274                        if eqString[-1] == '':
    12721275                            eqString[-1] += '%s '%(term[1])
    1273                             if term[0] != 0: first = term[0]
     1276                            first = term[0]
    12741277                        else:
    1275                             eqString[-1] += ' = %.3f*%s '%(term[0]/first,term[1])
     1278                            eqString[-1] += ' = %.3f*%s '%(first/term[0],term[1])
    12761279                    typeString = ' EQUIV   '
    12771280                else:
    … …  
    12901293                EqSizer.Add(wx.StaticText(pageDisplay,-1,s),0,wx.ALIGN_CENTER_VERTICAL)
    12911294            constSizer.Add(EqSizer,0,wx.ALIGN_CENTER_VERTICAL)
    1292             if item[-1] == 'f':
    1293                 constRef = wx.CheckBox(pageDisplay,-1,label=' Refine?')
    1294                 constRef.SetValue(item[-2])
    1295                 constRef.Bind(wx.EVT_CHECKBOX,OnConstRef)
    1296                 Indx[constRef.GetId()] = item
    1297                 constSizer.Add(constRef)
    1298             else:
    1299                 constSizer.Add((5,5),0)
     1295            # if item[-1] == 'f':
     1296            #     constRef = wx.CheckBox(pageDisplay,-1,label=' Refine?')
     1297            #     constRef.SetValue(item[-2])
     1298            #     constRef.Bind(wx.EVT_CHECKBOX,OnConstRef)
     1299            #     Indx[constRef.GetId()] = item
     1300            #     constSizer.Add(constRef)
     1301            # else:
     1302            #     constSizer.Add((5,5),0)
    13001303        return constSizer
    13011304               
    1302     def OnConstRef(event):
    1303         Obj = event.GetEventObject()
    1304         Indx[Obj.GetId()][-2] = Obj.GetValue()
     1305    # def OnConstRef(event):
     1306    #     Obj = event.GetEventObject()
     1307    #     Indx[Obj.GetId()][-2] = Obj.GetValue()
    13051308       
    13061309    def OnConstDel(event):
    … …  
    13141317        Obj = event.GetEventObject()
    13151318        Id,name = Indx[Obj.GetId()]
     1319        sep = '*'
    13161320        if data[name][Id][-1] == 'f':
    13171321            items = data[name][Id][:-3]+[[],]
    … …  
    13201324        elif data[name][Id][-1] == 'c':
    13211325            items = data[name][Id][:-3]+[
    1322                 [data[name][Id][-3],'= fixed value'],[]]
     1326                [data[name][Id][-3],'fixed value ='],[]]
    13231327            constType = 'Constraint'
    13241328            lbl = 'Edit value for each term in constant constraint sum'
    … …  
    13261330            items = data[name][Id][:-3]+[[],]
    13271331            constType = 'Equivalence'
    1328             lbl = 'Variable ' + items[0][1] + ' is equal to: '
     1332            lbl = 'The following terms are set to be equal:'
     1333            sep = '/'
    13291334        else:
    13301335            return
    1331         dlg = G2frame.ConstraintDialog(G2frame,constType,lbl,items)
     1336        dlg = G2frame.ConstraintDialog(G2frame.dataFrame,constType,lbl,items,sep)
    13321337        try:
    13331338            if dlg.ShowModal() == wx.ID_OK:

  • trunk/GSASIImapvars.py

    r484 r567  
    99and parameter simplification contraints.
    1010
    11 Parameter redefinition is done by creating one or more relationships
     11Parameter redefinition (new vars) is done by creating one or more relationships
    1212between a set of parameters
    1313   Mx1 * Px + My1 * Py +...
    … …  
    1515where Pj is a parameter name and Mjk is a constant.
    1616
    17 Relations are typically supplied as input to InputParse, where either
    18 of two keywords can be attached to a relationship:
    19  * VARY which means the generated parameter name will be included in
    20   the list of variables to be refined (varyList)
    21  * VARYFREE which means that all generated parameter names for a
    22   group, including the ones for generated relationships (see below)
    23   will be included in the list of variables to be refined (varyList)
    24   The case of VARY and VARYFREE is immaterial.
    25 
    26 Relations can also be supplied in the form of an equation:
     17Constant constraint Relations can also be supplied in the form of an equation:
    2718  nx1 * Px + ny1 * Py +... = C1
    2819where Cn is a constant. These equations define an algebraic
    29 constrant. The keyword VARY makes no sense when used with a constraint
    30 equation, but VARYFREE can be used.  #RVD??? is VARYFREE required or default??
    31 
    32 Unconstrained relations describe a new, independent, parameter, which
     20constrant.
     21
     22Parameters can also be "fixed" (held), which prevents them from being refined.
     23
     24All of the above three cases is supplied as input using routines
     25GroupConstraints and GenerateConstraints. The input consists of a list of
     26relationship dictionaries:
     27    constrDict = [
     28         {'0:12:Scale': 2.0, '0:14:Scale': 4.0, '0:13:Scale': 3.0, '0:0:Scale': 0.5},
     29         {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
     30         {'0::A0': 0.0}]
     31
     32    fixedList = ['5.0', None, '0']
     33
     34Where the dictionary defines the first part of an expression and the corresponding fixedList
     35item is either None (for parameter redefinition) or the constant value, for a constant
     36constraint equation. A dictionary that contains a single term defines a variable that
     37will be fixed (held). The multiplier and the fixedList value in this case are ignored.
     38
     39Parameters can also be equivalenced or "slaved" to another parameter, such that one
     40(independent) parameter is equated to several (now dependent) parameters. In
     41algebraic form this is:
     42   P0 = M1 * P1 = M2 * P2 = ...
     43Thus parameters P0, P1 and P2,... are linearly equivalent. Routine StoreEquivalence is
     44used to specify these equivalences.
     45
     46Parameter redefinition (new vars) describes a new, independent, parameter, which
    3347is defined in terms of dependent parameters that are defined in the
    34 Model, while constrained relations effectively reduce the complexity
    35 of the Model by removing a degree of freedom.
     48Model, while fixed constrained relations effectively reduce the complexity
     49of the Model by removing a degree of freedom. It is possible for a parameter to appear
     50in both a parameter redefinition expression and a fixed constraint equation, but a
     51parameter cannot be used a parameter equivalance cannot be used elsewhere (not fixed,
     52constrained or redefined). Likewise a fixed parameter cannot be used elsewhere (not
     53equivalanced, constrained or redefined).
    3654
    3755Relationships are grouped so that a set of dependent parameters appear
    … …  
    4260the output from InputParse and GroupConstraints generates the
    4361"missing" relationships and saves that information in the module's
    44 global variables. Each generated parameter is named
    45 sequentially using paramPrefix.
    46 
    47 Parameter redefinition is done by equating one (independent) parameter to several
    48 (now dependent) parameters, in algebraic form:
    49   P1 = n2 * P2 = n3 * P3 ...
    50 Each equality in the relationship reduces the complexity of the model
    51 by one potentially variable parameter. Input is provided to routine
    52 StoreEquivalence in the form of an independent parameter and a list of
    53 dependent parameters, optionally with a multiplier.
    54 
    55 Note that none of the dependent parameters in any constraint or
    56 reformulation can be refined (see dependentParmList, below).
     62global variables. Each generated parameter is named sequentially using paramPrefix.
     63
     64A list of parameters that will be varied is specified as input to GenerateConstraints
     65(varyList). A fixed parameter will simply be removed from this list preventing that
     66parameter from being varied. Note that all parameters in a relationship must specified as
     67varied (appear in varyList) or none can be varied. This is checked in GenerateConstraints
     68(as well as for generated relationships in SetVaryFlags).
     69* If all parameters in a parameter redefinition (new var) relationship are varied, the
     70  parameter assigned to this expression (::constr:n, see paramPrefix) newly generated
     71  parameter is varied. Note that any generated "missing" relations are not varied. Only
     72  the input relations are varied.
     73* If all parameters in a fixed constraint equation are varied, the generated "missing"
     74  relations in the group are all varied. This provides the N-C degrees of freedom.
    5775
    5876External Routines:
    … …  
    109127
    110128import numpy as np
    111 import re
    112129# data used for constraints;
    113130debug = False # turns on printing as constraint input is processed
    … …  
    123140fixedVarList = [] # List of variables that should not be refined
    124141
    125 # compile regular expressions used for parsing input
    126 rex_mult = re.compile('[+-]?[0-9.]+[eE]?[+-]?[0-9]*')
    127 rex_star = re.compile('\s*\*\s*')
    128 rex_var = re.compile('[a-zA-Z0-9_:]+')
    129 rex_plusminus = re.compile('\s*[+-=]\s*')
    130 rex_vary = re.compile('\s*Vary\s*', re.IGNORECASE)
    131 rex_varyfree = re.compile('(.*)\s*VaryFree\s*', re.IGNORECASE)
    132 
    133142# prefix for parameter names
    134143paramPrefix = "::constr:"
    … …  
    145154    consNum = 0 # number of the next constraint to be created
    146155    fixedVarList = []
    147 
    148 def InputParse(mapList):
    149     '''Converts a set relationships used to remap parameters into composite
    150     parameters or to one or more constants.
    151        
    152     Input:
    153       mapList: a list of character strings where each string defines
    154     a relationship in form:
    155     ('<const11>*<prm1> [+-] <const12>*<prm2> [+-] ... = <value> [VaryFree]',
    156     '<const21>*<prm1> [+-] <const22>*<prm2> [+-] ... [Vary/VaryFree]',
    157     ...)
    158     these define either a constraint or a new independent parameter,
    159     where <constXX> is a constant containing an optional sign, digits,
    160     optionally a decimal point and/or an exponent, prefixed by e or E,
    161     and <prmN> is the name of a parameter defined in the Model.
    162     Parameters can be included in any order. Multiplying a parameter
    163     by 0 causes that parameter to be included in a group of
    164     relationships.  Note that if the same parameter is included twice
    165     in a relationship/equation, the coefficients are added.
    166    
    167     When an equality is defined, a constant, <value>, is
    168     included. This then describes a constraint equation.
    169    
    170     Vary is an optional keyword that indicates the indicated remapped
    171     parameter should be varied (note, case insensitive). Should be the
    172     last item on the line. Makes no sense with an equation.
    173 
    174     VaryFree is an optional keyword that indicates all possible
    175     remapped parameters should be varied (note, case
    176     insensitive). Makes most sense with a constraint equation.  Should
    177     be the last item on the line.
    178 
    179     returns
    180       constrDict: a list with a dict for each item in mapList that
    181         defines the relationship. The keys are parameter names and the
    182         values are the multiplier for the parameter name
    183       constrFlag: a list for each item in mapList with a list contains
    184         'Vary' and/or 'VaryFree'
    185       fixedList: a list for each item in mapList. Contains the value
    186         (as a string) for each contraint equation or None for a
    187         constraint relationship.
    188     '''
    189     i = 0
    190     constrDict = []
    191     fixedList = []
    192     constrFlag = []
    193     if debug: print 50*'-','\n(Input)'
    194     for line in mapList:
    195         inputline = line[:]
    196         line = line.strip()
    197         if len(line) == 0: continue # empty lines are ignored
    198         constrDict.append({})
    199         constrFlag.append([])
    200         fixedList.append(None)
    201         i += 1
    202         fixedval = None
    203         if debug: print '\t',line
    204         try:
    205             j = 0
    206             sign = 1.0
    207             while len(line):
    208                 j += 1
    209                 m = sign * float(rex_mult.match(line).group())
    210                 line = line[rex_mult.match(line).end():]
    211                 j += 1
    212                 line = line[rex_star.match(line).end():]
    213                 j += 1
    214                 v = rex_var.match(line).group()
    215                 line = line[rex_var.match(line).end():]
    216                 #print m,'times',v
    217                 #if v not in varlist: varlist.append(v)
    218                 if constrDict[i-1].get(v) == None:
    219                     constrDict[i-1][v] = m
    220                 else:
    221                     constrDict[i-1][v] += m
    222                 if len(line.strip()) == 0: break
    223                 j += 1
    224                 # advance to next separator (+-=)
    225                 pm = rex_plusminus.match(line)
    226                 if pm != None:
    227                     line = line[pm.end():]
    228                     pm = pm.group()
    229                 else:
    230                     pm = ''
    231                 if pm.strip() == '+':
    232                     sign = 1.0
    233                 elif pm.strip() == '-':
    234                     sign = -1.0
    235                 elif pm.strip() == '=':
    236                     # found a fixed value, also check for a VaryFree flag
    237                     if rex_varyfree.match(line):
    238                         constrFlag[-1].append('VaryFree')
    239                         #fixedval = float(rex_varyfree.split(line)[1])
    240                         fixedval = rex_varyfree.split(line)[1].strip()
    241                     else:
    242                         #fixedval = float(line.strip())
    243                         fixedval = line.strip()
    244                     fixedList[-1] = fixedval
    245                     line = ""
    246                 elif rex_varyfree.match(line):
    247                     constrFlag[-1].append('VaryFree')
    248                     line = line[rex_varyfree.match(line).end():]
    249                 elif rex_vary.match(line):
    250                     constrFlag[-1].append('Vary')
    251                     line = line[rex_vary.match(line).end():]
    252                 else:
    253                     # something else is on the line, but not a keyword
    254                     raise Exception,'Error in line '+str(inputline)+'\nat '+str(line)
    255         except SyntaxError:
    256             if debug: print 'Error in line',i,'token',j,'@','"'+line+'"'
    257             raise Exception,'Error in line %d token %d, beginning with %s'% (i,j, line)
    258 
    259     if debug: # on debug, show what is parsed in a semi-readable
    260         print 50*'-','\n(parsed relationship/equation & flag)'
    261         for i in range(len(constrDict)):
    262             flags = ''
    263             for f in constrFlag[i]:
    264                 if flags != '':
    265                     flags = flags + ' + ' + f
    266                 else:
    267                     flags = f
    268             if fixedList[i] is None:
    269                 print '#',i+1,constrDict[i],'\t',constrFlag[i]
    270             else:
    271                 print '#',i+1,constrDict[i],'=',fixedList[i],'\t',constrFlag[i]
    272 
    273     return constrDict,constrFlag,fixedList
    274156
    275157def GroupConstraints(constrDict):
    … …  
    307189    return groups,ParmList
    308190
    309 def GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList):
    310     '''Takes a list of relationship entries comprising a group of constraints and
    311     builds the relationship lists and their inverse and stores them in global variables
     191def GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList):
     192    '''Takes a list of relationship entries comprising a group of
     193    constraints and builds the relationship lists and their inverse
     194    and stores them in global variables Also checks for internal
     195    conflicts or inconsistencies in parameter/variable definitions.
    312196    '''
    313197    global dependentParmList,arrayList,invarrayList,indParmList,consNum
     198    msg = ''
     199
     200    # process fixed (held) variables
     201    for cdict in constrDict:
     202        if len(cdict) == 1:
     203            fixedVarList.append(cdict.keys()[0])
     204    #print 'fixedVarList',fixedVarList
     205   
     206    # process equivalences: make a list of dependent and independent vars
     207    #    and check for repeated uses (repetition of a parameter as an
     208    #    independent var is OK)
     209    indepVarList = []
     210    depVarList = []
     211    multdepVarList = []
     212    for varlist,mapvars,multarr,invmultarr in zip(
     213        dependentParmList,indParmList,arrayList,invarrayList):
     214        if multarr is None:
     215            zeromult = False
     216            for mv in mapvars:
     217                varied = 0
     218                notvaried = ''
     219                if mv in varyList:
     220                    varied += 1
     221                else:
     222                    if notvaried: notvaried += ', '
     223                    notvaried += mv
     224                if mv not in indepVarList: indepVarList.append(mv)
     225                for v,m in zip(varlist,invmultarr):
     226                    if m == 0: zeromult = True
     227                    if v in varyList:
     228                        varied += 1
     229                    else:
     230                        if notvaried: notvaried += ', '
     231                        notvaried += v
     232                    if v in depVarList:
     233                        multdepVarList.append(v)
     234                    else:
     235                        depVarList.append(v)
     236            if varied > 0 and varied != len(varlist)+1:
     237                msg += "\nNot all variables refined in equivalence:\n\t"
     238                s = ""
     239                for v in varlist:
     240                    if s != "": s+= " & "
     241                    s += str(v)           
     242                msg += str(mv) + " => " + s
     243                msg += '\nNot refined: ' + notvaried + '\n'
     244            if zeromult:
     245                msg += "\nZero multiplier is invalid in equivalence:\n\t"
     246                s = ""
     247                for v in varlist:
     248                    if s != "": s+= " & "
     249                    s += str(v)           
     250                msg += str(mv) + " => " + s + '\n'
     251
     252    #print 'indepVarList',indepVarList
     253    #print 'depVarList',depVarList
     254    # check for errors:
     255    inboth = set(indepVarList).intersection(set(depVarList))
     256    if len(inboth) > 0:
     257        msg += "\nThe following parameters(s) are used as both dependent and independent variables in Equivalence relations:\n"
     258        s = ''
     259        for var in inboth:
     260            if s != "": s+= ", "
     261            s += str(v)
     262        msg += '\t'+ s + '\n'
     263    if len(multdepVarList) > 0:
     264        msg += "\nThe following parameters(s) are used in multiple Equivalence relations as dependent variables:\n"
     265        s = ''
     266        for var in multdepVarList:
     267            if s != "": s+= ", "
     268            s += str(v)           
     269        msg += '\t'+ s + '\n'
     270    equivVarList = list(set(indepVarList).union(set(depVarList)))
     271    #print 'equivVarList',equivVarList
     272    inboth = set(fixedVarList).intersection(set(equivVarList))
     273    if len(inboth) > 0:
     274        msg += "\nError! The following variables are used in both Equivalence and Fixed constraints:\n"
     275        s = ''
     276        for var in inboth:
     277            if s != "": s+= ", "
     278            s += str(v)
     279        msg += '\t'+ s + '\n'
     280
     281    # scan through parameters in each relationship. Are all varied? If only some are
     282    # varied, create an error message.
     283    # If all are varied and this is a constraint equation, then set VaryFree flag
     284    # so that newly created relationships will be varied
    314285    for group,varlist in zip(groups,parmlist):
     286        if len(varlist) == 1: continue
    315287        VaryFree = False
    316         for row in group:
    317             if 'VaryFree' in constrFlag[row]: VaryFree = True
    318         if len(varlist) == 1:
    319             #print "Debug: Fixing parameter (%s)" % (varlist[0])
    320             fixedVarList.append(varlist[0])
    321             continue
     288        for rel in group:
     289            varied = 0
     290            notvaried = ''
     291            for var in constrDict[rel]:
     292                if var in varyList:
     293                    varied += 1
     294                else:
     295                    if notvaried: notvaried += ', '
     296                    notvaried += var
     297                if var in fixedVarList:
     298                    msg += '\nError: parameter '+var+" is Fixed and used in a constraint:\n\t"
     299                    msg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
     300                if var in equivVarList:
     301                    msg += '\nError: parameter '+var+" is Equivalenced and used in a constraint:\n\t"
     302                    msg += FormatConstraint(constrDict[rel],fixedList[rel])+"\n"
     303            if varied > 0 and varied != len(constrDict[rel]):
     304                msg += "\nNot all variables refined in constraint:\n\t"
     305                msg += FormatConstraint(constrDict[rel],fixedList[rel])
     306                msg += '\nNot refined: ' + notvaried + '\n'
     307            if fixedList[rel] is not None and varied > 0:
     308                VaryFree = True
     309    # if there were errors found, go no farther
     310    if msg:
     311        print ' *** ERROR in constraint definitions! ***'
     312        print msg
     313        raise Exception
     314               
     315    # now process each group and create the relations that are needed to form
     316    # non-singular square matrix
     317    for group,varlist in zip(groups,parmlist):
     318        if len(varlist) == 1: continue
    322319        arr = MakeArray(constrDict, group, varlist)
    323320        constrArr = FillInMissingRelations(arr,len(group))
    324321        mapvar = []
    325322        group = group[:]
     323        # scan through all generated and input variables, add to the varied list
     324        # all the new parameters where VaryFree has been set or where all the
     325        # dependent parameters are varied. Check again for inconsistent variable use
     326        # for new variables -- where varied and unvaried parameters get grouped
     327        # together. I don't think this can happen when not flagged before, but
     328        # it does not hurt to check again.
    326329        for i in range(len(varlist)):
    327             vary = VaryFree
     330            varied = 0
     331            notvaried = ''
    328332            if len(group) > 0:
    329333                rel = group.pop(0)
    330334                fixedval = fixedList[rel]
    331                 if 'Vary' in constrFlag[rel]: vary = True
     335                for var in constrDict[rel]:
     336                    if var in varyList:
     337                        varied += 1
     338                    else:
     339                        if notvaried: notvaried += ', '
     340                        notvaried += var
    332341            else:
    333342                fixedval = None
    … …  
    336345                mapvar.append(varname)
    337346                consNum += 1
    338                 if vary: varyList.append(varname)
     347                if VaryFree or varied > 0:
     348                    varyList.append(varname)
    339349            else:
    340350                mapvar.append(fixedval)
     351            if varied > 0 and notvaried != '':
     352                msg += "\nNot all variables refined in generated constraint\n\t"
     353                msg += '\nNot refined: ' + notvaried + '\n'
    341354        dependentParmList.append(varlist)
    342355        arrayList.append(constrArr)
    343356        invarrayList.append(np.linalg.inv(constrArr))
    344357        indParmList.append(mapvar)
     358    if msg:
     359        print ' *** ERROR in constraint definitions! ***'
     360        print msg
     361        print VarRemapShow(varyList)
     362        raise Exception
    345363    # setup dictionary containing the fixed values
    346364    global fixedDict
    … …  
    395413    return
    396414
    397 def SetVaryFlags(varyList):
    398     '''Adds independent variables to the varyList provided that all
    399     dependent variables are being varied.
    400     Ignores independent variables where no dependent variables are
    401     being varied.
    402     Returns a non-empty text message where some but not all dependent
    403     variables are being varied.
    404     '''
    405     global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
    406     msg = ""
    407     for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
    408         for mv in mapvars:
    409             varied = []
    410             notvaried = []
    411             if mv in fixedDict: continue
    412             if multarr is None:
    413                 if mv in varyList:
    414                     varied.append(mv)
    415                 else:
    416                     notvaried.append(mv)
    417             for v in varlist:
    418                 if v in varyList:
    419                     varied.append(v)
    420                 else:
    421                     notvaried.append(v)
    422             if len(varied) > 0 and len(notvaried) > 0:
    423                 if msg != "": msg += "\n"
    424                 msg += "Error: inconsistent use of parameter " + mv
    425                 msg += "\n  varied: "
    426                 for v in varied: msg += v
    427                 msg += "\n  not varied: "
    428                 for v in notvaried: msg += v   
    429             elif len(varied) > 0 and multarr is not None:
    430                 varyList.append(mv)
    431     return msg
     415# def SetVaryFlags(varyList,fixedList):
     416#     '''Adds independent variables to the varyList provided that all
     417#     dependent variables are being varied.
     418#     Ignores independent variables where no dependent variables are
     419#     being varied.
     420#     Returns a non-empty text message where some but not all dependent
     421#     variables are being varied.
     422#     '''
     423#     global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
     424#     msg = ""
     425#     for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
     426#         for mv in mapvars:
     427#             varied = []
     428#             notvaried = []
     429#             if mv in fixedDict: continue
     430#             if multarr is None:
     431#                 if mv in varyList:
     432#                     varied.append(mv)
     433#                 else:
     434#                     notvaried.append(mv)
     435#             for v in varlist:
     436#                 if v in varyList:
     437#                     varied.append(v)
     438#                 else:
     439#                     notvaried.append(v)
     440#             if len(varied) > 0 and len(notvaried) > 0:
     441#                 if msg != "": msg += "\n"
     442#                 msg += "Error: inconsistent use of parameter " + mv
     443#                 msg += "\n  varied: "
     444#                 for v in varied: msg += v
     445#                 msg += "\n  not varied: "
     446#                 for v in notvaried: msg += v   
     447#             #elif len(varied) > 0 and multarr is not None:
     448#             #    varyList.append(mv)
     449#     return msg
    432450
    433451def GetDependentVars():
    … …  
    489507    returns a dictionary containing the esd values for dependent parameters
    490508    '''
     509    # I think this fails for equivalencies: the ESD for the master (independent variable)
     510    # needs to be adjusted.
    491511    sigmaDict = {}
    492512    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
    493         if invmultarr is None: continue # probably not needed
     513        #if invmultarr is None: continue # probably not needed
    494514        valuelist = [parmDict[var] for var in mapvars]
    495515        # get the v-covar matrix for independent parameters
    … …  
    506526            sigmaDict[v] = np.sqrt(np.inner(vec.T,np.inner(vcov,vec)))
    507527    return sigmaDict
    508        
     528
     529def FormatConstraint(RelDict,RelVal):
     530    '''Formats a Constraint or Function for use in a convenient way'''
     531    s = [""]
     532    for var,val in RelDict.items():
     533        if len(s[-1]) > 60: s.append(' ')
     534        m = val
     535        if s[-1] != "" and m >= 0:
     536            s[-1] += ' + '
     537        elif s[-1] != "":
     538            s[-1] += ' - '
     539            m = abs(m)
     540        s[-1] += '%.3f*%s '%(m,var)
     541    if len(s[-1]) > 60: s.append(' ')
     542    if RelVal is None:
     543        s[-1] += ' = New variable'
     544    else:
     545        s[-1] += ' = ' + RelVal
     546    s1 = ''
     547    for s2 in s:
     548        if s1 != '': s1 += '\n\t'
     549        s1 += s2
     550    return s1
     551
    509552def VarRemapShow(varyList):
    510553    '''List out the saved relationships.
    … …  
    522565        for i,mv in enumerate(mapvars):
    523566            if multarr is None:
    524                 s += '  ' + str(mv) + ' defines parameter(s): '
     567                s += '  ' + str(mv) + ' is equivalent to parameter(s): '
    525568                j = 0
    526569                for v,m in zip(varlist,invmultarr):
    527                     print v,m[0]
     570                    #print v,m[0]
    528571                    if j > 0: s += '  & '
    529572                    j += 1
    … …  
    560603    '''
    561604    global dependentParmList,arrayList,invarrayList,indParmList,invarrayList
    562     for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
     605    for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
    563606        for i,name in enumerate(mapvars):
    564607            # grouped variables: need to add in the derv. w/r
    565608            # dependent variables to the independent ones
    566             if name not in varyList: continue # if independent var not varied
    567             for m,v in zip(invmultarr[:,i],varlist):
    568 #                print 'add derv',v,'*',m,'to derv',name
    569                 if m == 0: continue
    570                 dMdv[varyList.index(name)] += m * derivDict[v]
     609            if name not in varyList: continue # skip if independent var not varied
     610            if multarr is None:
     611                for v,m in zip(varlist,invmultarr):
     612                    #print 'add derv',v,'/',m[0],'to derv',name
     613                    if m == 0: continue
     614                    dMdv[varyList.index(name)] += derivDict[v] / m[0]
     615            else:
     616                for v,m in zip(varlist,multarr[i,:]):
     617                    #print 'add derv',v,'*',m,'to derv',name
     618                    if m == 0: continue
     619                    dMdv[varyList.index(name)] += m * derivDict[v]
    571620
    572621def Map2Dict(parmDict,varyList):
    … …  
    604653
    605654def Dict2Map(parmDict,varyList):
    606     #I think this needs fixing to update parmDict with new values
    607     #   from the constraints based on what is in varyList - RVD
    608655    '''Convert the remapped values back to the original parameters
    609656   
    610657    This should be done to apply constraints to parameter values (use
    611     Map2Dict first). It should be done as part of the Model function
     658    Map2Dict once first). It should be done as part of the Model function
    612659    evaluation, before any computation is done
    613660    '''
     661    #I think this needs fixing to update parmDict with new values
     662    #   from the constraints based on what is in varyList - RVD. Don't think so -- BHT
    614663    global dependentParmList,arrayList,invarrayList,indParmList,fixedDict
    615664    # reset fixed values (should not be needed, but very quick)
    616665    # - this seems to update parmDict with {'0':0.0} & {'1.0':1.0} - probably not what was intended
     666    # not needed, but also not a problem - BHT
    617667    parmDict.update(fixedDict)
    618668    for varlist,mapvars,invmultarr in zip(dependentParmList,indParmList,invarrayList):
    619         if invmultarr is None: continue
     669        #if invmultarr is None: continue
    620670        valuelist = [parmDict[var] for var in mapvars]
    621671        parmDict.update(zip(varlist,
    … …  
    786836if __name__ == "__main__":
    787837    import sys
    788     #remap = MapParameters() # create an object (perhaps better as a module)
    789     #remap.debug = 1
    790     #debug = 1
    791838    parmdict = {}
     839    constrDict = [
     840        {'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},
     841        {'0:0:eA': 0.0},
     842        {'2::C(10,6,1)': 1.0, '1::C(10,6,1)': 1.0},
     843        {'1::C(10,0,1)': 1.0, '2::C(10,0,1)': 1.0},
     844        {'1::AUiso:0': 1.0, '0::AUiso:0': 1.0},
     845        {'0::A0': 0.0}
     846        ]
     847    fixedList = ['5.0', '0', None, None, '1.0', '0']
    792848    StoreEquivalence('2::atomx:3',('2::atomy:3', ('2::atomz:3',2,), ))
    793     varylist = ['2::atomx:3',]
    794     #print VarRemapShow(varylist)
    795     #msg = SetVaryFlags(varylist)
    796     #if msg != "": print msg
    797     varylist = ['2::atomx:3', '2::atomy:3', '2::atomz:3',]
     849    #StoreEquivalence('1::atomx:3',('2::atomx:3', ('2::atomz:3',2,), )) # error: dependent & independent vars mixed
     850    #StoreEquivalence('1::atomx:3',('2::atomy:3', ('2::atomz:3',2,), )) # error: dependent vars repeated
     851    #StoreEquivalence('0:1:eA',('0:0:eA',)) # error: equiv & fixed
     852    #StoreEquivalence('0:99:Scale',('0:12:Scale',)) # error: equiv & constrained
     853    #StoreEquivalence('0:12:Scale',('0:99:Scale',)) # error: equiv & constrained
     854    varylist = ['2::atomx:3',
     855                '2::C(10,6,1)', '1::C(10,6,1)',
     856                '2::atomy:3', '2::atomz:3',
     857                '0:12:Scale', '0:11:Scale', '0:14:Scale', '0:13:Scale', '0:0:Scale']
     858    groups,parmlist = GroupConstraints(constrDict)
     859    GenerateConstraints(groups,parmlist,varylist,constrDict,fixedList)
    798860    print VarRemapShow(varylist)
    799     msg = SetVaryFlags(varylist)
    800     if msg != "": print msg
    801     #parmdict = {
    802     #    '2::atomx:3':.1 ,
    803     #    '2::atomy:3':.2 , # conflicts with constraint
    804     #    '2::atomz:3':.3 ,
    805     #}
    806 
    807     varylist = []
    808     mapList1 = [
    809         "1 * 2::atomz:3 = 0",
    810         "1*p1 + 2e0*p2 - 1.0*p3",
    811         "1*p4 + 1*p7",
    812         "1*p1+2*p2-3.0*p2 VARY",
    813         ]
    814     mapList = [
    815         "1 * 2::atomz:3 = 0",
    816         "1*p1 + 2e0*p2 - 1.0*p3",
    817         "1*p4 + 1*p7",
    818         "1*p1+2*p2-3.0*p2 VARY",
    819         "1*p21 + 0*p22 - 0*p23 + 0*p24 varyfree",
    820         "1*p5 + 2*p6 = 1 varyfree",
    821         "-1 * p6 + 1*p5",
    822         "-10e-1 * p1 - -2*p2 + 3.0*p4",
    823         ]
    824     constrDict,constrFlag,fixedList = InputParse(mapList1)
    825     #print constrDict
    826     #print constrFlag
    827     #print fixedList
    828     groups,parmlist = GroupConstraints(constrDict)
    829     GenerateConstraints(groups,parmlist,varylist,constrDict,constrFlag,fixedList)
    830     print VarRemapShow(varylist)
    831     sys.exit()
    832     parmdict.update( {'p1':1,'p2':2,'p3':3,'p4':4,
    833                       'p6':6,'p5':5,  # conflicts with constraint
    834                       'p7':7,
    835                       "p21":.1,"p22":.2,"p23":.3,"p24":.4,
    836                       })
    837     #print 'parmdict start',parmdict
     861    parmdict.update( {
     862        '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,
     863        '0:0:eA': 0.0,
     864        '2::C(10,6,1)': 0.2, '1::C(10,6,1)': 0.3,
     865        '1::C(10,0,1)': 0.2, '2::C(10,0,1)': 0.3,
     866        '1::AUiso:0': 0.02, '0::AUiso:0': 0.03,
     867        '0::A0': 0.0,
     868        '2::atomx:3':0.23,'2::atomy:3':-.23, '2::atomz:3':-0.11,
     869        })
     870    print 'parmdict start',parmdict
     871    print 'varylist start',varylist
    838872    before = parmdict.copy()
    839     Map2Dict(parmdict,[])
     873    Map2Dict(parmdict,varylist)
    840874    print 'parmdict before and after Map2Dict'
    841875    print '  key / before / after'
    842876    for key in sorted(parmdict.keys()):
    843877        print '  '+key,'\t',before.get(key),'\t',parmdict[key]
    844 
     878    print 'varylist after',varylist
    845879    before = parmdict.copy()
    846     Dict2Map(parmdict,[])
     880    Dict2Map(parmdict,varylist)
    847881    print 'after Dict2Map'
    848882    print '  key / before / after'
    849883    for key in sorted(parmdict.keys()):
    850884        print '  '+key,'\t',before.get(key),'\t',parmdict[key]
     885    dMdv = len(varylist)*[0]
     886    deriv = {}
     887    for i,v in enumerate(parmdict.keys()): deriv[v]=i
     888    Dict2Deriv(varylist,deriv,dMdv)

  • trunk/GSASIIstruct.py

    r553 r567  
    5555    return Controls
    5656   
    57 def GetConstraints(GPXfile,varyList,parmDict):
     57def GetConstraints(GPXfile):
    5858    '''Read the constraints from the GPX file and interpret them
    59     Note that this needs an update to match the new-style constraints
    6059    '''
    6160    constList = []
    … …  
    7372    fl.close()
    7473    constDict = []
    75     constFlag = []
    7674    fixedList = []
     75    ignored = 0
    7776    for item in constList:
    7877        if item[-1] == 'h':
    7978            # process a hold
    8079            fixedList.append('0')
    81             constFlag.append([])
    8280            constDict.append({item[0][1]:0.0})
    8381        elif item[-1] == 'f':
    8482            # process a new variable
    8583            fixedList.append(None)
    86             constFlag.append([])
    8784            constDict.append({})
    8885            for term in item[:-3]:
    8986                constDict[-1][term[1]] = term[0]
    90             if item[-2]:
    91                 constFlag[-1] = ['Vary']
     87            #constFlag[-1] = ['Vary']
    9288        elif item[-1] == 'c':
    9389            # process a contraint relationship
    94             const = item[-3]
    95             vars = 0
    96             used = {}
     90            fixedList.append(str(item[-3]))
     91            constDict.append({})
    9792            for term in item[:-3]:
    98                 if term[1] in varyList:
    99                     used[term[1]] = term[0]
    100                     vars += 1
    101                 else:
    102                     const -= parmDict[term[1]]*term[0]
    103             if vars == 0:
    104                 print 'Contraint is not used',item[:-3],'=',item[-3]
    105                 print 'no refined variables'
    106                 continue
    107             fixedList.append(str(const))
    108             constFlag.append(['VaryFree'])
    109             constDict.append(used)
     93                constDict[-1][term[1]] = term[0]
     94            #constFlag[-1] = ['VaryFree']
    11095        elif item[-1] == 'e':
    11196            # process an equivalence
    112             first = None
     97            firstmult = None
    11398            eqlist = []
    11499            for term in item[:-3]:
    115                 if first is None:
    116                     first = term[0]
    117                     if first == 0: first = 1.0
    118                     firstvar = term[1]
     100                if term[0] == 0: term[0] = 1.0
     101                if firstmult is None:
     102                    firstmult,firstvar = term
    119103                else:
    120                     eqlist.append([term[1],term[0]/first])
    121             G2mv.StoreEquivalence(firstvar,eqlist)               
     104                    eqlist.append([term[1],firstmult/term[0]])
     105            G2mv.StoreEquivalence(firstvar,eqlist)
    122106        else:
    123             pass
    124     return constDict,constFlag,fixedList
     107            ignored += 1
     108    if ignored:
     109        print ignored,'old-style Constraints were rejected'
     110    return constDict,fixedList
    125111   
    126112def GetPhaseNames(GPXfile):
    … …  
    28482834    calcControls = {}
    28492835    calcControls.update(Controls)           
     2836    constrDict,fixedList = GetConstraints(GPXfile)
    28502837    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
    28512838    if not Phases:
    … …  
    28712858    GetFprime(calcControls,Histograms)
    28722859    # do constraint processing
    2873     constrDict,constrFlag,fixedList = GetConstraints(GPXfile,varyList,parmDict)
    28742860    try:
    28752861        groups,parmlist = G2mv.GroupConstraints(constrDict)
    2876         G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
     2862        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
    28772863    except:
    28782864        print ' *** ERROR - your constraints are internally inconsistent ***'
     2865        # traceback for debug
     2866        #import traceback
     2867        #print traceback.format_exc()
    28792868        raise Exception(' *** Refine aborted ***')
    2880     # check to see which generated parameters are fully varied
    2881     msg = G2mv.SetVaryFlags(varyList)
    2882     if msg:
    2883         print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
    2884         print msg
    2885         raise Exception(' *** Refine aborted ***')
     2869    # # check to see which generated parameters are fully varied
     2870    # msg = G2mv.SetVaryFlags(varyList)
     2871    # if msg:
     2872    #     print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
     2873    #     print msg
     2874    #     raise Exception(' *** Refine aborted ***')
    28862875    print G2mv.VarRemapShow(varyList)
    28872876    G2mv.Map2Dict(parmDict,varyList)
    … …  
    29912980    Controls = GetControls(GPXfile)
    29922981    ShowControls(Controls)           
     2982    constrDict,fixedList = GetConstraints(GPXfile)
    29932983    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
    29942984    if not Phases:
    … …  
    30533043        GetFprime(calcControls,Histo)
    30543044        # do constraint processing
    3055         constrDict,constrFlag,fixedList = GetConstraints(GPXfile,varyList,parmDict)
    3056         groups,parmlist = G2mv.GroupConstraints(constrDict)
    3057         G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
     3045        try:
     3046            groups,parmlist = G2mv.GroupConstraints(constrDict)
     3047            G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
     3048        except:
     3049            print ' *** ERROR - your constraints are internally inconsistent ***'
     3050            raise Exception(' *** Refine aborted ***')
     3051        # check to see which generated parameters are fully varied
     3052        # msg = G2mv.SetVaryFlags(varyList)
     3053        # if msg:
     3054        #     print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
     3055        #     print msg
     3056        #     raise Exception(' *** Refine aborted ***')
     3057        #print G2mv.VarRemapShow(varyList)
    30583058        G2mv.Map2Dict(parmDict,varyList)
    30593059        Rvals = {}
Note: See TracChangeset for help on using the changeset viewer.