Changeset 5118


Ignore:
Timestamp:
Dec 23, 2021 6:00:10 PM (8 months ago)
Author:
toby
Message:

revisions for constraint processing & display

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIconstrGUI.py

    r5082 r5118  
    750750            errmsg,warnmsg = CheckConstraints(G2frame,Phases,Histograms,data,[],reqVaryList,seqhst=seqhistnum,seqmode=seqmode)
    751751        except Exception as msg:
     752            if GSASIIpath.GetConfigValue('debug'):
     753                import traceback
     754                print (traceback.format_exc())
    752755            return 'CheckConstraints error retrieving parameter\nError='+str(msg),''
    753756        return errmsg,warnmsg
    … …  
    12581261            panel.delBtn.checkboxList.append([constDel,Id,name])
    12591262            if refineflag:
    1260                 ch = G2G.G2CheckBox(panel,'vary ',item,-2)
     1263                refresh = lambda event: wx.CallAfter(UpdateConstraints, G2frame, data, G2frame.constr.GetSelection(), True)
     1264                ch = G2G.G2CheckBox(panel,'vary ',item,-2,OnChange=refresh)
    12611265                constSizer.Add(ch,0,wx.LEFT|wx.RIGHT|WACV|wx.ALIGN_CENTER,1)
    12621266            else:
    … …  
    13841388            btn.Bind(wx.EVT_BUTTON,lambda event:
    13851389                         G2G.ShowScrolledInfo(panel,
    1386                         'Constraints after processing\n\n'+G2mv.VarRemapShow(),
     1390                        '*** Constraints after processing ***\n'+G2mv.VarRemapShow(),
    13871391                         header='Generated constraints'))
    13881392            panel.delBtn = wx.Button(panel, wx.ID_ANY, 'Delete selected')

  • trunk/GSASIIdataGUI.py

    r5106 r5118  
    52975297        for key in ('parmMinDict','parmMaxDict','parmFrozen'):
    52985298            if key not in Controls: Controls[key] = {}
     5299        G2mv.Map2Dict(parmValDict,G2mv.saveVaryList)
    52995300        wx.EndBusyCursor()
    53005301        # # check for limits on dependent vars
    … …  
    53185319        #    imp.reload(G2G)
    53195320        # end debug stuff   
    5320         dlg = G2G.ShowLSParms(self,'Least Squares Parameters',parmValDict,varyList,reqVaryList,Controls)
     5321        dlg = G2G.ShowLSParms(self,'Least Squares Parameters',parmValDict,G2mv.saveVaryList,reqVaryList,Controls)
    53215322        dlg.CenterOnParent()
    53225323        dlg.ShowModal()

  • trunk/GSASIImapvars.py

    r5110 r5118  
    507507
    508508This is done by creating a square matrix from the group using
    509 :func:`_FillArray` with parameter ``FillDiagonals=False`` (the default). Any
    510 unspecified rows are left as all zero. The first Nc rows in the
    511 array are then coverted to row-echelon form using :func:`_RowEchelon`. This
    512 will create an Exception if any two rows are linearly dependent (which means
    513 that no matter what values are used for the remaining rows, that the matrix
    514 will be singular). :func:`_FillArray` is then called with parameter
    515 ``FillDiagonals=True``, which again creates a square matrix but where
    516 unspecified rows are zero except for the diagonal elements. The 
     509:func:`_FillArray`. The top Nc rows in the matrix are filled
     510as described above. Then :func:`_RowEchelon` is used to see if
     511those entries in the matrix can be coverted to row-echelon form. This
     512will raise an Exception there is linear dependence between the initial Nc rows
     513(which means that no matter what values are used for any remaining rows, that the matrix
     514will be singular). If that is not the case and Nc<Np then any remaining rows that
     515were not specified are filled in. For each of these rows, first only the
     516diagonal element in that row of the matrix is set to 1
     517and the upper portion of the matrix is again tested with :func:`_RowEchelon`
     518to check for linear independence. This is likely to be non-singular,
     519but should :func:`_RowEchelon` fail,
     520:func:`_FillArray` will then try setting each other element in that row to either
     5211 or -1. One of those options should be linearly independent from every other
     522row of the matrix.
     523
     524The 
    517525`Gram-Schmidt process <http://en.wikipedia.org/wiki/Gram-Schmidt>`_,
    518526implemented  in :func:`GramSchmidtOrtho`, is used to find orthonormal unit
    519 vectors for the remaining Np-Nc rows of the matrix. This will fail with
    520 a ``ConstraintException`` if this is not possible (singular matrix) or
    521 the result is placed in :data:`constrArr` as a numpy array.
     527vectors which are used to replace the remaining Np-Nc rows of the matrix. This will fail with
     528a ``ConstraintException`` if this is not possible (singular matrix), but that would be
     529unexpected since the matrix could be converted to row-echelon form. The
     530Gram-Schmidt result is placed in :data:`constrArr` as a numpy array.
    522531
    523532Rows in the matrix corresponding to "New Var" constraints and those that
    524 were generated by the Gram-Schmidt process are provided with parameter names
    525 (this can be specified if a "New Var" entry by using a ``"_name"`` element
    526 in the constraint dict, but at present this is not implemented.) Names are
    527 generated using :data:`paramPrefix` which is set to ``"::constr"``, plus a
    528 number to make the new parameter name unique.
     533were generated by the Gram-Schmidt process are provided with parameter names.
     534These names are generated using :data:`paramPrefix`, which is set to ``"::constr"``,
     535plus a number to make the new parameter name unique,
     536unless a name was specified for the
     537"New Var" entry by using a ``"_name"`` element in the constraint dict.
    529538
    530539Finally the parameters used as input to the constraint are placed in
    … …  
    609618**Mixing Hold (Fixed) parameters in equivalences**
    610619
    611 * If one parameter is designated as a "Hold" in an equivalence, then all parameters in that
     620* If one parameter (or more) is designated as a "Hold" in an equivalence, then all parameters in that
    612621  equivalence cannot be varied. Considering this equivalence::
    613622
    … …  
    660669These globals are expected to be used only by this module's (:mod:`GSASIImapvars`) internal routines.
    661670
     671Lists with information from Constraint Equation and New Var constraints. Each entry
     672in these variables is related to a group of constraints.
    662673
    663674.. tabularcolumns:: |l|p{4.5in}|
    … …  
    667678=============================  ===================================================================
    668679:data:`dependentParmList`        a list containing group of lists of
    669                                  parameters used in the group. Note that parameters listed in
     680                                 parameters used in the group.
     681                                 The columns of the matrices in :data:`arrayList` match
     682                                 the order of parameters here.
     683                                 Note that parameters listed in
    670684                                 dependentParmList will not be included in the Hessian as their
    671685                                 derivatives will not affect the model
    672686
    673 :data:`indParmList`              a list containing groups of Independent parameters defined in
    674                                  each group. This contains both parameters used in parameter
    675                                  redefinitions as well as names of generated new parameters.
     687:data:`indParmList`              a list containing groups of variables or constants matching
     688                                 the columns of the matrices in :data:`invarrayList`.
    676689
    677690:data:`arrayList`                a list containing group of relationship matrices to relate
    … …  
    681694                                 parameters in indParmList to those in dependentParmList.
    682695                                 Unlikely to be used externally.
    683 
    684 :data:`holdParmList`             a list of parameters that have been marked as "Hold".
    685                                  Unlikely to be used externally. Holds have a single
    686                                  entry in the indParmList parameter list.
    687                                  Set in :func:`StoreHold`. (Also see :data:`newHolds`)
    688696
    689697:data:`symGenList`               a list of boolean values that will be True to indicate
    … …  
    691699                                 meaning it is generated based on symmetry, twining
    692700                                 or Pawley overlap.
     701=============================  ===================================================================
     702
     703Lists with information from Hold and Equivalence constraints. Each entry
     704in these variables is related to a group of constraints.
     705
     706.. tabularcolumns:: |l|p{4.5in}|
     707
     708=============================  ===================================================================
     709  variable                      explanation
     710=============================  ===================================================================
     711
     712:data:`holdParmList`             a list of parameters that have been marked as "Hold".
     713                                 Unlikely to be used externally.
     714                                 
     715                                 Set in :func:`StoreHold`. (Also see :data:`newHolds`)
    693716
    694717:data:`dependentVars`            a list of dependent variables in equivalences, compiled
    … …  
    696719                                 Used within :func:`GetDependentVars`.
    697720
    698 :data:`independentVars`          a list of dependent variables in equivalences, compiled
    699                                  from (:data:`indParmList`).
     721:data:`independentVars`          a list of dependent variables in equivalences.
    700722                                 Used within :func:`GetIndependentVars`.
    701723
    … …  
    772794#------------------------------------------------------------------------------------
    773795# Global vars used for storing constraints and equivalences after processing
    774 #   note that constraints are stored listed by contraint groups,
    775 #   where each constraint
    776 #   group contains those parameters that must be handled together
     796#   note that new var and constraint equations are stored together in groups,
     797#   where each constraint group contains those parameters that must be handled
     798#   together. Equivalences are also stored in these
    777799
    778800dependentParmList = []
    … …  
    781803indParmList = [] # a list of names for the new parameters
    782804'''a list of lists where each item contains a list for each constraint group with
    783 fixed values for constraint equations and names of generated (New Var) parameters.
     805fixed values for constraint equations and names of generated/New Var parameters.
     806In the case of equivalences, the name of a single independent parameter is stored.
    784807'''
    785808arrayList = []
    … …  
    793816(in :data:`dependentParmList`).
    794817'''
     818symGenList = []
     819'''A list of flags that if True indicates a constraint was generated by symmetry
     820'''
    795821holdParmList = []
    796822'''List of parameters that should not be refined ("Hold"s).
    797823Set in :func:`StoreHold`. Initialized in :func:`InitVars`.
    798 '''
    799 symGenList = []
    800 '''A list of flags that if True indicates a constraint was generated by symmetry
    801824'''
    802825dependentVars = []
    … …  
    825848
    826849#------------------------------------------------------------------------------------
    827 # global variables is used in :func:`getConstrError` to store error and warning information.
     850# global variables used in :func:`getConstrError` to store error and warning information.
    828851# set in CheckEquivalences and in GenerateConstraints
    829852convVarList = []
    … …  
    905928    groups = [] # contains a list of grouplists
    906929    ParmList = []
    907     for i,consi in enumerate(constrDict):
     930    for i,constrI in enumerate(constrDict):
    908931        if i in assignedlist: continue # already in a group, skip
    909932        # starting a new group
    910933        grouplist = [i,]
    911934        assignedlist.append(i)
    912         groupset = set(VarKeys(consi))
     935        groupset = set(VarKeys(constrI))
    913936        changes = True # always loop at least once
    914937        while(changes): # loop until we can't find anything to add to the current group
    915938            changes = False # but don't loop again unless we find something
    916             for j,consj in enumerate(constrDict):
     939            for j,constrJ in enumerate(constrDict):
    917940                if j in assignedlist: continue # already in a group, skip
    918                 if len(set(VarKeys(consj)) & groupset) > 0: # true if this needs to be added
     941                if len(set(VarKeys(constrJ)) & groupset) > 0: # true if this needs to be added
    919942                    changes = True
    920943                    grouplist.append(j)
    921944                    assignedlist.append(j)
    922                     groupset = groupset | set(VarKeys(consj))
     945                    groupset = groupset | set(VarKeys(constrJ))
    923946        group = sorted(grouplist)
    924947        varlist = sorted(list(groupset))
    … …  
    931954    :func:`ProcessConstraints` into lists ``constrDict`` and ``fixedList``
    932955
    933     This routine then calls :func:`CheckEquivalences`
    934     for internal consistency. This includes converting equivalenced variables into
     956    This routine then calls :func:`CheckEquivalences` for internal
     957    consistency. This includes converting equivalenced variables into
    935958    constraints when a variable is used in both.
    936959
    … …  
    9941017    global saveVaryList
    9951018    saveVaryList = copy.copy(varyList)
     1019   
     1020    global dependentVars # List of dependent and independent variables for all constraints
     1021    global independentVars
    9961022
    9971023    errmsg = ''  # save error messages here. If non-blank, constraints cannot be used.
    9981024    warning = '' # save informational text messages here.
    9991025
    1000     # tabulate a list of all parameters as defined initially setting the following global variables:
    1001     global depVarList # parameters used in equivalences as dependent parameters
    1002     global indepVarList # parameters used in equivalences as independent parameters
    1003     global constrVarList  # list of parameters used in other constraints
     1026    # Process the equivalences; If there are conflicting parameters, move them into constraints
     1027    warning = CheckEquivalences(constrDict,varyList,fixedList,parmDict)
     1028   
     1029    # find parameters used in constraint equations & new var assignments (all are dependent)
     1030    global constrVarList
    10041031    constrVarList = []
    10051032    for cdict in constrDict:
    10061033        constrVarList += [i for i in cdict if i not in constrVarList and not i.startswith('_')]
    1007     # tabulate all parameters in equivalences by type and look for repeated dependent vars
    1008     global depVarList,indepVarList
    1009     depVarList = []  # list of all dependent parameters in equivalences
    1010     indepVarList = []  # list of all independent parameters in equivalences
    1011     for cnum,(varlist,mapvars,multarr,invmultarr) in enumerate(zip(
    1012             dependentParmList,indParmList,arrayList,invarrayList)):
    1013         if multarr is not None: continue # equivalence
    1014         indepVarList += [mv for mv in mapvars if mv not in indepVarList]
    1015         depVarList += [v for v in varlist if v not in depVarList]
    1016    
    1017     # Translate & generate lookup table for wildcard parameter names (sequential fits only)
    1018     # TODO: where should this be done in the sequence of events? Probably here, so that
    1019     # variable names can be looked up in varyList properly. Also, should names be changed in place
    1020     # or stored in translateTable and then looked up later?     
    1021     #translateTable = {} # lookup table for wildcard-referenced parameters
    1022     if SeqHist is not None:
    1023         for varlist,mapvars,multarr,invmultarr in zip(
    1024             dependentParmList,indParmList,arrayList,invarrayList):
    1025             for i,mv in enumerate(mapvars):
    1026                 if mv.split(':')[1] == '*':
    1027                     # convert wildcard var to reference current histogram; save translation in table
    1028                     sv = mv.split(':')
    1029                     sv[1] = str(SeqHist)
    1030                     #mapvars[i][1] = translateTable[mv] = ':'.join(sv)
    1031                     mapvars[i][1] = ':'.join(sv)
    1032             for i,(v,m) in enumerate(zip(varlist,invmultarr)):
    1033                 if v.split(':')[1] == '*':
    1034                     # convert wildcard var to reference current histogram; save translation in table
    1035                     sv = v.split(':')
    1036                     sv[1] = str(SeqHist)
    1037                     #varlist[i] = translateTable[v] = ':'.join(sv)
    1038                     varlist[i] = ':'.join(sv)
    1039         for rel in constrDict:
    1040             for i,var in enumerate(rel):
    1041                 if var.startswith('_'): continue
    1042                 if var.split(':')[1] == '*': # convert wildcard var to reference current histogram
    1043                     sv = var.split(':')
    1044                     sv[1] = str(SeqHist)
    1045                     #rel[i][1] = translateTable[var] = ':'.join(sv)   # add to translation in table for later use
    1046                     rel[i][1] = ':'.join(sv)
    1047 
    1048     # Process the equivalences; If there are conflicting parameters, move them into constraints
    1049     warning = CheckEquivalences(constrDict,varyList,fixedList,parmDict)
    10501034
    10511035    # look through "Constr" and "New Var" constraints looking for zero multipliers and
    10521036    # Hold, Unvaried & Undefined parameters
    10531037    skipList = []
    1054     for cnum,(cdict,val) in enumerate(zip(constrDict,fixedList)):
    1055         valid = 0          # count of good parameter
     1038    for cnum,(cdict,fixVal) in enumerate(zip(constrDict,fixedList)):
     1039        valid = 0          # count of good parameters
    10561040        # error reporting
    10571041        zeroList = []      # parameters with zero multipliers
    … …  
    10631047        problem = False    # constraint must be dropped
    10641048        dropList = []      # parameters to remove
    1065         setAsHoldList = [] # parameters that will be held if the constraint is invalid
    1066         for i,var in enumerate(cdict):
    1067             if var.startswith('_'): continue
     1049        for var in VarKeys(cdict):   # assemble warning info
    10681050            if cdict[var] == 0:    # zero multiplier
    10691051                if var not in zeroList: zeroList.append(var)
    1070                 setAsHoldList.append(var)
    1071                 dropList.append(var)
    1072             elif var in holdParmList:   # already hold
    1073                 #if var not in newHolds: newHolds.append(var)
     1052                if var not in dropList: dropList.append(var)
     1053            elif var in holdParmList:   # hold invalid in New Var, drop from constraint eqn
    10741054                holdList.append(var)
    1075                 dropList.append(var)
    1076             elif ':*:' in var and SeqHist is None:  # unvaried
     1055                if fixVal is None:
     1056                    problem = True
     1057                else:
     1058                    if var not in dropList: dropList.append(var)
     1059            elif ':*:' in var :  # wildcard still present should be treated as undefined
     1060                if var not in undefinedVars: undefinedVars.append(var)
    10771061                noWildcardList.append(var)
    1078                 dropList.append(var)
    1079             elif var not in varyList:  # unvaried
    1080                 if var not in unvariedParmsList: unvariedParmsList.append(var)
    1081                 noVaryList.append(var)
    1082                 setAsHoldList.append(var)
    1083                 dropList.append(var)
     1062                problem = True
    10841063            elif parmDict is not None and var not in parmDict: # not defined, constraint will not be used
    10851064                if var not in undefinedVars: undefinedVars.append(var)
    10861065                notDefList.append(var)
    1087                 if ':dAx:' in var or ':dAy:' in var or ':dAz:' in var: # undefined atoms
    1088                     dropList.append(var) # these can be ignored
     1066                if ':dAx:' in var or ':dAy:' in var or ':dAz:' in var: # coordinates from undefined atoms
     1067                    if fixVal is None:
     1068                        problem = True  # invalid in New Var
     1069                    else:
     1070                        if var not in dropList: dropList.append(var) # ignore in constraint eqn
    10891071                else:
    10901072                    problem = True
     1073            elif var not in varyList and fixVal is not None:  # unvaried, constraint eq. only
     1074                if var not in unvariedParmsList: unvariedParmsList.append(var)
     1075                noVaryList.append(var)
     1076                dropList.append(var)
    10911077            else:
    10921078                valid += 1
    … …  
    11011087                    if i != 0: msg += ', '
    11021088                    msg += v
    1103                 warn(msg,cdict,val)
    1104         if not valid and (problem or len(dropList) > 0): # no valid entries
    1105             warn('Ignoring constraint',cdict,val)
     1089                warn(msg,cdict,fixVal)
     1090        if valid == 0: # no valid entries
     1091            warn('Ignoring constraint, no valid entries',cdict)
    11061092            skipList.append(cnum)
    11071093        elif problem: # mix of valid & refined and undefined items, cannot use this
    1108             warn('Holding varied items & Dropping constraint',cdict,val)
     1094            warn('Invalid constraint',cdict)
    11091095            skipList.append(cnum)
    1110             newHolds += [i for i in holdList if i not in newHolds]
     1096            newHolds += [i for i in VarKeys(cdict) if i not in newHolds]
    11111097        elif len(dropList) > 0: # mix of valid and problematic items, drop problem vars, but keep rest
    11121098            if GSASIIpath.GetConfigValue('debug'):
    … …  
    11151101                    if msg: msg += ' ,'
    11161102                    msg += v
    1117                 warn('removing: '+msg,cdict,val)
     1103                warn('removing: '+msg,cdict)
    11181104            value = fixedList[cnum]
    11191105            for var in dropList:   # do cleanup
    1120                 # TODO: evaluate expressions in constraint multipliers here? I am assuming no
    1121                 # as I think G2mv.EvaluateMultipliers does that already. For now assuming
    1122                 # otherwise; note SubfromParmDict
    1123                 if ':dAx:' in var or ':dAy:' in var or ':dAz:' in var: # undefined atoms can be ignored
    1124                     pass
     1106                # NB expressions in constraint multipliers have already been evaluated
     1107                if ':dAx:' in var or ':dAy:' in var or ':dAz:' in var:
     1108                    pass # treat delta coords as 0; no change in total is needed
    11251109                elif cdict[var] != 0:
    11261110                    value = float(value) - cdict[var]*parmDict[var]
    11271111                del cdict[var]
    1128             if float(value) != float(fixedList[cnum]): fixedList[cnum] = str(np.round(value,12))
     1112            if float(value) != float(fixedList[cnum]): fixedList[cnum] = float(np.round(value,12))
    11291113            if GSASIIpath.GetConfigValue('debug'):
    11301114                warn('revised as: '+_FormatConstraint(constrDict[cnum],fixedList[cnum]))
    … …  
    11451129    # and all newly created relationships will be varied. For NewVar constraints,
    11461130    # vary if the vary flag was set.
    1147     for group,varlist in zip(groups,parmlist):
    1148         if len(varlist) < len(group): # too many relationships -- no can do
     1131    for group,depPrmList in zip(groups,parmlist):
     1132        if len(depPrmList) < len(group): # too many relationships -- no can do
    11491133            if errmsg: errmsg += '\n'
    11501134            errmsg += "Over-constrained input. "
    11511135            errmsg += "There are more constraints (" + str(len(group))
    1152             errmsg += ") than parameters (" + str(len(varlist)) + ")\nin these constraints:"
     1136            errmsg += ") than parameters (" + str(len(depPrmList)) + ")\nin these constraints:"
    11531137            for rel in group:
    11541138                errmsg += '\n\t'+ _FormatConstraint(constrDict[rel],fixedList[rel])
    1155             groupErrors += varlist
     1139            groupErrors += depPrmList
    11561140            continue # go on to next group
    11571141
    1158         VaryFree = False
    1159         for rel in group:
    1160             varied = 0
    1161             for var in VarKeys(constrDict[rel]):
    1162                 #var = translateTable.get(var,var) # replace wildcards
    1163                 #if parmDict is not None and var not in parmDict:
    1164                 if var in varyList: varied += 1
    1165             if fixedList[rel] is not None and varied > 0:
    1166                 VaryFree = True
    1167                
    1168         # fill in additional degrees of freedom
    11691142        try:
    1170             arr = _FillArray(group,constrDict,varlist)
    1171             _RowEchelon(len(group),arr,varlist)
    1172         except:
     1143            constrArr = _FillArray(group,constrDict,depPrmList)
     1144        except Exception as err:
    11731145            if errmsg: errmsg += '\n'
    1174             errmsg += "\nSingular input. "
    1175             errmsg += "There are internal inconsistencies in these constraints:"
    1176             for rel in group:
    1177                 errmsg += '\n\t' + _FormatConstraint(constrDict[rel],fixedList[rel])
    1178             groupErrors += varlist
     1146            if 'Initial' in str(err):
     1147                errmsg += "\nSingular input. "
     1148                errmsg += "There are internal inconsistencies in these constraints:"
     1149            else:
     1150                errmsg += "\nError expanding matrix with these constraints:"
     1151                for rel in group:
     1152                    errmsg += '\n\t' + _FormatConstraint(constrDict[rel],fixedList[rel])
     1153            groupErrors += depPrmList
    11791154            continue
    11801155
    11811156        try:
    1182             constrArr = _FillArray(group,constrDict,varlist,FillDiagonals=True)
    11831157            GramSchmidtOrtho(constrArr,len(group))
    11841158        except:
    … …  
    11871161            for rel in group:
    11881162                errmsg += '\n\t' + _FormatConstraint(constrDict[rel],fixedList[rel])
    1189             groupErrors += varlist
     1163            groupErrors += depPrmList
    11901164            continue
    11911165       
    1192         try:
    1193             np.linalg.inv(constrArr)
     1166        try: 
     1167            invConstrArr = np.linalg.inv(constrArr)
    11941168        except:
    11951169            if errmsg: errmsg += '\n'
    … …  
    12011175            for rel in group:
    12021176                errmsg += '\n\t' + _FormatConstraint(constrDict[rel],fixedList[rel])
    1203             groupErrors += varlist
     1177            groupErrors += depPrmList
    12041178            continue
    1205        
    1206         mapvar = []
    1207         group = group[:]
    1208         # scan through all generated and input relationships, we need to add to the varied list
    1209         # all the new parameters where VaryFree has been set or where a New Var is varied.
    1210         #
    1211         # If a group does not contain any fixed values (constraint equations)
    1212         # and nothing in the group is varied, drop this group, so that the
    1213         # dependent parameters can be refined individually.
    1214         unused = True
    1215         for i in range(len(varlist)):
    1216             if len(group) > 0: # get the original equation reference
    1217                 rel = group.pop(0)
    1218                 fixedval = fixedList[rel]
    1219                 varyflag = constrDict[rel].get('_vary',False)
    1220                 varname = constrDict[rel].get('_name','')
    1221             else: # this relationship has been generated
    1222                 varyflag = False
    1223                 varname = ''
    1224                 fixedval = None
    1225             if fixedval is None: # this is a new parameter, not a constraint
    1226                 if not varname:
    1227                     varname = paramPrefix + str(consNum) # no assigned name, create one
     1179
     1180        # scan through current group looking for new var assignments
     1181        hasNewVar = False
     1182        for i,rel in enumerate(group):
     1183            if fixedList[rel] is None:
     1184                hasNewVar = True # there a New Var relationship in this group
     1185                break
     1186        else: # only constraint equations, check for unvaried parameters in
     1187            unvaried = False
     1188            # this should not happen as they should have been removed
     1189            for var in depPrmList:
     1190                if var not in varyList:
     1191                    unvaried = True
     1192                    break
     1193            if unvaried: # something is not varied: skip group & remove all parameters from varyList
     1194                for var in depPrmList:
     1195                    if var not in newHolds: newHolds.append(var)
     1196                if GSASIIpath.GetConfigValue('debug'):
     1197                    print('Unexpected: Constraint group ignored (some parameters unvaried)')
     1198                    for rel in group:
     1199                        print ('   '+_FormatConstraint(constrDict[rel],fixedList[rel]))                   
     1200                continue
     1201        maplist = []   # value or param name mapped to each row in expression matrix           
     1202        if not hasNewVar: # constraint equations; all remaining entries varied, vary generated
     1203            for i in range(len(depPrmList)):
     1204                if i >= len(group): # tag generated degrees of freedom with names and vary them
     1205                    varname = paramPrefix + str(consNum) # assign a unique name
    12281206                    consNum += 1
    1229                 mapvar.append(varname)
    1230                 #genVarLookup[varname] = varlist # save list of parameters related to this new var
    1231                 # vary the new relationship if it is a degree of freedom in
    1232                 # a set of contraint equations or if a New Var is flagged to be varied.
    1233                 if VaryFree or varyflag:
    1234                     unused = False
     1207                    maplist.append(varname)
    12351208                    varyList.append(varname)
    1236                     # fix (prevent varying) of all the parameters inside the constraint group
    1237                     # (dependent vars)
    1238                     for var in varlist:
    1239                         if var in varyList: varyList.remove(var)
    1240             else:
    1241                 unused = False
    1242                 mapvar.append(float(fixedval))
    1243         if unused: # skip over constraints that don't matter (w/o fixed value or any refined parameters)
    1244             if GSASIIpath.GetConfigValue('debug'):
    1245                 print('Unexpected: Constraint ignored (all parameters unvaried)')
    1246                 print ('   '+_FormatConstraint(constrDict[rel],fixedList[rel]))
    1247             continue
    1248         #dependentParmList.append([translateTable.get(var,var) for var in varlist])
    1249         dependentParmList.append(varlist)
     1209                else:
     1210                    maplist.append(fixedList[rel])
     1211        else:   # ------------------------- groups with new var assignments, vary only NV's w/flags set
     1212            for i,rel in enumerate(group):
     1213                if fixedList[rel] is None:
     1214                    varname = constrDict[rel].get('_name','::?????')
     1215                    maplist.append(varname)
     1216                    if  constrDict[rel].get('_vary',False):
     1217                        varyList.append(varname)
     1218                else:
     1219                   maplist.append(fixedList[rel])   # constraint equation
     1220            for i in range(len(depPrmList)):
     1221                if i >= len(group): # name generated degrees of freedom
     1222                    varname = paramPrefix + str(consNum) # assign a unique name
     1223                    consNum += 1
     1224                    maplist.append(varname)
     1225            for var in depPrmList:
     1226                if var not in newHolds: newHolds.append(var)
     1227        # keep this group
     1228        dependentParmList.append(depPrmList)
    12501229        arrayList.append(constrArr)
    1251         invarrayList.append(np.linalg.inv(constrArr))
    1252         indParmList.append(mapvar)
     1230        invarrayList.append(invConstrArr)
     1231        indParmList.append(maplist)
    12531232        symGenList.append(False)
     1233           
    12541234    if errmsg and SeqHist is not None:
    12551235        print (' *** ERROR in constraint definitions! ***')
    … …  
    12621242        return errmsg,warning,None,None
    12631243
    1264     # Make list of dependent and independent variables (after possibly dropping unused vars)
    1265     global dependentVars
    1266     global independentVars
     1244    # Make list of dependent and independent variables for all constraints
    12671245    dependentVars = []
    12681246    independentVars = []
    … …  
    12731251        for mv in varlist:
    12741252            if mv not in dependentVars: dependentVars.append(mv)
    1275     saveVaryList = copy.copy(varyList)
     1253            if mv not in newHolds: newHolds.append(mv)
     1254    saveVaryList = copy.copy(varyList)  # save varyList so it can be used within module
    12761255
    12771256    # if equivMoved:
    … …  
    13311310    removeList = []  # equivalences that are not needed
    13321311    convertList = [] # equivalences that should be converted to "Const" constraints
     1312   
     1313    # tabulate parameters used in equivalences by type
     1314    depVarList = []  # list of all dependent parameters in equivalences
     1315    indepVarList = []  # list of all independent parameters in equivalences
     1316    for cnum,(varlist,mapvars,multarr,invmultarr) in enumerate(zip(
     1317            dependentParmList,indParmList,arrayList,invarrayList)):
     1318        #if multarr is not None: continue # equivalence
     1319        indepVarList += [mv for mv in mapvars if mv not in indepVarList]
     1320        depVarList += [v for v in varlist if v not in depVarList]
    13331321
    13341322    # process equivalences: make a list of dependent and independent vars
    … …  
    14201408            else:
    14211409                msg = 'All parameters set as "Hold". '
    1422             msg += " Ignoring equivalence"
     1410            msg += "\n Ignoring equivalence"
    14231411            warn(msg,cnum)
    14241412            removeList.append(cnum)
    … …  
    14461434            else:
    14471435                msg = 'No parameters varied. '               
    1448             msg += " Ignoring equivalence"
     1436            msg += "\n Ignoring equivalence"
    14491437            warn(msg,cnum)
    14501438            removeList.append(cnum)
    … …  
    15451533    :param str seqmode: one of 'use-all', 'wildcards-only' or 'auto-wildcard'.
    15461534       When seqmode=='wildcards-only' then any constraint with a numerical
    1547        constraint number is skipped. With seqmode=='auto-wildcard',
     1535       histogram number is skipped. With seqmode=='auto-wildcard',
    15481536       any non-null constraint number is set to the selected histogram.
    15491537    :param int seqhst: number for current histogram (used for
    … …  
    15641552    namedVarList = []
    15651553    for constr in constList:
    1566         terms = copy.deepcopy(constr[:-3])
     1554        terms = copy.deepcopy(constr[:-3]) # don't change the tree contents
     1555        # deal with wildcards in sequential fits
    15671556        if seqmode == 'wildcards-only' and seqhst is not None:
    15681557            skip = False
    … …  
    15761565            for term in terms:
    15771566                term[1] = term[1].varname(seqhst)
    1578         else:
     1567        elif seqhst is not None:
    15791568            for term in terms:
    1580                 if term[1].histogram == '*' and seqhst is not None:
     1569                if term[1].histogram == '*':
    15811570                    term[1] = term[1].varname(seqhst)
    1582                 else:
    1583                     term[1] = term[1].varname()
    1584                
     1571                # else:
     1572                #     term[1] = term[1].varname()   # does this change anything???
     1573        # separate processing by constraint type
    15851574        if constr[-1] == 'h':
    15861575            # process a hold
    1587             var = str(constr[0][1])
     1576            var = str(term[0][1])
    15881577            if '?' not in var:
    15891578                StoreHold(var)
    … …  
    16021591            if len(D) > 1:
    16031592                # add extra dict terms for input variable name and vary flag
    1604                 if varname is not None:
     1593                if varname is None: # no assigned name, create one
     1594                    global consNum
     1595                    varname = str(consNum)
     1596                    consNum += 1
     1597                else:
    16051598                    varname = str(varname) # in case this is a G2VarObj
    1606                     if varname.startswith(':'):
    1607                         D['_name'] = varname
    1608                     else:
    1609                         D['_name'] = '::nv-' + varname
    1610                     D['_name'] = G2obj.MakeUniqueLabel(D['_name'],namedVarList)
     1599                if '::' in varname:
     1600                    D['_name'] = varname.replace('::','::nv-')
     1601                else:
     1602                    D['_name'] = '::nv-' + varname
     1603                D['_name'] = G2obj.MakeUniqueLabel(D['_name'],namedVarList)
    16111604                D['_vary'] = varyFlag == True # force to bool
    16121605                constrDict.append(D)
    … …  
    16151608            #constFlag[-1] = ['Vary']
    16161609        elif constr[-1] == 'c':
    1617             # process a contraint relationship
     1610            # process a constraint equation
    16181611            D = {}
    16191612            for term in terms:
    … …  
    16221615                    D[var] = term[0]
    16231616            if len(D) >= 1:
    1624                 fixedList.append(str(constr[-3]))
     1617                fixedList.append(float(constr[-3]))
    16251618                constrDict.append(D)
    16261619            else:
    … …  
    20272020            elif v in undefinedVars:
    20282021                undef.append(str(v))
    2029             elif v in unvariedParmsList:
     2022            elif v in unvariedParmsList and not (constrLst[-1] == 'f' and constrLst[-2]):
    20302023                unvar.append(str(v))
    20312024            elif v in holdParmList:
    … …  
    21722165    s = ''
    21732166    if len(holdParmList) > 0:
    2174         s += 'User-supplied Fixed Parameters:\n'
     2167        s += '\nUser-supplied Fixed Parameters:\n'
    21752168        for v in holdParmList:
    21762169            s += '    ' + v + '\n'
    21772170    if len(newHolds) > 0:
    2178         s += 'Additional Fixed Parameters:\n'
     2171        s1 = ''
     2172        s2 = ''
    21792173        for v in newHolds:
    2180             s += '    ' + v + '\n'
     2174            if v in dependentVars:
     2175                s1 += '    ' + v + '\n'
     2176            else:
     2177                s2 += '    ' + v + '\n'
     2178        if s2:
     2179            s += '\nParameters to be Fixed based on constraints:\n'
     2180            s += s2
     2181        if s1:
     2182            s += '\nDependent parameters with values to be set from constraints:\n'
     2183            s += s1
    21812184    userOut = ''
    21822185    symOut = ''
    21832186    consOut = ''
    21842187    varOut = ''
     2188    freeOut = ''
    21852189    global dependentParmList,arrayList,invarrayList,indParmList,symGenList
    21862190
    … …  
    22082212                    userOut += s1 + '\n'
    22092213                continue
    2210             if mv in varyList:
    2211                 lineOut = '  {} (**Varied**) = '.format(mv)
    2212             else:
    2213                 lineOut = '  {} = '.format(mv)
     2214            ln = ''
     2215            # if mv in varyList:
     2216            #     lineOut = '  (V)* {} = '.format(mv)
     2217            # else:
     2218            #     lineOut = '  {} = '.format(mv)
    22142219            j = 0
     2220            lineOut = '  {} = '.format(mv)
    22152221            for (m,v) in zip(multarr[i,:],varlist):
    22162222                if m == 0: continue
    … …  
    22212227                    lineOut += ' + '
    22222228                j += 1
    2223                 if len(lineOut) > 60 and type(mv) is float:
    2224                     consOut += lineOut
    2225                     lineOut = '\n\t'
    2226                 elif len(lineOut) > 60:
    2227                     varOut += lineOut
     2229                if len(lineOut) > 60:
     2230                    ln += lineOut
    22282231                    lineOut = '\n\t'
    22292232                if m == 1:
    … …  
    22312234                else:
    22322235                    lineOut += '({:.4g} * {})'.format(m,v)
    2233             if type(mv) is not float:
    2234                 varOut += lineOut + '\n'
     2236            if mv in varyList:
     2237                lineOut += ' **VARIED**'
     2238            ln += lineOut
     2239           
     2240            if type(mv) is float:
     2241                consOut += ln + '\n'
     2242            elif '::nv-' in mv:
     2243                varOut += ln + '\n'
    22352244            else:
    2236                 consOut += lineOut + '\n'
     2245                freeOut += ln + '\n'
    22372246    if userOut:
    2238         s += '\nUser-supplied equivalences:\n' + userOut
     2247        s += '\nEquivalences:\n' + userOut
    22392248    if consOut:
    2240         s += '\nUser-supplied Constraints:\n' + consOut
     2249        s += '\nConstraint Equations:\n' + consOut
    22412250    if varOut:
    2242         s += '\nUser-input generated New Variables:\n' + varOut
     2251        s += '\nNew Variable assignments:\n' + varOut
     2252    if freeOut:
     2253        s += '\nGenerated Degrees of Freedom:\n' + freeOut
    22432254    if symOut:
    22442255        s += '\nSymmetry-generated equivalences:\n' + symOut
    … …  
    23712382                    if debug: print ('start dMdv',dMdv[varyList.index(name)])
    23722383                    if debug: print ('add derv',v,'*',m,'to derv',name,'add=',m * derivDict[v])
    2373                     if m == 0: continue
     2384                    if m == 0 or name not in varyList: continue
    23742385                    dMdv[varyList.index(name)] += m * derivDict[v]
    23752386
    … …  
    24022413    global dependentParmList,arrayList,invarrayList,indParmList
    24032414    for varlist,mapvars,multarr,invmultarr in zip(dependentParmList,indParmList,arrayList,invarrayList):
    2404         for item in varlist:
     2415        for item in varlist: # TODO: is this still needed?
    24052416            if item in varyList: varyList.remove(item)
    24062417        if multarr is None:
    … …  
    24172428    for varlist,mapvars,multarr in zip(dependentParmList,indParmList,arrayList):
    24182429        if multarr is None: continue
    2419         # evaluate all constraints in the forward direction
    2420         for key,value in zip(mapvars,np.dot(multarr,np.array([parmDict[var] for var in varlist]))):
    2421             if type(key) is float: continue
    2422             #print('parmDict set',key,':',value)
    2423             parmDict[key] = value
     2430        # evaluate constraints in the forward direction
     2431        z = zip(mapvars,np.dot(multarr,np.array([parmDict[var] for var in varlist])))
     2432        # add/replace in parameter dict
     2433        parmDict.update([i for i in z if type(i[0]) is not float])
    24242434    global saveVaryList
    24252435    saveVaryList = copy.copy(varyList)
    … …  
    24722482    return a
    24732483
    2474 def _FillArray(sel,dict,collist,FillDiagonals=False):
    2475     '''Construct a n by n matrix (n = len(collist)
    2476     filling in the rows using the relationships defined
    2477     in the dictionaries selected by sel
    2478 
    2479     If FillDiagonals is True, diagonal elements in the
    2480     array are set to 1.0
     2484def _FillArray(sel,d,collist):
     2485    '''Construct a n by n matrix [n = len(collist)]
     2486    with the initial m rows [m = len(sel)] using the
     2487    relationships defined in the expressions dict, d.
     2488    Since m may be smaller than n, the remaining rows
     2489    are filled with rows that are tested to not create
     2490    a singular matrix.
     2491
     2492    :param list sel:     a list of indices in dict d
     2493    :param list d:       a list of dict's where each dict describes an
     2494      expression from a constraint equation or a new var
     2495    :param list collist: a list parameter names.
     2496    :returns: an n by n numpy.array matrix
    24812497    '''
    24822498    n = len(collist)
    2483     if FillDiagonals:
    2484         arr = np.eye(n)
    2485     else:
    2486         arr = np.zeros(2*[n,])
     2499    m = len(sel)
     2500    arr = np.zeros(2*[n,])
    24872501    # fill the top rows
    24882502    for i,cnum in enumerate(sel):
    24892503        for j,var in enumerate(collist):
    2490             arr[i,j] = dict[cnum].get(var,0)
     2504            arr[i,j] = d[cnum].get(var,0)
     2505    try:
     2506        _RowEchelon(m,copy.copy(arr),collist)
     2507    except:
     2508        raise Exception('Initial constraints singular')
     2509    for i in range(m,n):
     2510        arr[i][i] = 1   # add a diagonal element
     2511        try:
     2512            _RowEchelon(i+1,copy.copy(arr),collist)
     2513            continue
     2514        except:
     2515            pass
     2516        for j in range(n):
     2517            if j == i: continue
     2518            arr[i][j] = 1   # add another element
     2519            try:
     2520                _RowEchelon(i+1,copy.copy(arr),collist)
     2521                break
     2522            except:
     2523                pass
     2524            arr[i][j] = -1   # try a different valuefor this element
     2525            try:
     2526                _RowEchelon(i+1,copy.copy(arr),collist)
     2527                break
     2528            except:
     2529                arr[i][j] = 0   # reset to add another element
     2530        else:
     2531            raise Exception('Unable to create non-singular matrix')
    24912532    return arr
    24922533
    … …  
    25212562
    25222563if __name__ == "__main__":
     2564    d = [{'a': 0, 'b': 1.5,'d':0}, {'d': -1}]
     2565    lbls = ['a','c','b','d']
     2566    sel = (0,1)
     2567    try:
     2568        arr2 = _FillArray(sel,d,lbls)
     2569    except Exception as err:
     2570        if 'Initial' in str(err):
     2571            print('initial error')
     2572        else:
     2573            print('unable to extend matrix error')
     2574        import sys; sys.exit()
     2575    print(arr2)
     2576
     2577    d = [{'a': 1, 'b': 1,}]
     2578    lbls = ['a','b']
     2579    sel = (0,)
     2580    arr1 = _FillArray(sel,d,lbls)
     2581    print(arr1)
     2582    print(GramSchmidtOrtho(arr1,1))
     2583   
     2584    import sys; sys.exit()
     2585   
    25232586    parmdict = {}
    25242587    constrDict = [

  • trunk/GSASIIobj.py

    r5110 r5118  
    19911991        'RBR([TLS])([123AB][123AB])' : 'Residue rigid body group disp. param.',
    19921992        'RBRU' : 'Residue rigid body group Uiso param.',
    1993         'constr([0-9]*)' : 'Parameter from constraint',
    1994         'nv-([^_]+)_*' : 'New variable constraint parameter named \\1',
     1993        'constr([0-9]*)' : 'Generated degree of freedom from constraint',
     1994        'nv-(.+)' : 'New variable assignment with name \\1',
    19951995        # supersymmetry parameters  p::<var>:a:o 'Flen','Fcent'?
    19961996        'mV([0-2])$' : 'Modulation vector component \\1',
Note: See TracChangeset for help on using the changeset viewer.