Changeset 4076


Ignore:
Timestamp:
Aug 2, 2019 1:33:41 PM (4 years ago)
Author:
vondreele
Message:

fix axes amnesia for PDF plots
add plotting of P(R) pair distribution for small angle data
force element selection if none for 1st pdf calc
Size & pair dist plots now hollow bars
implement Moore method for making pair distribution P(R) from protein small angle data

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIplot.py

    r4074 r4076  
    588588                limits = [Plot.get_xlim(),Plot.get_ylim()] # save previous limits
    589589                if len(Axes)>1 and Plot.get_aspect() == 'auto':  # aspect will be equal for image plotting
    590                     limits[1] = Axes[1].get_ylim()
     590                    if Axes[1].get_aspect() != 'auto':      #not colorbars!
     591                        limits[1] = Axes[0].get_ylim()
     592                    else:
     593                        limits[1] = Axes[1].get_ylim()
    591594                if newImage:
    592595                    Page.figure.clf()
    … …  
    28032806            kwargs['PickName'] = G2frame.GPXtree.GetItemText(G2frame.PickId)
    28042807        G2frame.G2plotNB.RegisterRedrawRoutine(G2frame.G2plotNB.lastRaisedPlotTab,ReplotPattern,
    2805                                             (G2frame,newPlot,plotType),kwargs)
     2808            (G2frame,newPlot,plotType),kwargs)
    28062809    # now start plotting
    28072810    G2frame.G2plotNB.status.DestroyChildren() #get rid of special stuff on status bar
    … …  
    45614564    if not new:
    45624565        if not newPlot:
    4563             xylim = lim
     4566            xylim = copy.copy(lim)
    45644567    else:
    45654568        newPlot = True
    … …  
    52205223        elif 'Size' in PlotText:
    52215224            PlotSASDSizeDist(G2frame)
     5225        elif 'Pair' in PlotText:
     5226            PlotSASDPairDist(G2frame)
    52225227   
    52235228    def OnMotion(event):
    … …  
    52295234                G2frame.G2plotNB.status.SetStatusText('diameter =%9.3f f(D) =%9.3g'%(xpos,ypos),1)                   
    52305235            except TypeError:
    5231                 G2frame.G2plotNB.status.SetStatusText('Select Strain pattern first',1)
     5236                G2frame.G2plotNB.status.SetStatusText('Select Size pattern first',1)
    52325237
    52335238    new,plotNum,Page,Plot,lim = G2frame.G2plotNB.FindPlotTab('Size Distribution','mpl')
    … …  
    52455250        Plot.set_xscale("log",nonposx='mask')
    52465251        Plot.set_xlim([np.min(2.*Bins)/2.,np.max(2.*Bins)*2.])
    5247     Plot.bar(2.*Bins-Dbins,BinMag,2.*Dbins,facecolor='green')       #plot diameters
     5252    Plot.bar(2.*Bins-Dbins,BinMag,2.*Dbins,facecolor='white',edgecolor='green')       #plot diameters
    52485253    colors=['b','r','c','m','k']
    52495254    if 'Size Calc' in data:
    … …  
    52525257            if len(Rbins[i]):
    52535258                Plot.plot(2.*Rbins[i],Dist[i],color=colors[i%5])       #plot diameters
     5259    Page.canvas.draw()
     5260
     5261################################################################################
     5262##### PlotSASDPairDist
     5263################################################################################
     5264           
     5265def PlotSASDPairDist(G2frame):
     5266    'Needs a description'
     5267   
     5268    def OnPageChanged(event):
     5269        PlotText = G2frame.G2plotNB.nb.GetPageText(G2frame.G2plotNB.nb.GetSelection())
     5270        if 'Powder' in PlotText:
     5271            PlotPatterns(G2frame,plotType='SASD',newPlot=True)
     5272        elif 'Size' in PlotText:
     5273            PlotSASDSizeDist(G2frame)
     5274        elif 'Pair' in PlotText:
     5275            PlotSASDPairDist(G2frame)
     5276   
     5277    def OnMotion(event):
     5278        xpos = event.xdata
     5279        if xpos:                                        #avoid out of frame mouse position
     5280            ypos = event.ydata
     5281            SetCursor(Page)
     5282            try:
     5283                G2frame.G2plotNB.status.SetStatusText('pair dist =%9.3f f(D) =%9.3g'%(xpos,ypos),1)                   
     5284            except TypeError:
     5285                G2frame.G2plotNB.status.SetStatusText('Select Pair pattern first',1)
     5286
     5287    new,plotNum,Page,Plot,lim = G2frame.G2plotNB.FindPlotTab('Pair Distribution','mpl')
     5288    if new:
     5289        G2frame.G2plotNB.Bind(wx.aui.EVT_AUINOTEBOOK_PAGE_CHANGED,OnPageChanged)
     5290        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
     5291    Page.Choice = None
     5292    PatternId = G2frame.PatternId
     5293    data = G2frame.GPXtree.GetItemPyData(G2gd.GetGPXtreeItemId(G2frame,PatternId, 'Models'))
     5294    Bins,Dbins,BinMag = data['Pair']['Distribution']
     5295    Plot.set_title('Pair Distribution')
     5296    Plot.set_xlabel(r'$R, \AA$',fontsize=14)
     5297    Plot.set_ylabel(r'$Pair\ distribution,\ P(R)$',fontsize=14)
     5298    Plot.bar(Bins-Dbins,BinMag,Dbins,facecolor='white',edgecolor='green')       #plot diameters
     5299    colors=['b','r','c','m','k']
    52545300    Page.canvas.draw()
    52555301

  • trunk/GSASIIpwdGUI.py

    r4068 r4076  
    351351def SetDefaultSASDModel():
    352352    'Fills in default items for the SASD Models dictionary'   
    353     return {'Back':[0.0,False],'Size':{'MinDiam':50,'MaxDiam':10000,'Nbins':100,'logBins':True,'Method':'MaxEnt','Distribution':[],
    354         'Shape':['Spheroid',1.0],'MaxEnt':{'Niter':100,'Precision':0.01,'Sky':-3},
    355         'IPG':{'Niter':100,'Approach':0.8,'Power':-1},'Reg':{},},           
     353    return {'Back':[0.0,False],
     354        'Size':{'MinDiam':50,'MaxDiam':10000,'Nbins':100,'logBins':True,'Method':'MaxEnt',
     355                'Distribution':[],'Shape':['Spheroid',1.0],
     356                'MaxEnt':{'Niter':100,'Precision':0.01,'Sky':-3},
     357                'IPG':{'Niter':100,'Approach':0.8,'Power':-1},'Reg':{},},
     358        'Pair':{'Method':'Regularization','MaxRadius':100.,'NBins':100,'Errors':'User',
     359                'Percent error':2.5,'Background':[0,False],'Distribution':[],
     360                'Moore':20,'Dist G':100.,},           
    356361        'Particle':{'Matrix':{'Name':'vacuum','VolFrac':[0.0,False]},'Levels':[],},
    357362        'Current':'Size dist.','BackFile':'',
    … …  
    52165221    if 'BackFile' not in data:
    52175222        data['BackFile'] = ''
     5223    if 'Pair' not in data:
     5224        data['Pair'] = {'Method':'Regularization','MaxRadius':100.,'NBins':100,'Errors':'User',
     5225            'Percent error':2.5,'Background':[0,False],'Distribution':[],'Moore':[20,False],'Dist G':100.,}           
     5226
    52185227    #end patches
    52195228   
    … …  
    53985407            wx.CallAfter(UpdateModelsGrid,G2frame,data)
    53995408           
     5409        elif data['Current'] == 'Pair distance':
     5410            SaveState()
     5411            G2sasd.PairDistFxn(Profile,ProfDict,Limits,Sample,data)
     5412            G2plt.PlotSASDPairDist(G2frame)
     5413            RefreshPlots(True)
     5414            wx.CallAfter(UpdateModelsGrid,G2frame,data)
     5415           
    54005416    def OnUnDo(event):
    54015417        DoUnDo()
    … …  
    55425558        return sizeSizer
    55435559       
     5560    def PairSizer():
     5561#        'Pair':{'Method':'Regularization','MaxRadius':[100.,False],'NBins':101,'Errors':'User',
     5562#                'Percent error':2.5,'Background':[0,False],'Distribution':[],
     5563#                'Moore':20,},           
     5564               
     5565        def OnMethod(event):
     5566            data['Pair']['Method'] = method.GetValue()
     5567            wx.CallAfter(UpdateModelsGrid,G2frame,data)
     5568           
     5569        def OnError(event):
     5570            data['Pair']['Errors'] = error.GetValue()
     5571            wx.CallAfter(UpdateModelsGrid,G2frame,data)
     5572           
     5573        def OnMaxRadEst(event):
     5574            Results = G2sasd.RgFit(Profile,ProfDict,Limits,Sample,data)
     5575            if not Results[0]:
     5576                    G2frame.ErrorDialog('Failed refinement',
     5577                        ' Msg: '+Results[-1]+'\nYou need to rethink your selection of parameters\n'+    \
     5578                        ' Model restored to previous version')
     5579#            G2sasd.PairFxn(Profile,ProfDict,Limits,Sample,data)
     5580            RefreshPlots(True)
     5581            wx.CallAfter(UpdateModelsGrid,G2frame,data)                   
     5582                       
     5583        def OnMooreTerms(event):
     5584            data['Pair']['Moore'] = int(round(Limits[1][1]*data['Pair']['MaxRadius']/np.pi))-1
     5585            wx.CallAfter(UpdateModelsGrid,G2frame,data)                   
     5586           
     5587           
     5588        pairSizer = wx.BoxSizer(wx.VERTICAL)
     5589        pairSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Pair distribution parameters: '),0,WACV)
     5590        binSizer = wx.FlexGridSizer(0,6,5,5)
     5591        binSizer.Add(wx.StaticText(G2frame.dataWindow,label=' No. R bins: '),0,WACV)
     5592        bins = ['50','100','150','200']
     5593        nbins = wx.ComboBox(G2frame.dataWindow,value=str(data['Pair']['NBins']),choices=bins,
     5594            style=wx.CB_READONLY|wx.CB_DROPDOWN)
     5595        Indx[nbins.GetId()] = [data['Pair'],'NBins',0]
     5596        nbins.Bind(wx.EVT_COMBOBOX,OnIntVal)       
     5597        binSizer.Add(nbins,0,WACV)
     5598        binSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Max diam.: '),0,WACV)
     5599        maxdiam = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data['Pair'],'MaxRadius',min=10.,nDig=(10,1))
     5600        binSizer.Add(maxdiam,0,WACV)
     5601        maxest = wx.Button(G2frame.dataWindow,label='Make estimate')
     5602        maxest.Bind(wx.EVT_BUTTON,OnMaxRadEst)
     5603        binSizer.Add(maxest,0,WACV)
     5604        pairSizer.Add(binSizer,0)
     5605        pairSizer.Add((5,5),0)
     5606        fitSizer = wx.BoxSizer(wx.HORIZONTAL)
     5607        methods = ['Regularization','Moore']
     5608        fitSizer.Add(wx.StaticText(G2frame.dataWindow,label='Fitting method: '),0,WACV)
     5609        method = wx.ComboBox(G2frame.dataWindow,value=data['Pair']['Method'],choices=methods,
     5610            style=wx.CB_READONLY|wx.CB_DROPDOWN)
     5611        method.Bind(wx.EVT_COMBOBOX,OnMethod)
     5612        fitSizer.Add(method,0,WACV)
     5613        pairSizer.Add(fitSizer,0,WACV)
     5614        if 'Moore' in data['Pair']['Method']:
     5615            mooreSizer = wx.BoxSizer(wx.HORIZONTAL)
     5616            mooreSizer.Add(wx.StaticText(G2frame.dataWindow,label='Number of functions: '),0,WACV)
     5617            moore = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data['Pair'],'Moore',min=10)
     5618            mooreSizer.Add(moore,0,WACV)
     5619            mooreterms = wx.Button(G2frame.dataWindow,label = 'Auto determine?')
     5620            mooreterms.Bind(wx.EVT_BUTTON,OnMooreTerms)
     5621            mooreSizer.Add(mooreterms,0,WACV)
     5622            pairSizer.Add(mooreSizer,0,WACV)
     5623        errorSizer = wx.BoxSizer(wx.HORIZONTAL)
     5624        errorSizer.Add(wx.StaticText(G2frame.dataWindow,label='Error method: '),0,WACV)
     5625        errors = ['User','Sqrt','Percent']
     5626        error = wx.ComboBox(G2frame.dataWindow,value=data['Pair']['Errors'],choices=errors,
     5627            style=wx.CB_READONLY|wx.CB_DROPDOWN)
     5628        error.Bind(wx.EVT_COMBOBOX,OnError)
     5629        if 'Percent' in data['Pair']['Errors']:
     5630            percent = G2G.ValidatedTxtCtrl(G2frame.dataWindow,data['Pair'],'Percent error',min=0.5,nDig=(10,1))
     5631            errorSizer.Add(percent,0,WACV)
     5632        errorSizer.Add(error,0,WACV)
     5633        pairSizer.Add(errorSizer,0,WACV)
     5634        return pairSizer   
     5635   
     5636    def ShapeSizer():
     5637        shapeSizer = wx.BoxSizer(wx.VERTICAL)
     5638        shapeSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Shape parameters: TBD'),0,WACV)
     5639       
     5640       
     5641        return shapeSizer
     5642       
     5643       
     5644
    55445645    def PartSizer():
    55455646       
    … …  
    57505851                            parmSldr.Bind(wx.EVT_SLIDER,OnParmSlider)
    57515852                            parmSizer.Add(parmSldr,1,wx.EXPAND)
    5752             return parmSizer               
    5753            
     5853            return parmSizer
     5854
    57545855        Indx = {}
    57555856        partSizer = wx.BoxSizer(wx.VERTICAL)
    … …  
    58395940    mainSizer = G2frame.dataWindow.GetSizer()
    58405941    topSizer = wx.BoxSizer(wx.HORIZONTAL)
    5841     models = ['Size dist.','Particle fit']
     5942    models = ['Size dist.','Particle fit','Pair distance',]
     5943    if len(data['Pair']['Distribution']):
     5944        models += ['Shape',]
    58425945    topSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Modeling by: '),0,WACV)
    58435946    fitSel = wx.ComboBox(G2frame.dataWindow,value=data['Current'],choices=models,
    … …  
    58615964    elif 'Particle' in data['Current']:
    58625965        mainSizer.Add(PartSizer(),1,wx.ALIGN_LEFT|wx.EXPAND)
     5966    elif 'Pair' in data['Current']:
     5967        mainSizer.Add(PairSizer(),1,wx.ALIGN_LEFT|wx.EXPAND)
     5968    elif 'Shape' in data['Current']:
     5969        mainSizer.Add(ShapeSizer(),1,wx.ALIGN_LEFT|wx.EXPAND)
    58635970    G2G.HorizontalLine(mainSizer,G2frame.dataWindow)   
    58645971    backSizer = wx.BoxSizer(wx.HORIZONTAL)
    … …  
    70027109        if not data['ElList']:
    70037110            G2frame.ErrorDialog('PDF error','Chemical formula not defined')
    7004             return
     7111            OnAddElement(event)
    70057112        auxPlot = computePDF(G2frame,data)
    70067113        if auxPlot is None: return

  • trunk/GSASIIrestrGUI.py

    r4019 r4076  
    18821882    def UpdateGeneralRestr(generalRestData):
    18831883        '''Display any generalized restraint expressions'''
     1884       
    18841885        def OnEditGenRestraint(event):
    18851886            '''Edit a restraint expression'''
    … …  
    18871888            parmDict = SetupParmDict(G2frame)
    18881889            dlg = G2exG.ExpressionDialog(G2frame,parmDict,
    1889                             exprObj=generalRestData['General'][n][0],
    1890                             header="Edit a restraint expression",
    1891                             fit=False,wildCard=G2frame.testSeqRefineMode())
     1890                exprObj=generalRestData['General'][n][0],
     1891                header="Edit a restraint expression",
     1892                fit=False,wildCard=G2frame.testSeqRefineMode())
    18921893            restobj = dlg.Show(True)
    18931894            if restobj:
    18941895                generalRestData['General'][n][0] = restobj
    18951896                wx.CallAfter(UpdateGeneralRestr,restrData['General'])
    1896         def OnDelGenRestraint(event):
     1897               
     1898        def OnDelGenRestraint(event):              #does this work??
    18971899            '''Delete a restraint expression'''
    18981900            n = event.GetEventObject().index
    … …  
    19081910        hSizer = wx.BoxSizer(wx.HORIZONTAL)
    19091911        hSizer.Add(wx.StaticText(GeneralRestr,wx.ID_ANY,'Weight factor: '))
    1910         hSizer.Add(
    1911                     G2G.ValidatedTxtCtrl(GeneralRestr,generalRestData,
    1912                                 'wtFactor',nDig=(10,1),typeHint=float)
    1913                     )
     1912        hSizer.Add(G2G.ValidatedTxtCtrl(GeneralRestr,generalRestData,
     1913            'wtFactor',nDig=(10,1),typeHint=float))
    19141914        btn = G2G.G2CheckBox(GeneralRestr,'Use?',generalRestData,'Use')
    19151915        hSizer.Add(btn)

  • trunk/GSASIIsasd.py

    r4072 r4076  
    10781078    print (' Final chi^2: %.3f'%(chisq))
    10791079    data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
     1080   
     1081################################################################################
     1082#### Modelling
     1083################################################################################
     1084
     1085def PairDistFxn(Profile,ProfDict,Limits,Sample,data):
     1086    print('create P(R) fit to data - TBD')
     1087   
     1088    def CalcMoore():
     1089       
     1090#        def MooreRg(cw,q):      #lines 1400-1409
     1091#            n = 0
     1092#            nmax = len(cw)-5
     1093#            POR = 0.
     1094#            dmax = cw[nmax+2]
     1095#            while n < nmax:
     1096#                POR += cw[n]*((-1**n)*((np.pi*(n+1))**2-6)/((n+1)**3))
     1097#                n += 1
     1098#            POR *= 4*(dmax**4)/((np.pi**2)*cw[nmax+3])
     1099#            return POR
     1100       
     1101        def MoorePOR(cw,r,dmax):    #lines 1417-1428
     1102            n = 0
     1103            nmax = len(cw)
     1104            POR = np.zeros(len(r))
     1105            while n < nmax:
     1106                POR += 4.*r*cw[n]*np.sin((n+1.)*np.pi*r/dmax)
     1107                n += 1
     1108            return POR
     1109       
     1110        def MooreIOREFF(cw,q,dmax):      #lines 1437-1448
     1111            n = 0
     1112            nmax = len(cw)
     1113            POR = np.zeros(len(q))
     1114            dq = dmax*q
     1115            fpd = 8.*(np.pi**2)*dmax*np.sin(dq)/q
     1116            while n < nmax:
     1117                POR += cw[n]*(n+1.)*((-1)**n)*fpd/(((n+1.)*np.pi)**2-dq**2)
     1118                n += 1
     1119            return POR
     1120       
     1121#        def MooreINaught(cw,q,dmax):         #lines 1457-1466
     1122#            n = 0
     1123#            nmax = len(cw)
     1124#            POR = 0.
     1125#            while n < nmax:
     1126#                POR += cw[n]*8*(dmax**2)*((-1.)**n)/(n+1)
     1127#                n += 1
     1128#            return POR
     1129       
     1130        def calcSASD(values,Q,Io,wt,Ifb,dmax,ifBack):
     1131            if ifBack:
     1132                M = np.sqrt(wt)*(MooreIOREFF(values[:-1],Q,dmax)+Ifb+values[-1]-Io)
     1133            else:
     1134                M = np.sqrt(wt)*(MooreIOREFF(values,Q,dmax)+Ifb-Io)
     1135            return M
     1136
     1137        Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
     1138        Qmin = Limits[1][0]
     1139        Qmax = Limits[1][1]
     1140        wtFactor = ProfDict['wtFactor']
     1141        Back,ifBack = data['Back']
     1142        Ibeg = np.searchsorted(Q,Qmin)
     1143        Ifin = np.searchsorted(Q,Qmax)+1    #include last point
     1144        Ic[Ibeg:Ifin] = 0
     1145        Bins = np.linspace(0.,pairData['MaxRadius'],pairData['NBins']+1,True)
     1146        Dbins = np.diff(Bins)
     1147        Bins = Bins[:-1]+Dbins/2.
     1148        N = pairData['Moore']
     1149        if ifBack:
     1150            N += 1
     1151        MPV = np.zeros(N)
     1152        MPS = np.zeros(N)
     1153        MPV[0] = Q[Ibeg]
     1154        dmax = pairData['MaxRadius']
     1155        result = so.leastsq(calcSASD,MPV,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
     1156            args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],dmax,ifBack))
     1157        if ifBack:
     1158            MPVR = result[0][:-1]
     1159            data['Back'][0] = result[0][-1]
     1160            Back = data['Back'][0]
     1161        else:       
     1162            MPVR = result[0]
     1163            Back = 0.
     1164        chisq = np.sum(result[2]['fvec']**2)
     1165        Ic[Ibeg:Ifin] = MooreIOREFF(MPVR,Q[Ibeg:Ifin],dmax)+Ifb+Back
     1166        covM = result[1]
     1167        GOF = chisq/(Ifin-Ibeg-N)
     1168        MPS = np.sqrt(np.diag(covM)*GOF)
     1169        BinMag = MoorePOR(MPVR,Bins,dmax)/2.
     1170        return Bins,Dbins,BinMag
     1171   
     1172    def CalcRegular():
     1173        Bins = np.linspace(0.,pairData['MaxRadius'],pairData['NBins']+1,True)
     1174        Dbins = np.diff(Bins)
     1175        Bins = Bins[:-1]+Dbins/2.
     1176
     1177#test       
     1178        midBin = pairData['MaxRadius']/4.
     1179        wid = midBin/4.
     1180        BinMag  = 1./np.sqrt(2.*np.pi*wid*2)*np.exp(-(midBin-Bins)**2/(2.*wid**2))
     1181        return Bins,Dbins,BinMag
     1182               
     1183    pairData = data['Pair']
     1184   
     1185    if pairData['Method'] == 'Regularization':
     1186        print('Regularization calc; dummy Gaussian - TBD')
     1187        Bins,Dbins,BinMag = CalcRegular()
     1188       
     1189       
     1190    elif pairData['Method'] == 'Moore':
     1191        Bins,Dbins,BinMag = CalcMoore()       
     1192   
     1193    data['Pair']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
     1194   
     1195   
    10801196       
    10811197################################################################################
    … …  
    12921408    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
    12931409        return False,0,0,0,0,0,0,Msg
     1410
     1411def RgFit(Profile,ProfDict,Limits,Sample,Model):
     1412    print('unified fit single Rg to data to estimate a size')
     1413   
     1414    def GetModelParms():
     1415        parmDict = {}
     1416        varyList = []
     1417        values = []
     1418        Back = Model['Back']
     1419        if Back[1]:
     1420            varyList += ['Back',]
     1421            values.append(Back[0])
     1422        parmDict['Back'] = Back[0]
     1423        pairData = Model['Pair']
     1424        parmDict['G'] = pairData.get('Dist G',100.)
     1425        parmDict['Rg'] = pairData['MaxRadius']/2.5
     1426        varyList  += ['G','Rg',]
     1427        values += [parmDict['G'],parmDict['Rg'],]
     1428        return parmDict,varyList,values
     1429
     1430    def calcSASD(values,Q,Io,wt,Ifb,parmDict,varyList):
     1431        parmDict.update(zip(varyList,values))
     1432        M = np.sqrt(wt)*(getSASD(Q,parmDict)+Ifb-Io)
     1433        return M
     1434   
     1435    def getSASD(Q,parmDict):
     1436        Ic = np.zeros_like(Q)
     1437        Rg,G = parmDict['Rg'],parmDict['G']
     1438        Guin = G*np.exp(-(Q*Rg)**2/3)
     1439        Ic += Guin
     1440        Ic += parmDict['Back']  #/parmDict['Scale']
     1441        return Ic
     1442       
     1443    def SetModelParms():
     1444        print (' Refined parameters: ')
     1445        if 'Back' in varyList:
     1446            Model['Back'][0] = parmDict['Back']
     1447            print ('  %15s %15.4f esd: %15.4g'%('Background:',parmDict['Back'],sigDict['Back']))
     1448        pairData = Model['Pair']
     1449        pairData['Dist G'] = parmDict['G']
     1450        pairData['MaxRadius'] = parmDict['Rg']*2.5
     1451        for item in varyList:
     1452            print (' %15s: %15.4g esd: %15.4g'%(item,parmDict[item],sigDict[item]))
     1453
     1454    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
     1455    Qmin = Limits[1][0]
     1456    Qmax = Limits[1][1]
     1457    wtFactor = ProfDict['wtFactor']
     1458    Ibeg = np.searchsorted(Q,Qmin)
     1459    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
     1460    Ic[:] = 0
     1461    parmDict,varyList,values = GetModelParms()
     1462    result = so.leastsq(calcSASD,values,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
     1463        args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],parmDict,varyList))
     1464    parmDict.update(zip(varyList,result[0]))
     1465    chisq = np.sum(result[2]['fvec']**2)
     1466    ncalc = result[2]['nfev']
     1467    covM = result[1]
     1468    Rvals = {}
     1469    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
     1470    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
     1471    Ic[Ibeg:Ifin] = getSASD(Q[Ibeg:Ifin],parmDict)
     1472    Msg = 'Failed to converge'
     1473    try:
     1474        Nans = np.isnan(result[0])
     1475        if np.any(Nans):
     1476            name = varyList[Nans.nonzero(True)[0]]
     1477            Msg = 'Nan result for '+name+'!'
     1478            raise ValueError
     1479        Negs = np.less_equal(result[0],0.)
     1480        if np.any(Negs):
     1481            name = varyList[Negs.nonzero(True)[0]]
     1482            Msg = 'negative coefficient for '+name+'!'
     1483            raise ValueError
     1484        if len(covM):
     1485            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
     1486            sigDict = dict(zip(varyList,sig))
     1487        print (' Results of Rg fit:')
     1488        print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
     1489        print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
     1490        SetModelParms()
     1491        covMatrix = covM*Rvals['GOF']
     1492        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
     1493    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
     1494        return False,0,0,0,0,0,0,Msg
     1495
     1496    return [None,]
     1497
    12941498   
    12951499def ModelFxn(Profile,ProfDict,Limits,Sample,sasdData):

  • trunk/imports/G2phase.py

    r4073 r4076  
    7575        file = open(filename, 'Ur')
    7676        Phase = {}
    77         Title = ''
     77        Title = os.path.basename(filename)
     78       
    7879        Compnd = ''
    7980        Atoms = []
    … …  
    128129                    self.warnings += "Change this in phase's General tab."
    129130                    SGData = G2obj.P1SGData # P 1
    130                     cell = [1.0,1.0,1.0,90.,90.,90.]
     131                    cell = [20.0,20.0,20.0,90.,90.,90.]
    131132                    Volume = G2lat.calc_V(G2lat.cell2A(cell))
    132133                    AA,AB = G2lat.cell2AB(cell)                   
Note: See TracChangeset for help on using the changeset viewer.