Changeset 4840


Ignore:
Timestamp:
Mar 5, 2021 1:54:34 PM (3 years ago)
Author:
toby
Message:

LeBail? issues: bug on zero cycles fixed; new LeBail? fit menu item; HessianLSQ: remove extra func eval on 0 cycles; fix noprint options; misc cleanups

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIdataGUI.py

    r4837 r4840  
    753753       
    754754        item = parent.Append(wx.ID_ANY,'&Refine\tCTRL+R','Perform a refinement')
    755         if len(self.Refine): # extend state for new menus to match main (on mac)
     755        if len(self.Refine): # extend state for new menus to match main
    756756            state = self.Refine[0].IsEnabled()
    757757        else:
    … …  
    760760        self.Refine.append(item)
    761761        self.Bind(wx.EVT_MENU, self.OnRefine, id=item.GetId())
     762        item = parent.Append(wx.ID_ANY,'&LeBail fit\tCTRL+B','Fit LeBail intensities only')
     763        item.Enable(state)
     764        self.Refine.append(item)
     765        self.Bind(wx.EVT_MENU, self.OnLeBail, id=item.GetId())
    762766
    763767        item = parent.Append(wx.ID_ANY,'&Run Fprime','X-ray resonant scattering')
    … …  
    797801            seqSetting = None
    798802           
     803        for item in self.Refine:
     804            if 'LeBail' in item.GetItemLabel():
     805                if seqSetting: item.Enable(False)
     806            elif seqSetting:
     807                item.SetItemLabel('Se&quential refine\tCtrl+R')    #might fail on old wx
     808            else:
     809                item.SetItemLabel('&Refine\tCtrl+R')    #might fail on old wx
    799810        if seqSetting:
    800             for item in self.Refine:
    801                 item.SetItemLabel('Se&quential refine\tCtrl+R')    #might fail on old wx
    802811            seqMode = True
    803812        else:
    804             for item in self.Refine:
    805                 item.SetItemLabel('&Refine\tCtrl+R')    #might fail on old wx
    806813            seqMode = False
    807814        for menu,Id in self.ExportSeq:
    … …  
    810817            menu.Enable(Id,not seqMode)
    811818        return seqSetting
    812 
    813819       
    814820    def PreviewFile(self,filename):
    … …  
    51495155        dlg.ShowModal()
    51505156        dlg.Destroy()
    5151 
     5157       
    51525158    def OnRefine(self,event):
    51535159        '''Perform a refinement or a sequential refinement (depending on controls setting)
    … …  
    52125218            dlg2 = wx.MessageDialog(self,text,'Refinement results, Rw =%.3f'%(Rw),wx.OK|wx.CANCEL)
    52135219            dlg2.CenterOnParent()
    5214             dlg2.Raise()
     5220            #dlg2.Raise()  # crashes sometimes on Mac
    52155221            try:
    52165222                if dlg2.ShowModal() == wx.ID_OK:
    … …  
    52385244        else:
    52395245            self.ErrorDialog('Refinement error',Rvals['msg'])
     5246           
     5247    def OnLeBail(self,event):
     5248        seqList = self.testSeqRefineMode()
     5249        if seqList:
     5250            self.ErrorDialog('Not for Sequential Fits',
     5251                'This command is not yet implemented for sequential fitting')
     5252            return
     5253        item = GetGPXtreeItemId(self,self.root,'Covariance')
     5254        covData = self.GPXtree.GetItemPyData(item)
     5255        GOF0 = covData['Rvals']['GOF']
     5256       
     5257        dlg = wx.ProgressDialog('Residual','All data Rw =',101.0,
     5258            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_CAN_ABORT|wx.STAY_ON_TOP,parent=self)
     5259        Size = dlg.GetSize()
     5260        if 50 < Size[0] < 500: # sanity check on size, since this fails w/Win & wx3.0
     5261            dlg.SetSize((int(Size[0]*1.2),Size[1])) # increase size a bit along x
     5262        dlg.CenterOnParent()
     5263        dlg.Raise()
     5264        self.SaveTreeSetting() # save the current tree selection
     5265        self.GPXtree.SaveExposedItems()             # save the exposed/hidden tree items
     5266        if self.PatternId and self.GPXtree.GetItemText(self.PatternId).startswith('PWDR '):
     5267            refPlotUpdate = G2plt.PlotPatterns(self,refineMode=True) # prepare for plot updating
     5268        else:
     5269            refPlotUpdate = None
     5270
     5271        try:
     5272            OK,Rvals = G2stMn.DoLeBail(self.GSASprojectfile,dlg,refPlotUpdate=refPlotUpdate)
     5273        finally:
     5274            dlg.Update(101.) # forces the Auto_Hide; needed after move w/Win & wx3.0
     5275            dlg.Destroy()
     5276        if OK:
     5277            text = ''
     5278            rtext = 'LeBail+only fit done. '
     5279            Rwp = Rvals.get('Rwp','?')
     5280            if 'GOF' in Rvals:
     5281                txt = 'Final Reduced Chi^2: {:.3f} (before ref: {:.3f})\n'.format(
     5282                    Rvals['GOF']**2,GOF0**2)
     5283                text += txt
     5284                rtext += txt
     5285            text += '\nLoad new result?'
     5286            dlg2 = wx.MessageDialog(self,text,'LeBail fit: Rwp={:.3f}'
     5287                                            .format(Rwp),wx.OK|wx.CANCEL)
     5288            dlg2.CenterOnParent()
     5289            try:
     5290                if dlg2.ShowModal() == wx.ID_OK:
     5291                    if refPlotUpdate: refPlotUpdate({},restore=True)
     5292                    wx.CallAfter(self.reloadFromGPX,rtext)
     5293                else:
     5294                    if refPlotUpdate: refPlotUpdate({},restore=True)
     5295            finally:
     5296                dlg2.Destroy()
     5297        else:
     5298            self.ErrorDialog('LeBail error',Rvals['msg'])
    52405299           
    52415300    def reloadFromGPX(self,rtext=None):

  • trunk/GSASIImath.py

    r4839 r4840  
    359359    nfev += 1
    360360    try:
    361         M = func(x0,*args)
     361        if icycle == 0: # no parameter changes, skip recalc
     362            M = M2
     363        else:
     364            M = func(x0,*args)
    362365    except Exception as Msg:
    363366        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)

  • trunk/GSASIIplot.py

    r4837 r4840  
    67576757    G2frame.G2plotNB.status.SetStatusText('',1)
    67586758    G2frame.G2plotNB.status.SetStatusWidths([G2frame.G2plotNB.status.firstLen,-1])
    6759     acolor = mpl.cm.get_cmap(G2frame.VcovColor)
    6760     Img = Plot.imshow(Page.covArray,aspect='equal',cmap=acolor,interpolation='nearest',origin='lower',
    6761         vmin=-1.,vmax=1.)
    6762     imgAx = Img.axes
    6763     ytics = imgAx.get_yticks()
    6764     ylabs = [Page.varyList[int(i)] for i in ytics[:-1]]
    6765     imgAx.set_yticklabels(ylabs)
    6766     Page.figure.colorbar(Img)
    6767     Plot.set_title('V-Cov matrix'+title)
    6768     Plot.set_xlabel('Variable number')
    6769     Plot.set_ylabel('Variable name')
     6759    if Page.varyList:
     6760        acolor = mpl.cm.get_cmap(G2frame.VcovColor)
     6761        Img = Plot.imshow(Page.covArray,aspect='equal',cmap=acolor,interpolation='nearest',origin='lower',
     6762                          vmin=-1.,vmax=1.)
     6763        imgAx = Img.axes
     6764        ytics = imgAx.get_yticks()
     6765        ylabs = [Page.varyList[int(i)] for i in ytics[:-1]]
     6766        imgAx.set_yticklabels(ylabs)
     6767        Page.figure.colorbar(Img)
     6768        Plot.set_title('V-Cov matrix'+title)
     6769        Plot.set_xlabel('Variable number')
     6770        Plot.set_ylabel('Variable name')
    67706771    Page.canvas.draw()
    67716772   

  • trunk/GSASIIstrIO.py

    r4839 r4840  
    10761076            (RBModel['rbRef'][0],RBModel['rbRef'][1],RBModel['rbRef'][0],RBModel['rbRef'][2]))
    10771077           
     1078    if Print and pFile is None: raise Exception("specify pFile or Print=False")
    10781079    rbVary = []
    10791080    rbDict = {}
    … …  
    14321433                phaseVary += [name,]
    14331434                   
     1435    if Print and pFile is None: raise Exception("specify pFile or Print=False")
    14341436    if Print:
    14351437        pFile.write('\n Phases:\n')
    … …  
    21602162    ##########################################################################
    21612163    # SetPhaseData starts here
    2162     pFile.write('\n Phases:\n')
     2164    if pFile: pFile.write('\n Phases:\n')
    21632165    for phase in Phases:
    2164         pFile.write(' Result for phase: %s\n'%phase)
    2165         pFile.write(135*'='+'\n')
     2166        if pFile: pFile.write(' Result for phase: %s\n'%phase)
     2167        if pFile: pFile.write(135*'='+'\n')
    21662168        Phase = Phases[phase]
    21672169        General = Phase['General']
    … …  
    21782180            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
    21792181            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
    2180             pFile.write(' Reciprocal metric tensor: \n')
     2182            if pFile: pFile.write(' Reciprocal metric tensor: \n')
    21812183            ptfmt = "%15.9f"
    21822184            names = ['A11','A22','A33','A12','A13','A23']
    … …  
    21912193                else:
    21922194                    sigstr += 15*' '
    2193             pFile.write(namstr+'\n')
    2194             pFile.write(ptstr+'\n')
    2195             pFile.write(sigstr+'\n')
     2195            if pFile: pFile.write(namstr+'\n')
     2196            if pFile: pFile.write(ptstr+'\n')
     2197            if pFile: pFile.write(sigstr+'\n')
    21962198            cell[1:7] = G2lat.A2cell(A)
    21972199            cell[7] = G2lat.calc_V(A)
    2198             pFile.write(' New unit cell:\n')
     2200            if pFile: pFile.write(' New unit cell:\n')
    21992201            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
    22002202            names = ['a','b','c','alpha','beta','gamma','Volume']
    … …  
    22092211                else:
    22102212                    sigstr += 12*' '
    2211             pFile.write(namstr+'\n')
    2212             pFile.write(ptstr+'\n')
    2213             pFile.write(sigstr+'\n')
     2213            if pFile: pFile.write(namstr+'\n')
     2214            if pFile: pFile.write(ptstr+'\n')
     2215            if pFile: pFile.write(sigstr+'\n')
    22142216        ik = 6  #for Pawley stuff below
    22152217        if General.get('Modulated',False):
    … …  
    22172219            Vec,vRef,maxH = General['SuperVec']
    22182220            if vRef:
    2219                 pFile.write(' New modulation vector:\n')
     2221                if pFile: pFile.write(' New modulation vector:\n')
    22202222                namstr = '  names :'
    22212223                ptstr =  '  values:'
    … …  
    22292231                    else:
    22302232                        sigstr += 12*' '
    2231                 pFile.write(namstr+'\n')
    2232                 pFile.write(ptstr+'\n')
    2233                 pFile.write(sigstr+'\n')
     2233                if pFile: pFile.write(namstr+'\n')
     2234                if pFile: pFile.write(ptstr+'\n')
     2235                if pFile: pFile.write(sigstr+'\n')
    22342236           
    22352237        General['Mass'] = 0.
    … …  
    22472249            RRBIds = RBIds['Residue']
    22482250            RBModels = Phase['RBModels']
    2249             for irb,RBObj in enumerate(RBModels.get('Vector',[])):
    2250                 jrb = VRBIds.index(RBObj['RBId'])
    2251                 rbsx = str(irb)+':'+str(jrb)
    2252                 pFile.write(' Vector rigid body parameters:\n')
    2253                 PrintRBObjPOAndSig('RBV',rbsx)
    2254                 PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
    2255             for irb,RBObj in enumerate(RBModels.get('Residue',[])):
    2256                 jrb = RRBIds.index(RBObj['RBId'])
    2257                 rbsx = str(irb)+':'+str(jrb)
    2258                 pFile.write(' Residue rigid body parameters:\n')
    2259                 PrintRBObjPOAndSig('RBR',rbsx)
    2260                 PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
    2261                 PrintRBObjTorAndSig(rbsx)
     2251            if pFile:
     2252                for irb,RBObj in enumerate(RBModels.get('Vector',[])):
     2253                    jrb = VRBIds.index(RBObj['RBId'])
     2254                    rbsx = str(irb)+':'+str(jrb)
     2255                    pFile.write(' Vector rigid body parameters:\n')
     2256                    PrintRBObjPOAndSig('RBV',rbsx)
     2257                    PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
     2258                for irb,RBObj in enumerate(RBModels.get('Residue',[])):
     2259                    jrb = RRBIds.index(RBObj['RBId'])
     2260                    rbsx = str(irb)+':'+str(jrb)
     2261                    pFile.write(' Residue rigid body parameters:\n')
     2262                    PrintRBObjPOAndSig('RBR',rbsx)
     2263                    PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
     2264                    PrintRBObjTorAndSig(rbsx)
    22622265            atomsSig = {}
    22632266            wavesSig = {}
    … …  
    23252328                                if pfx+name in sigDict:
    23262329                                    wavesSig[name] = sigDict[pfx+name]
    2327                    
    2328             PrintAtomsAndSig(General,Atoms,atomsSig)
    2329             if General['Type'] == 'magnetic':
     2330
     2331            if pFile: PrintAtomsAndSig(General,Atoms,atomsSig)
     2332            if pFile and General['Type'] == 'magnetic':
    23302333                PrintMomentsAndSig(General,Atoms,atomsSig)
    2331             if General.get('Modulated',False):
     2334            if pFile and General.get('Modulated',False):
    23322335                PrintWavesAndSig(General,Atoms,wavesSig)
    23332336               
    23342337            density = G2mth.getDensity(General)[0]
    2335             pFile.write('\n Density: {:.4f} g/cm**3\n'.format(density))
     2338            if pFile: pFile.write('\n Density: {:.4f} g/cm**3\n'.format(density))
    23362339           
    23372340       
    … …  
    36143617            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
    36153618
    3616             if not seq:
     3619            if Print and not seq:
    36173620                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
    36183621                pFile.write(135*'='+'\n')
    3619             pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
    3620             pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
     3622            if Print:
     3623                pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
     3624                pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
    36213625                (Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
    3622             pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
     3626                pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
    36233627                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
    3624             if Print:
    36253628                pFile.write(' Instrument type: %s\n'%Sample['Type'])
    36263629                if calcControls != None:    #skipped in seqRefine

  • trunk/GSASIIstrMain.py

    r4825 r4840  
    288288                ReportProblems(result,Rvals,varyList)
    289289                # process singular variables
    290                 psing = result[2].get('psing',[])
    291290                if dlg: break # refining interactively
    292                 # non-interactive refinement
    293                 #for val in sorted(psing,reverse=True):
    294                 #    G2fil.G2Print ('Removing parameter: '+varyList[val])
    295                 #    del(varyList[val])
    296                 #if not psing: break    # removed variable(s), try again
    297291            else:
    298292                G2fil.G2Print ('**** Refinement failed - singular matrix ****',mode='error')
    … …  
    304298                        break
    305299    if IfOK:
     300        if CheckLeBail(Phases):   # only needed for LeBail extraction
     301            G2stMth.errRefine(values,[Histograms,Phases,restraintDict,rigidbodyDict],
     302                                parmDict,varyList,calcControls,pawleyLookup,dlg)
     303
    306304        G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
    307305    if chisq0 is not None:
    … …  
    483481            Rvals['msg'] = Msg.msg
    484482        return False,Rvals
    485 
     483   
    486484#for testing purposes, create a file for testderiv
    487485    if GSASIIpath.GetConfigValue('debug'):   # and IfOK:
    … …  
    501499        G2fil.G2Print ('Reported from refinement:',mode='warn')
    502500        G2fil.G2Print (Rvals['msg'],mode='warn')
     501
     502def CheckLeBail(Phases):
     503    '''Check if there is a LeBail extraction in any histogram
     504
     505    :returns: True if there is at least one LeBail flag turned on, False otherwise
     506    '''
     507    for key in Phases:
     508        phase = Phases[key]
     509        for h in phase['Histograms']:
     510            #phase['Histograms'][h]
     511            if not phase['Histograms'][h]['Use']: continue
     512            if phase['Histograms'][h]['LeBail']:
     513                 return True
     514    return False
     515       
     516def DoLeBail(GPXfile,dlg=None,cycles=3,refPlotUpdate=None):
     517    '''Fit LeBail intensities without changes to any other refined parameters.
     518    This is a stripped-down version of :func:`Refine` that does not perform
     519    any refinement cycles
     520    '''
     521    import GSASIImpsubs as G2mp
     522    G2mp.InitMP()
     523    import pytexture as ptx
     524    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
     525
     526    #varyList = []
     527    parmDict = {}
     528    Controls = G2stIO.GetControls(GPXfile)
     529    calcControls = {}
     530    calcControls.update(Controls)
     531    constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
     532    restraintDict = {}
     533    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
     534    if not Phases:
     535        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
     536        return False,{'msg':'No phases'}
     537    if not Histograms:
     538        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
     539        return False,{'msg':'No data'}
     540    if not CheckLeBail(Phases):
     541        msg = 'Warning: There are no histograms with LeBail extraction enabled'
     542        G2fil.G2Print ('*** '+msg+' ***')
     543        return False,{'msg': msg}
     544    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
     545    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
     546    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,Print=False)
     547    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
     548         maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,Print=False)
     549    calcControls['atomIndx'] = atomIndx
     550    calcControls['Natoms'] = Natoms
     551    calcControls['FFtables'] = FFtables
     552    calcControls['BLtables'] = BLtables
     553    calcControls['MFtables'] = MFtables
     554    calcControls['maxSSwave'] = maxSSwave
     555    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,Print=False)
     556    calcControls.update(controlDict)
     557    histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,Print=False)
     558    calcControls.update(controlDict)
     559    parmDict.update(rbDict)
     560    parmDict.update(phaseDict)
     561    parmDict.update(hapDict)
     562    parmDict.update(histDict)
     563    G2stIO.GetFprime(calcControls,Histograms)
     564    parmFrozenList = Controls['parmFrozen']['FrozenList']
     565    try:
     566        for i in range(cycles):
     567            M = G2stMth.errRefine([],[Histograms,Phases,restraintDict,rigidbodyDict],parmDict,[],calcControls,pawleyLookup,dlg)
     568            G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
     569            if refPlotUpdate is not None: refPlotUpdate(Histograms,i)
     570        Rvals = {}
     571        Rvals['chisq'] = np.sum(M**2)
     572        Rvals['Nobs'] = Histograms['Nobs']
     573        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
     574        Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs'])) # no variables
     575
     576        covData = {'variables':0,'varyList':[],'sig':[],'Rvals':Rvals,
     577                       'varyListStart':[],
     578                       'covMatrix':np.zeros([0,0]),'title':GPXfile,
     579                       #'newAtomDict':newAtomDict,'newCellDict':newCellDict,
     580                       'freshCOV':True}
     581       
     582        G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,parmFrozenList,True)
     583        G2fil.G2Print (' ***** LeBail fit completed *****')
     584        return True,Rvals
     585    except Exception as Msg:
     586        G2fil.G2Print (' ***** LeBail fit error *****')
     587        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
     588        if GSASIIpath.GetConfigValue('debug'):
     589            import traceback
     590            print(traceback.format_exc())       
     591        return False,{'msg':Msg.msg}
    503592
    504593def phaseCheck(phaseVary,Phases,histogram):

  • trunk/GSASIIstrMath.py

    r4834 r4840  
    496496        if 'General' == pName.split(':')[1]:
    497497            # initialize for General restraint(s) here
    498             GeneralInit = True
    499498            parmDict0 = parmDict.copy()
    500499            # setup steps for each parameter
    … …  
    39443943        else:
    39453944            dMdV = np.concatenate((dMdV.T,np.sqrt(wtFactor)*dMdvh.T)).T
    3946            
     3945
    39473946    GetFobsSq(Histograms,Phases,parmDict,calcControls)
    39483947    pNames,pVals,pWt,pWsum,pWnum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
Note: See TracChangeset for help on using the changeset viewer.