Changeset 3913


Ignore:
Timestamp:
Apr 21, 2019 10:05:23 PM (4 years ago)
Author:
toby
Message:

Update Rietveld cycle-by-cycle

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIdataGUI.py

    r3909 r3913  
    38803880            self.GPXtree.Expand(Id)
    38813881            SelectDataTreeItem(self,Id)
     3882            self.GPXtree.SelectItem(Id)  # needed on OSX or item is not selected in tree; perhaps not needed elsewhere
    38823883        elif phaseId:
    38833884            SelectDataTreeItem(self,phaseId)
     3885            self.GPXtree.SelectItem(phaseId) # as before for OSX
    38843886        self.CheckNotebook()
    38853887        if self.dirname: os.chdir(self.dirname)           # to get Mac/Linux to change directory!
    … …  
    45554557            return
    45564558        if warnmsg:
    4557             print('Conflict between refinment flag settings and constraints:\n'+
     4559            print('Conflict between refinement flag settings and constraints:\n'+
    45584560                warnmsg+'\nRefinement not possible')
    45594561            self.ErrorDialog('Refinement Flag Error',
    … …  
    45714573        self.SaveTreeSetting() # save the current tree selection
    45724574        self.GPXtree.SaveExposedItems()             # save the exposed/hidden tree items
     4575        refPlotUpdate = G2plt.PlotPatterns(self,refineMode=True) # prepare for plot updating
    45734576        try:
    4574             OK,Msg = G2stMn.Refine(self.GSASprojectfile,dlg)    #Msg is Rvals dict if Ok=True
     4577            OK,Msg = G2stMn.Refine(self.GSASprojectfile,dlg,refPlotUpdate=refPlotUpdate)    #Msg is Rvals dict if Ok=True
    45754578        finally:
    45764579            dlg.Update(101.) # forces the Auto_Hide; needed after move w/Win & wx3.0
    … …  
    45894592            try:
    45904593                if dlg2.ShowModal() == wx.ID_OK:
     4594                    if refPlotUpdate: refPlotUpdate({},restore=True)
    45914595                    wx.CallAfter(self.reloadFromGPX)
     4596                else:
     4597                    if refPlotUpdate: refPlotUpdate({},restore=True)
    45924598            finally:
    45934599                dlg2.Destroy()
    45944600        else:
    45954601            self.ErrorDialog('Refinement error',Msg)
    4596 
     4602           
    45974603    def reloadFromGPX(self):
    45984604        '''Deletes current data tree & reloads it from GPX file (after a
    … …  
    46444650                try:
    46454651                    win.replotFunction(*win.replotArgs,**win.replotKWargs)
    4646                 except:
     4652                except Exception as msg:
    46474653                    if GSASIIpath.GetConfigValue('debug'):
    4648                         print('Error with args',win.replotArgs,win.replotKWargs)
     4654                        print('Error calling',win.replotFunction,'with args',
     4655                                  win.replotArgs,win.replotKWargs)
     4656                    if GSASIIpath.GetConfigValue('debug'):
     4657                        print(msg)
     4658                        GSASIIpath.IPyBreak()
    46494659        # delete any remaining plots that have not been updated and need a refresh (win.plotRequiresRedraw)
    46504660        for lbl,win in zip(self.G2plotNB.plotList,self.G2plotNB.panelList):
    … …  
    46654675            self.OnRefine(event)
    46664676            return
     4677        plotHist = self.GPXtree.GetItemText(self.PatternId)
    46674678        Id = GetGPXtreeItemId(self,self.root,'Sequential results')
    46684679        if not Id:
    … …  
    46944705            dlg.SetSize((int(Size[0]*1.2),Size[1])) # increase size a bit along x
    46954706        dlg.CenterOnParent()
     4707        # find 1st histogram to be refined
     4708        if 'Seq Data' in Controls:
     4709            histNames = Controls['Seq Data']
     4710        else: # patch from before Controls['Seq Data'] was implemented
     4711            histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
     4712        if Controls.get('Reverse Seq'):
     4713            histNames.reverse()
     4714        # select it
     4715        self.PatternId = GetGPXtreeItemId(self,self.root,histNames[0])       
     4716        refPlotUpdate = G2plt.PlotPatterns(self,refineMode=True) # prepare for plot updating
    46964717        try:
    4697             OK,Msg = G2stMn.SeqRefine(self.GSASprojectfile,dlg,G2plt.SequentialPlotPattern,self)    #Msg is Rvals dict if Ok=True
     4718            OK,Msg = G2stMn.SeqRefine(self.GSASprojectfile,dlg,refPlotUpdate) #Msg is Rvals dict if Ok=True
    46984719        finally:
    46994720            dlg.Update(101.) # forces the Auto_Hide; needed after move w/Win & wx3.0
    … …  
    47044725            try:
    47054726                if dlg.ShowModal() == wx.ID_OK:
     4727                    if refPlotUpdate: refPlotUpdate({},restore=True)
    47064728                    self.PickIdText = None  #force reload of PickId contents
    47074729                    self.GPXtree.DeleteChildren(self.root)
    … …  
    47124734                    self.GPXtree.RestoreExposedItems()
    47134735                    self.ResetPlots()
     4736                    self.PatternId = GetGPXtreeItemId(self,self.root,plotHist)
     4737                    SelectDataTreeItem(self,self.PatternId)
    47144738                    sId = GetGPXtreeItemId(self,self.root,'Sequential results')
    47154739                    SelectDataTreeItem(self,sId)
     4740                    self.GPXtree.SelectItem(sId)
     4741                else:
     4742                    if refPlotUpdate: refPlotUpdate({},restore=True)
    47164743            finally:
    47174744                dlg.Destroy()

  • trunk/GSASIImath.py

    r3900 r3913  
    9999    return res,nzero
    100100
    101 def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False):
    102    
     101def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False,refPlotUpdate=None):
    103102    """
    104103    Minimize the sum of squares of a function (:math:`f`) evaluated on a series of
    … …  
    209208            print (' Cycle: %d, Time: %.2fs, Chi**2: %.5g for %d obs., Lambda: %.3g,  Delta: %.3g'%(
    210209                icycle,time.time()-time0,chisq1,Nobs,lamMax,deltaChi2))
     210        Histograms = args[0][0]
     211        if refPlotUpdate is not None: refPlotUpdate(Histograms,icycle)   # update plot
    211212        if deltaChi2 < ftol:
    212213            ifConverged = True
    … …  
    233234        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1}]         
    234235           
    235 def HessianSVD(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False):
     236def HessianSVD(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False,refPlotUpdate=None):
    236237   
    237238    """
    … …  
    319320            print (' Cycle: %d, Time: %.2fs, Chi**2: %.5g, Delta: %.3g'%(
    320321                icycle,time.time()-time0,chisq1,deltaChi2))
     322        Histograms = args[0][0]
     323        if refPlotUpdate is not None: refPlotUpdate(Histograms,icycle)   # update plot
    321324        if deltaChi2 < ftol:
    322325            ifConverged = True

  • trunk/GSASIIplot.py

    r3909 r3913  
    17441744##### PlotPatterns
    17451745################################################################################
    1746 def SequentialPlotPattern(G2frame,refdata,histogram):
    1747     '''This is passed into :func:`GSASIIstrMain.SeqRefine` where it is used to
    1748     provide a plot of the current powder histogram just after a refinement. It
    1749     takes the old refinement information (Rfactors, curve locations, etc.) and
    1750     combines it with the refinement results in refdata and passes that to
    1751     :func:`PlotPatterns`
    1752     '''
    1753     if not histogram.startswith('PWDR'): return
    1754     pickId = G2frame.PickId
    1755     G2frame.PickId = G2frame.PatternId = G2gd.GetGPXtreeItemId(G2frame, G2frame.root, histogram)
    1756     treedata = G2frame.GPXtree.GetItemPyData(G2frame.PatternId)
    1757     PlotPatterns(G2frame,newPlot=True,plotType='PWDR',data=[treedata[0],refdata])
    1758     wx.Yield() # force a plot update (needed on Windows?)
    1759     G2frame.PickId = pickId
    1760    
    17611746def ReplotPattern(G2frame,newPlot,plotType,PatternName=None,PickName=None):
    17621747    '''This does the same as PlotPatterns except that it expects the information
    … …  
    17771762
    17781763def PlotPatterns(G2frame,newPlot=False,plotType='PWDR',data=None,
    1779                      extraKeys=[]):
     1764                     extraKeys=[],refineMode=False):
    17801765    '''Powder pattern plotting package - displays single or multiple powder patterns as intensity vs
    17811766    2-theta, q or TOF. Can display multiple patterns as "waterfall plots" or contour plots. Log I
    … …  
    25942579        else:
    25952580            Plot.set_ylim(CurLims['ylims'])
     2581        Page.toolbar.push_current()
    25962582        Plot.figure.canvas.draw()
    25972583        #GSASIIpath.IPyBreak()
    25982584       
    25992585    def onPlotFormat(event):
     2586        '''Change the appearance of the current plot'''
    26002587        changePlotSettings(G2frame,Plot)
     2588       
     2589    def refPlotUpdate(Histograms,cycle=None,restore=False):
     2590        '''called to update an existing plot during a Rietveld fit
     2591        '''
     2592        if restore:
     2593            (G2frame.SinglePlot,G2frame.Contour,G2frame.Weight,
     2594                G2frame.plusPlot,G2frame.SubBack,Page.plotStyle['logPlot']) = savedSettings
     2595            return
     2596
     2597        if plottingItem not in Histograms:
     2598            histoList = [i for i in Histograms.keys() if i.startswith('PWDR ')]
     2599            if len(histoList) == 0:
     2600                print('Skipping plot, no PWDR item found!')
     2601                return
     2602            plotItem = histoList[0]
     2603        else:
     2604            plotItem = plottingItem
     2605        xye = np.array(ma.getdata(Histograms[plotItem]['Data'])) # strips mask
     2606        xye0 = Histograms[plotItem]['Data'][0]
     2607        if Page.plotStyle['qPlot']:
     2608            X = 2.*np.pi/G2lat.Pos2dsp(Parms,xye0)
     2609            Ibeg = np.searchsorted(X,2.*np.pi/G2lat.Pos2dsp(Parms,limits[1][0]))
     2610            Ifin = np.searchsorted(X,2.*np.pi/G2lat.Pos2dsp(Parms,limits[1][1]))
     2611        elif Page.plotStyle['dPlot']:
     2612            X = G2lat.Pos2dsp(Parms,xye0)
     2613            Ibeg = np.searchsorted(X,G2lat.Pos2dsp(Parms,limits[1][1]))
     2614            Ifin = np.searchsorted(X,G2lat.Pos2dsp(Parms,limits[1][0]))
     2615        else:
     2616            X = copy.deepcopy(xye0)
     2617            Ibeg = np.searchsorted(X,limits[1][0])
     2618            Ifin = np.searchsorted(X,limits[1][1])
     2619        if Page.plotStyle['sqrtPlot']:
     2620            olderr = np.seterr(invalid='ignore') #get around sqrt(-ve) error
     2621            Y = np.where(xye[1]>=0.,np.sqrt(xye[1]),-np.sqrt(-xye[1]))
     2622            Z = np.where(xye[3]>=0.,np.sqrt(xye[3]),-np.sqrt(-xye[3]))
     2623            W = np.where(xye[4]>=0.,np.sqrt(xye[4]),-np.sqrt(-xye[4]))
     2624            #D = np.where(xye[5],(Y-Z),0.)-Page.plotStyle['delOffset']
     2625            np.seterr(invalid=olderr['invalid'])
     2626        else:
     2627            Y = copy.copy(xye[1])
     2628            Z = copy.copy(xye[3])
     2629            W = copy.copy(xye[4])
     2630            #D = xye[5]-Page.plotStyle['delOffset']  #powder background
     2631        DZ = (xye[1]-xye[3])*np.sqrt(xye[2])
     2632        DifLine[0].set_xdata(X[Ibeg:Ifin])
     2633        DifLine[0].set_ydata(DZ[Ibeg:Ifin])
     2634        Plot1.set_ylim((min(DZ[Ibeg:Ifin]),max(DZ[Ibeg:Ifin])))
     2635        CalcLine[0].set_xdata(X)
     2636        ObsLine[0].set_xdata(X)
     2637        BackLine[0].set_xdata(X)
     2638        CalcLine[0].set_ydata(Z)
     2639        ObsLine[0].set_ydata(Y)
     2640        BackLine[0].set_ydata(W)
     2641        if cycle:
     2642            Title = '{} cycle #{}'.format(plotItem,cycle)
     2643        else:
     2644            Title = plotItem
     2645        if Page.plotStyle['sqrtPlot']:
     2646            Plot.set_title(r'$\sqrt{I}$ for '+Title)
     2647        else:
     2648            Plot.set_title(Title)
     2649        Page.canvas.draw()
     2650
    26012651    #=====================================================================================
    26022652    # beginning PlotPatterns execution
    … …  
    26352685        G2frame.UseLimits = {'xlims':[False,False],'ylims':[False,False],
    26362686                                       'dylims':[False,False]}
     2687    #=====================================================================================
     2688    # code to setup for plotting Rietveld results. Turns off multiplot,
     2689    # sqrtplot, turn on + and weight plot, but sqrtPlot qPlot and dPlot are not changed.
     2690    # Magnification regions are ignored.
     2691    # the last-plotted histogram (from G2frame.PatternId) is used for this plotting
     2692    #    (except in seq. fitting)
     2693    # Returns a pointer to refPlotUpdate, which is used to update the plot when this
     2694    # returns
     2695    if refineMode:
     2696        plottingItem = G2frame.GPXtree.GetItemText(G2frame.PatternId)
     2697        # save settings to be restored after refinement with repPlotUpdate({},restore=True)
     2698        savedSettings = (G2frame.SinglePlot,G2frame.Contour,G2frame.Weight,
     2699                            G2frame.plusPlot,G2frame.SubBack,Page.plotStyle['logPlot'])
     2700        G2frame.SinglePlot = True
     2701        G2frame.Contour = False
     2702        G2frame.Weight = True
     2703        G2frame.plusPlot = True
     2704        G2frame.SubBack = False
     2705        Page.plotStyle['logPlot'] = False
     2706    #=====================================================================================
    26372707    if not new:
    26382708        G2frame.xylim = limits
    … …  
    29072977        if 'PWDR' in plottype and G2frame.SinglePlot and not (
    29082978                Page.plotStyle['logPlot'] or Page.plotStyle['sqrtPlot'] or G2frame.Contour):
    2909             magLineList = data[0].get('Magnification',[])
     2979            if not refineMode:
     2980                magLineList = data[0].get('Magnification',[])
    29102981            if ('C' in ParmList[0]['Type'][0] and Page.plotStyle['dPlot']) or \
    29112982                ('T' in ParmList[0]['Type'][0] and Page.plotStyle['qPlot']): # reversed regions relative to data order
    … …  
    30543125                    else:
    30553126                        Plot.set_ylim(bottom=np.min(np.trim_zeros(YB))/2.,top=np.max(Y)*2.)
     3127                # Matplotlib artist lists used for refPlotUpdate
     3128                ObsLine = None
     3129                CalcLine = None
     3130                BackLine = None
     3131                DifLine = None
    30563132                if G2frame.Weight:
    30573133                    Plot1.set_yscale("linear")                                                 
    … …  
    30953171                    if G2frame.SubBack:
    30963172                        if 'PWDR' in plottype:
    3097                             Plot.plot(Xum,Y-W,colors[0]+pP,picker=False,clip_on=Clip_on,label='_obs')  #Io-Ib
    3098                             Plot.plot(X,Z-W,colors[1],picker=False,label='_calc')               #Ic-Ib
     3173                            ObsLine = Plot.plot(Xum,Y-W,colors[0]+pP,picker=False,clip_on=Clip_on,label='_obs')  #Io-Ib
     3174                            CalcLine = Plot.plot(X,Z-W,colors[1],picker=False,label='_calc')               #Ic-Ib
    30993175                        else:
    31003176                            Plot.plot(X,YB,colors[0]+pP,picker=3.,clip_on=Clip_on,label='_obs')
    … …  
    31033179                        if 'PWDR' in plottype:
    31043180                            ObsLine = Plot.plot(Xum,Y,colors[0]+pP,picker=3.,clip_on=Clip_on,label='_obs')    #Io
    3105                             Plot.plot(X,Z,colors[1],picker=False,label='_calc')                 #Ic
     3181                            CalcLine = Plot.plot(X,Z,colors[1],picker=False,label='_calc')                 #Ic
    31063182                        else:
    31073183                            Plot.plot(X,YB,colors[0]+pP,picker=3.,clip_on=Clip_on,label='_obs')
    31083184                            Plot.plot(X,ZB,colors[1],picker=False,label='_calc')
    31093185                    if 'PWDR' in plottype and (G2frame.SinglePlot or G2frame.plusPlot):
    3110                         Plot.plot(X,W,colors[2],picker=False,label='_bkg')                 #Ib
     3186                        BackLine = Plot.plot(X,W,colors[2],picker=False,label='_bkg')                 #Ib
    31113187                        if not G2frame.Weight: DifLine = Plot.plot(X,D,colors[3],picker=1.,label='_diff')                 #Io-Ic
    31123188                    Plot.axhline(0.,color='k',label='_zero')
    … …  
    31873263                    Plot.axvline(hkl[-2],color=clr,dashes=(5,5))
    31883264        elif G2frame.GPXtree.GetItemText(PickId) in ['Reflection Lists'] or \
    3189             'PWDR' in G2frame.GPXtree.GetItemText(PickId):
     3265            'PWDR' in G2frame.GPXtree.GetItemText(PickId) or refineMode:
    31903266            Phases = G2frame.GPXtree.GetItemPyData(G2gd.GetGPXtreeItemId(G2frame,PatternId,'Reflection Lists'))
    31913267            l = GSASIIpath.GetConfigValue('Tick_length',8.0)
    … …  
    32813357        if DifLine[0]:
    32823358            G2frame.dataWindow.moveDiffCurve.Enable(True)
     3359    if refineMode: return refPlotUpdate
    32833360           
    32843361def PublishRietveldPlot(G2frame,Pattern,Plot,Page):

  • trunk/GSASIIstrMain.py

    r3855 r3913  
    4848
    4949def RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
    50     calcControls,pawleyLookup,ifSeq,printFile,dlg):
     50    calcControls,pawleyLookup,ifSeq,printFile,dlg,refPlotUpdate=None):
    5151    '''Core optimization routines, shared between SeqRefine and Refine
    5252
    … …  
    8080                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
    8181            ncyc = int(result[2]['nfev']/2)
     82            if refPlotUpdate is not None: refPlotUpdate(Histograms)   # update plot after completion
    8283        elif 'analytic Hessian' in Controls['deriv type']:
    8384            Lamda = Controls.get('Marquardt',-3)
    8485            maxCyc = Controls['max cyc']
    8586            result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,lamda=Lamda,
    86                 args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
     87                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
     88                refPlotUpdate=refPlotUpdate)
    8789            ncyc = result[2]['num cyc']+1
    8890            Rvals['lamMax'] = result[2]['lamMax']
    … …  
    9193            maxCyc = Controls['max cyc']
    9294            result = G2mth.HessianSVD(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,
    93                 args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
     95                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
     96                refPlotUpdate=refPlotUpdate)
    9497            if result[1] is None:
    9598                IfOK = False
    … …  
    104107            if len(varyList):
    105108                ncyc = int(result[2]['nfev']/len(varyList))
     109            if refPlotUpdate is not None: refPlotUpdate(Histograms)   # update plot
    106110#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
    107111#        for item in table: print item,table[item]               #useful debug - are things shifting?
    … …  
    173177    return IfOK,Rvals,result,covMatrix,sig
    174178
    175 def Refine(GPXfile,dlg=None,makeBack=True):
     179def Refine(GPXfile,dlg=None,makeBack=True,refPlotUpdate=None):
    176180    'Global refinement -- refines to minimize against all histograms'
    177181    import GSASIImpsubs as G2mp
    … …  
    249253        covData = {}
    250254        IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histograms,Phases,restraintDict,
    251             rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg)
     255            rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
     256            refPlotUpdate=refPlotUpdate)
    252257        if IfOK:
    253258            sigDict = dict(zip(varyList,sig))
    … …  
    315320    return NewVary
    316321
    317 def SeqRefine(GPXfile,dlg,PlotFunction=None,G2frame=None):
     322def SeqRefine(GPXfile,dlg,refPlotUpdate=None):
    318323    '''Perform a sequential refinement -- cycles through all selected histgrams,
    319324    one at a time
    … …  
    480485#        try:
    481486            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
    482                 rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg)
    483             if PlotFunction:
    484                 PlotFunction(G2frame,Histo[histogram]['Data'],histogram)
     487                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
     488                refPlotUpdate=refPlotUpdate)
    485489            print ('  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f'%(
    486490                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2']))

  • trunk/makeBat.py

    r3680 r3913  
    1313
    1414'''
     15version = "$Id$"
    1516# creates Windows files to aid in running GSAS-II
    1617#   creates RunGSASII.bat and a desktop shortcut to that file
Note: See TracChangeset for help on using the changeset viewer.