Changeset 2614


Ignore:
Timestamp:
Jan 6, 2017 1:04:03 PM (5 years ago)
Author:
toby
Message:

add PDF correction optimizer; fix vdWaals bug after adding element in structure plot

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIgrid.py

    r2611 r2614  
    150150
    151151[ wxID_PDFCOPYCONTROLS, wxID_PDFSAVECONTROLS, wxID_PDFLOADCONTROLS,
    152     wxID_PDFCOMPUTE, wxID_PDFCOMPUTEALL, wxID_PDFADDELEMENT, wxID_PDFDELELEMENT,
    153 ] = [wx.NewId() for item in range(7)]
     152    wxID_PDFCOMPUTE, wxID_PDFCOMPUTEALL, wxID_PDFADDELEMENT, wxID_PDFDELELEMENT, wxID_PDFOPT,
     153] = [wx.NewId() for item in range(8)]
    154154
    155155[ wxID_MCRON,wxID_MCRLIST,wxID_MCRSAVE,wxID_MCRPLAY,
     
    21772177        self.PDFEdit.Append(help='Compute all PDFs', id=wxID_PDFCOMPUTEALL, kind=wx.ITEM_NORMAL,
    21782178            text='Compute all PDFs')
     2179        self.PDFEdit.Append(help='Optimize PDF', id=wxID_PDFOPT, kind=wx.ITEM_NORMAL,
     2180            text='Optimize corrections for r<Rmin section of current G(r)')
    21792181        self.PostfillDataMenu()
    21802182       
  • trunk/GSASIIplot.py

    r2612 r2614  
    60306030#            print caller,generalData['Name']
    60316031# end of useful debug
     6032        vdWRadii = generalData['vdWRadii']
    60326033        mapData = generalData['Map']
    60336034        D4mapData = generalData.get('4DmapData',{})
  • trunk/GSASIIpwd.py

    r2600 r2614  
    280280           
    281281def CalcPDF(data,inst,limits,xydata):
    282     'needs a doc string'
     282    '''Computes I(Q), S(Q) & G(r) from Sample, Bkg, etc. diffraction patterns loaded into
     283    dict xydata; results are placed in xydata.
     284    Calculation parameters are found in dicts data and inst and list limits.
     285    The return value is at present an empty list.
     286    '''
    283287    auxPlot = []
    284288    import copy
  • trunk/GSASIIpwdGUI.py

    r2612 r2614  
    46944694                       
    46954695    def OnResetQ(event):
    4696         resetQ.SetValue(False)
    46974696        data['QScaleLim'][1] = qLimits[1]
    46984697        SQmax.SetValue(data['QScaleLim'][1])
     
    48054804        wx.CallAfter(UpdatePDFGrid,G2frame,data)
    48064805        #UpdatePDFGrid(G2frame,data)
     4806
     4807    def OnOptimizePDF(event):
     4808        '''Optimize Flat Bkg, BackRatio & Ruland corrections to remove spurious
     4809        "intensity" from portion of G(r) with r<Rmin.
     4810        Invoked by Optimize PDF button and from menu command.
     4811        '''
     4812        import scipy.optimize as opt
     4813        wx.BeginBusyCursor()
     4814        Min,Init,Done = SetupPDFEval()
     4815        xstart = Init()
     4816        rms = Min(xstart)
     4817        print('Optimizing corrections to improve G(r) at low r')
     4818        print('start: Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} (RMS:{:.2f})'.format(
     4819                data['Flat Bkg'],data['BackRatio'],data['Ruland'],rms))
     4820       
     4821        res = opt.minimize(Min,xstart,bounds=([0,None],[0,1],[0.01,1]),
     4822                           method='L-BFGS-B',options={'maxiter':5})
     4823        Done(res['x'])     
     4824        print('end:   Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} (RMS:{:.2f})\n'.format(
     4825                data['Flat Bkg'],data['BackRatio'],data['Ruland'],res['fun']))
     4826        wx.EndBusyCursor()
     4827        wx.CallAfter(UpdatePDFGrid,G2frame,data)
     4828        OnComputePDF(event)
     4829       
     4830    def SetupPDFEval():
     4831        '''Create functions needed to optimize the PDF at low r
     4832        '''
     4833        if not data['ElList']:
     4834            return None
     4835        Data = copy.deepcopy(data)
     4836        xydata = {}
     4837        for key in ['Sample','Sample Bkg.','Container','Container Bkg.']:
     4838            name = Data[key]['Name']
     4839            if name:
     4840                xydata[key] = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,G2frame.root,name))
     4841        powName = Data['Sample']['Name']
     4842        powId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,powName)
     4843        limits = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,powId,'Limits'))[1]
     4844        inst = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,powId,'Instrument Parameters'))[0]
     4845        numbDen = G2pwd.GetNumDensity(Data['ElList'],Data['Form Vol'])
     4846        BkgMax = 1.
     4847        def EvalLowPDF(arg):
     4848            '''Objective routine -- evaluates the RMS deviations in G(r)
     4849            from -4(pi)r for for r<Rmin
     4850            arguments are ['Flat Bkg','BackRatio','Ruland'] scaled so that
     4851            the min & max values are between 0 and 1.
     4852            '''
     4853            F,B,R = arg
     4854            Data['Flat Bkg'] = F*BkgMax
     4855            Data['BackRatio'] = B
     4856            Data['Ruland'] = R/10.
     4857            #for key in 'Flat Bkg','BackRatio','Ruland':
     4858            #    print key,Data[key],'; ',
     4859            G2pwd.CalcPDF(Data,inst,limits,xydata)
     4860            # test low r computation
     4861            g = xydata['GofR'][1][1]
     4862            r = xydata['GofR'][1][0]
     4863            g0 = g[r < Data['Rmin']] + 4*np.pi*r[r < Data['Rmin']]*numbDen
     4864            #print ' Res=',sum(g0**2)
     4865            return sum(g0**2)
     4866        def GetCurrentVals():
     4867            '''Get the current ['Flat Bkg','BackRatio','Ruland'] with scaling
     4868            '''
     4869            try:
     4870                F = data['Flat Bkg']/BkgMax
     4871            except:
     4872                F = 0
     4873            return [F,data['BackRatio'],max(10*data['Ruland'],.05)]
     4874        def SetFinalVals(arg):
     4875            '''Set the 'Flat Bkg', 'BackRatio' & 'Ruland' values from the
     4876            scaled, refined values and plot corrected region of G(r)
     4877            '''
     4878            F,B,R = arg
     4879            data['Flat Bkg'] = F*BkgMax
     4880            data['BackRatio'] = B
     4881            data['Ruland'] = R/10.
     4882            G2pwd.CalcPDF(data,inst,limits,xydata)
     4883            g = xydata['GofR'][1][1]
     4884            r = xydata['GofR'][1][0]
     4885            g0 = g[r < Data['Rmin']] + 4*np.pi*r[r < Data['Rmin']]*numbDen
     4886            G2plt.PlotXY(G2frame,[[r[r < Data['Rmin']], g0]],Title='G(r)+4pi*r',
     4887                         labelX=r'r, $\AA$',labelY=r'G(r)$+4\pi r$')           
     4888        EvalLowPDF(GetCurrentVals())
     4889        BkgMax = max(xydata['IofQ'][1][1])/50.
     4890        return EvalLowPDF,GetCurrentVals,SetFinalVals
    48074891               
    48084892    def ComputePDF(Data):
     4893        '''Calls :func:`GSASIIpwd.CalcPDF` to compute and the PDF. Called from OnComputePDF and OnComputeAllPDF
     4894        '''
     4895
    48094896        xydata = {}
    48104897        for key in ['Sample','Sample Bkg.','Container','Container Bkg.']:
     
    48254912       
    48264913    def OnComputePDF(event):
     4914        '''Compute and plot PDF, in response to a menu command or a change to a
     4915        computation parameter.
     4916        '''
    48274917        if not data['ElList']:
    48284918            G2frame.ErrorDialog('PDF error','Chemical formula not defined')
     
    48394929            G2plt.PlotISFG(G2frame,newPlot=True,plotType='S(Q)')
    48404930            G2plt.PlotISFG(G2frame,newPlot=True,plotType='F(Q)')
    4841             G2plt.PlotISFG(G2frame,newPlot=True,plotType='G(R)')
     4931            #G2plt.PlotISFG(G2frame,newPlot=True,plotType='G(R)')
     4932            G2plt.PlotISFG(G2frame,newPlot=False,plotType='G(R)')
    48424933        else:
    48434934            G2plt.PlotISFG(G2frame,newPlot=False)
     
    48994990    if 'IofQmin' not in data:
    49004991        data['IofQmin'] = 1.0       
     4992    if 'Rmin' not in data:
     4993        data['Rmin'] = 1.5
    49014994    if G2frame.dataDisplay:
    49024995        G2frame.dataFrame.Clear()
     
    49125005    G2frame.dataFrame.Bind(wx.EVT_MENU, OnComputePDF, id=G2gd.wxID_PDFCOMPUTE)
    49135006    G2frame.dataFrame.Bind(wx.EVT_MENU, OnComputeAllPDF, id=G2gd.wxID_PDFCOMPUTEALL)
     5007    G2frame.dataFrame.Bind(wx.EVT_MENU, OnOptimizePDF, id=G2gd.wxID_PDFOPT)
    49145008    mainSizer = wx.BoxSizer(wx.VERTICAL)
    49155009
     
    49735067       
    49745068    G2G.HorizontalLine(mainSizer,G2frame.dataDisplay)
    4975     mainSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' S(Q)->F(Q)->G(r) controls: '),0,WACV)
     5069    sqBox = wx.BoxSizer(wx.HORIZONTAL)
     5070    sqBox.Add(wx.StaticText(G2frame.dataDisplay,label=' S(Q)->F(Q)->G(r) controls: '),0,WACV)
     5071    sqBox.Add((1,1),1,wx.EXPAND,1)
     5072    optB = wx.Button(G2frame.dataDisplay,label='Optimize PDF',style=wx.BU_EXACTFIT)
     5073    optB.Bind(wx.EVT_BUTTON, OnOptimizePDF)
     5074    sqBox.Add(optB,0,WACV|wx.ALIGN_RIGHT)
     5075    mainSizer.Add(sqBox,0,wx.EXPAND)
     5076   
    49765077    mainSizer.Add((5,5),0)
    49775078    sqBox = wx.BoxSizer(wx.HORIZONTAL)
     
    49965097    flatSpin.Bind(wx.EVT_SPIN, OnFlatSpin)
    49975098    sqBox.Add(flatSpin,0,WACV)
    4998     mainSizer.Add(sqBox,0)
     5099    sqBox.Add((1,1),1,wx.EXPAND,1)
     5100    sqBox.Add(wx.StaticText(G2frame.dataDisplay,label='Rmin: '),0,WACV|wx.ALIGN_RIGHT)
     5101    rmin = G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data,'Rmin',nDig=(5,1),
     5102            typeHint=float,size=wx.Size(50,20))
     5103    sqBox.Add(rmin,0,WACV|wx.ALIGN_RIGHT)
     5104    mainSizer.Add(sqBox,0,wx.EXPAND)
    49995105       
    50005106    bkBox = wx.BoxSizer(wx.HORIZONTAL)
     
    50315137                                 typeHint=float,OnLeave=NewQmax)
    50325138    sqBox.Add(SQmax,0,WACV)
    5033     resetQ = wx.CheckBox(parent=G2frame.dataDisplay,label='Reset?')
     5139    resetQ = wx.Button(G2frame.dataDisplay,label='Reset?',style=wx.BU_EXACTFIT)
    50345140    sqBox.Add(resetQ,0,WACV)
    5035     resetQ.Bind(wx.EVT_CHECKBOX, OnResetQ)
     5141    resetQ.Bind(wx.EVT_BUTTON, OnResetQ)
    50365142    sqBox.Add(wx.StaticText(G2frame.dataDisplay,label=' Rmax: '),0,WACV)
    50375143    rmax = G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data,'Rmax',nDig=(10,1),min=10.,max=200.,
Note: See TracChangeset for help on using the changeset viewer.