Changeset 3802


Ignore:
Timestamp:
Jan 27, 2019 4:30:39 PM (3 years ago)
Author:
toby
Message:

generalized restraints completed

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIobj.py

    r3711 r3802  
    26832683    def UpdateDict(self,parmDict):
    26842684        '''Update the dict for the expression with values in a dict
    2685         :param list parmDict: a dict of values some of which may be in use here
     2685        :param dict parmDict: a dict of values, items not in use are ignored
    26862686        '''
    26872687        if self.eObj.expression.startswith('Dist') or self.eObj.expression.startswith('Angle'):
  • trunk/GSASIIrestrGUI.py

    r3780 r3802  
    18791879        mainSizer.Add((5,5),0)
    18801880        hSizer = wx.BoxSizer(wx.HORIZONTAL)
    1881         hSizer.Add(wx.StaticText(GeneralRestr,wx.ID_ANY,'(not implemented yet)'))
    1882         hSizer.Add((5,5),0)
    18831881        hSizer.Add(wx.StaticText(GeneralRestr,wx.ID_ANY,'Weight factor: '))
    18841882        hSizer.Add(
     
    19071905                    0,wx.ALIGN_CENTER)
    19081906            for i,rest in enumerate(generalRestData['General']):
    1909                 txt = rest[0].expression
     1907                eq = rest[0]
     1908                txt = eq.expression
    19101909                if len(txt) > 50:
    19111910                    txt = txt[:47]+'... '
     
    19181917                    calcobj = G2obj.ExpressionCalcObj(rest[0])
    19191918                    calcobj.SetupCalc(parmDict)
    1920                     txt = '{:f}'.format(calcobj.EvalExpression())
     1919                    txt = ' {:f} '.format(calcobj.EvalExpression())
    19211920                except:
    1922                     txt = '(error)'
     1921                    txt = ' (error) '
    19231922                GridSiz.Add(wx.StaticText(GeneralRestr,wx.ID_ANY,txt))
    19241923                GridSiz.Add(
     
    19331932                btn.Bind(wx.EVT_BUTTON,OnDelGenRestraint)
    19341933                GridSiz.Add(btn)
    1935                 txt = str(rest[0].assgnVars)[1:-1].replace("'","")
     1934                txt = ''
     1935                for i in eq.assgnVars:
     1936                    if txt: txt += '; '
     1937                    txt += str(i) + '=' + str(eq.assgnVars[i])
    19361938                if len(txt) > 50:
    19371939                    txt = txt[:47]+'...'
     
    21412143    Pages.append(txt)
    21422144   
    2143     if GSASIIpath.GetConfigValue('debug'):
    2144         txt = 'General'
    2145         GeneralRestr = wx.ScrolledWindow(G2frame.restrBook)
    2146         G2frame.restrBook.AddPage(GeneralRestr,txt)
    2147         Pages.append(txt)
     2145    txt = 'General'
     2146    GeneralRestr = wx.ScrolledWindow(G2frame.restrBook)
     2147    G2frame.restrBook.AddPage(GeneralRestr,txt)
     2148    Pages.append(txt)
    21482149   
    21492150    if General['SH Texture']['Order']:
  • trunk/GSASIIstrMath.py

    r3801 r3802  
    303303
    304304def penaltyFxn(HistoPhases,calcControls,parmDict,varyList):
    305     'Needs a doc string'
     305    'Compute user-supplied and built-in restraint functions'
    306306    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
    307307    pNames = []
     
    327327        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
    328328            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
    329             ['ChemComp','Sites'],['Texture','HKLs'],]
     329            ['ChemComp','Sites'],['Texture','HKLs'],['General','General'],]
    330330        for name,rest in names:
    331331            pWsum[name] = 0.
     
    335335            itemRest = phaseRest[name]
    336336            if itemRest[rest] and itemRest['Use']:
    337                 wt = itemRest['wtFactor']
     337                wt = itemRest.get('wtFactor',1.)
    338338                if name in ['Bond','Angle','Plane','Chiral']:
    339339                    for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
     
    405405                                pWsum[name] += wt*(Z2/esd2)**2
    406406                                pWnum[name] += 1
     407                elif name == 'General':
     408                    for i,(eq,obs,esd) in enumerate(itemRest[rest]):
     409                        pNames.append(str(pId)+':'+name+':'+str(i))
     410                        calcobj = G2obj.ExpressionCalcObj(eq)
     411                        calcobj.SetupCalc(parmDict)
     412                        calc = calcobj.EvalExpression()
     413                        pVals.append(obs-calc)
     414                        pWt.append(wt/esd**2)                   
     415                        pWsum[name] += wt*((obs-calc)/esd)**2
     416                        pWnum[name] += 1
    407417       
    408418    for phase in Phases:
     
    458468   
    459469def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList):
    460     'Needs a doc string'
     470    '''Compute derivatives on user-supplied and built-in restraint
     471    (penalty) functions
     472
     473    where pNames is list of restraint labels
     474
     475    returns pDerv with partial derivatives by variable# in varList and
     476       restraint# in pNames (pDerv[variable#][restraint#])
     477    '''
    461478    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
    462479    pDerv = np.zeros((len(varyList),len(pVal)))
     480    for pName in pNames: # loop over restraints
     481        if 'General' == pName.split(':')[1]:
     482            # initialize for General restraint(s) here
     483            GeneralInit = True
     484            parmDict0 = parmDict.copy()
     485            # setup steps for each parameter
     486            stepDict = {}
     487            for parm in varyList:
     488                if parm.split(':')[2].startswith('dA'):
     489                    stepDict[parm] = 1e-5
     490                elif True:
     491                    stepDict[parm] = 1e-4
     492            break
    463493    for phase in Phases:
    464494#        if phase not in restraintDict:
     
    484514            'ChemComp':'Sites','Texture':'HKLs'}
    485515        lasthkl = np.array([0,0,0])
    486         for ip,pName in enumerate(pNames):
     516        for ip,pName in enumerate(pNames): # loop over restraints
    487517            pnames = pName.split(':')
    488518            if pId == int(pnames[0]):
     
    546576                            dNames += [str(pId)+'::'+SHname]
    547577                            deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
     578                elif name == 'General':
     579                    deriv = []
     580                    dNames = []
     581                    eq,obs,esd = itemRest[name][id]
     582                    calcobj = G2obj.ExpressionCalcObj(eq)
     583                    parmlist = list(eq.assgnVars.values()) # parameters used in this expression
     584                    for parm in parmlist: # expand list if any parms are determined by constraints
     585                        if parm in G2mv.dependentVars:
     586                            parmlist += G2mv.independentVars
     587                            break
     588                    for ind,var in enumerate(varyList):
     589                        drv = 0
     590                        if var in parmlist:
     591                            step = stepDict.get(var,1e-5)
     592                            calc = []
     593                            # apply step to parameter
     594                            oneparm = True
     595                            for s in -step,2*step:
     596                                parmDict[var] += s
     597                                # extend shift if needed to other parameters
     598                                if var in G2mv.independentVars:
     599                                    G2mv.Dict2Map(parmDict,[])
     600                                    oneparm = False
     601                                elif var in G2mv.dependentVars:
     602                                    G2mv.Map2Dict(parmDict,[])
     603                                    oneparm = False
     604                                if 'RB' in var:
     605                                    ApplyRBModels(parmDict,Phases,rigidbodyDict)
     606# test
     607                                    oneparm = False
     608                                calcobj.SetupCalc(parmDict)
     609                                calc.append(calcobj.EvalExpression())
     610                            drv = (calc[1]-calc[0])*.5/step
     611                            # restore the dict
     612                            if oneparm:
     613                                parmDict[var] = parmDict0[var]
     614                            else:
     615                                parmDict = parmDict0.copy()
     616                        else:
     617                            drv = 0
     618                        pDerv[ind][ip] = drv
     619                # Add derivatives into matrix, if needed
    548620                for dName,drv in zip(dNames,deriv):
    549621                    try:
     
    821893    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
    822894        Flack = 1.-2.*parmDict[phfx+'Flack']
    823     time0 = time.time()
     895#    time0 = time.time()
    824896#reflection processing begins here - big arrays!
    825897    iBeg = 0
     
    14731545        for iel,El in enumerate(refDict['FF']['El']):
    14741546            refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
    1475     time0 = time.time()
     1547#    time0 = time.time()
    14761548#reflection processing begins here - big arrays!
    14771549    iBeg = 0
     
    16251697            for iel,El in enumerate(refDict['FF']['El']):
    16261698                refDict['FF']['FF'].T[iel] = G2el.ScatFac(FFtables[El],SQ)
    1627     time0 = time.time()
     1699#    time0 = time.time()
    16281700#reflection processing begins here - big arrays!
    16291701    iBeg = 0
Note: See TracChangeset for help on using the changeset viewer.