Changeset 1277


Ignore:
Timestamp:
Apr 16, 2014 12:09:16 PM (10 years ago)
Author:
vondreele
Message:

1st working SASD fit version. Implement UnDo? option for SASD fitting.

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIgrid.py

    r1269 r1277  
    137137[ wxID_MODELCOPY,wxID_MODELFIT,wxID_MODELADD,wxID_ELEMENTADD,wxID_ELEMENTDELETE,
    138138    wxID_ADDSUBSTANCE,wxID_LOADSUBSTANCE,wxID_DELETESUBSTANCE,wxID_COPYSUBSTANCE,
    139 ] = [wx.NewId() for item in range(9)]
     139    wxID_MODELUNDO,
     140] = [wx.NewId() for item in range(10)]
    140141
    141142[ wxID_SELECTPHASE,wxID_PWDHKLPLOT,
    … …  
    26142615        self.ModelEdit.Append(id=wxID_MODELFIT, kind=wx.ITEM_NORMAL,text='Fit',
    26152616            help='Fit model parameters to data')
     2617        self.SasdUndo = self.ModelEdit.Append(id=wxID_MODELUNDO, kind=wx.ITEM_NORMAL,text='Undo',
     2618            help='Undo model fit')
     2619        self.SasdUndo.Enable(False)           
    26162620        self.ModelEdit.Append(id=wxID_MODELCOPY, kind=wx.ITEM_NORMAL,text='Copy',
    26172621            help='Copy model parameters to other histograms')

  • trunk/GSASIIpwd.py

    r1274 r1277  
    11181118def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,oneCycle=False,controls=None,dlg=None):
    11191119    'needs a doc string'
    1120    
    11211120       
    11221121    def GetBackgroundParms(parmList,Background):

  • trunk/GSASIIpwdGUI.py

    r1274 r1277  
    14001400        G2sasd.SetScale(Data,refData)
    14011401        UpdateSampleGrid(G2frame,data)       
    1402         G2plt.PlotPatterns(G2frame,plotType='SASD',newPlot=True)
     1402        G2plt.PlotPatterns(G2frame,plotType='SASD',newPlot=False)
    14031403       
    14041404    def OnSampleCopy(event):
    … …  
    26642664    #end patches
    26652665   
    2666     def OnCopyModel(event):
    2667         print 'copy model'
    2668         print data
    2669        
    26702666    def OnAddModel(event):
    26712667        if data['Current'] == 'Particle fit':
    … …  
    26772673                    'FFargs':{},'NumPoints':50,'Cutoff':0.01,'AutoDist':True,'logR':True,
    26782674                    'StrFact':'Dilute'},
    2679                 'LogNormal':{'Volume':[0.05,False],'MinSize':[10.,False],'Mean':[1000.,False],'StdDev':[0.5,False]},
     2675                'LogNormal':{'Volume':[0.05,False],'Mean':[1000.,False],'StdDev':[0.5,False],'MinSize':[10.,False],},
    26802676                'Gaussian':{'Volume':[0.05,False],'Mean':[1000.,False],'StdDev':[300.,False],},
    26812677                'LSW':{'Volume':[0.05,False],'Mean':[1000.0,False],},
    … …  
    26922688        wx.CallAfter(UpdateModelsGrid,G2frame,data)
    26932689       
     2690    def OnCopyModel(event):
     2691        print 'copy model'
     2692        print data
     2693       
    26942694    def OnFitModel(event):
    2695         print 'fit model for '+data['Current']
    26962695        if not any(Sample['Contrast']):
    26972696            G2frame.ErrorDialog('No contrast; your sample is a vacuum!',
    … …  
    27052704           
    27062705        elif data['Current'] == 'Particle fit':
     2706            SaveState()
    27072707            G2sasd.ModelFit(Profile,ProfDict,Limits,Substances,Sample,data)
    27082708            G2plt.PlotPatterns(G2frame,plotType='SASD',newPlot=True)
     2709            wx.CallAfter(UpdateModelsGrid,G2frame,data)
    27092710           
     2711    def OnUnDo(event):
     2712        DoUnDo()
     2713        data = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,
     2714            G2frame.PatternId,'Models'))
     2715        G2frame.dataFrame.SasdUndo.Enable(False)
     2716        UpdateModelsGrid(G2frame,data)
     2717        G2sasd.ModelFxn(Profile,ProfDict,Limits,Substances,Sample,data)
     2718        G2plt.PlotPatterns(G2frame,plotType='SASD',newPlot=True)
     2719
     2720    def DoUnDo():
     2721        print 'Undo last refinement'
     2722        file = open(G2frame.undosasd,'rb')
     2723        PatternId = G2frame.PatternId
     2724        for item in ['Models']:
     2725            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, item),cPickle.load(file))
     2726            if G2frame.dataDisplay.GetName() == item:
     2727                if item == 'Background':
     2728                    UpdateBackground(G2frame,G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, item)))
     2729                elif item == 'Instrument Parameters':
     2730                    UpdateInstrumentGrid(G2frame,G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, item)))
     2731                elif item == 'Peak List':
     2732                    UpdatePeakGrid(G2frame,G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, item)))
     2733            print item,' recovered'
     2734        file.close()
     2735       
     2736    def SaveState():
     2737        G2frame.undosasd = os.path.join(G2frame.dirname,'GSASIIsasd.save')
     2738        file = open(G2frame.undosasd,'wb')
     2739        PatternId = G2frame.PatternId
     2740        for item in ['Models']:
     2741            cPickle.dump(G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,item)),file,1)
     2742        file.close()
     2743        G2frame.dataFrame.SasdUndo.Enable(True)
     2744       
    27102745    def OnSelectFit(event):
    27112746        data['Current'] = fitSel.GetValue()
    … …  
    27232758        if itemKey == 'Back':
    27242759            Profile[4][:] = value
    2725         G2plt.PlotPatterns(G2frame,plotType='SASD',newPlot=True)
     2760        G2plt.PlotPatterns(G2frame,plotType='SASD',newPlot=False)
    27262761       
    27272762    def OnCheckBox(event):
    … …  
    28672902            PlotText = G2frame.G2plotNB.nb.GetPageText(G2frame.G2plotNB.nb.GetSelection())
    28682903            if 'Powder' in PlotText:
    2869                 G2plt.PlotPatterns(G2frame,plotType='SASD',newPlot=True)
     2904                G2plt.PlotPatterns(G2frame,plotType='SASD',newPlot=False)
    28702905            elif 'Size' in PlotText:
    28712906                G2plt.PlotSASDSizeDist(G2frame)
    … …  
    28732908        def OnValue(event):
    28742909            Obj = event.GetEventObject()
    2875             item,parm,sldrObj = Indx[Obj.GetId()]
     2910            item,key,id,sldrObj = Indx[Obj.GetId()]
    28762911            try:
    28772912                value = float(Obj.GetValue())
    … …  
    28792914                    raise ValueError
    28802915            except ValueError:
    2881                 value = item[0]
    2882             item[0] = value
     2916                value = item[key][id]
     2917            item[key][id] = value
    28832918            Obj.SetValue('%.3g'%(value))
    2884             sldrObj.SetRange(1000.*(np.log10(value)-2),1000.*(np.log10(value)+2))
    2885             sldrObj.SetValue(1000.*np.log10(value))
     2919            if key == 'P':
     2920                sldrObj.SetValue(1000.*value)
     2921            else:
     2922                sldrObj.SetRange(1000.*(np.log10(value)-2),1000.*(np.log10(value)+2))
     2923                sldrObj.SetValue(1000.*np.log10(value))
    28862924            G2sasd.ModelFxn(Profile,ProfDict,Limits,Substances,Sample,data)
    28872925            RefreshPlots()
    … …  
    29082946        def OnParmSlider(event):
    29092947            Obj = event.GetEventObject()
    2910             item,key,pvObj = Indx[Obj.GetId()]
     2948            item,key,id,pvObj = Indx[Obj.GetId()]
    29112949            slide = Obj.GetValue()
    2912             value = 10.**float(slide/1000.)
    2913             item[key] = value
    2914             pvObj.SetValue('%.3g'%(item[key]))
     2950            if key == 'P':
     2951                value = float(slide/1000.)
     2952            else:
     2953                value = 10.**float(slide/1000.)
     2954            item[key][id] = value
     2955            pvObj.SetValue('%.3g'%(item[key][id]))
    29152956            G2sasd.ModelFxn(Profile,ProfDict,Limits,Substances,Sample,data)
    29162957            RefreshPlots()
    … …  
    30053046            Parms = level[level['Controls']['DistType']]
    30063047            FFargs = level['Controls']['FFargs']
    3007             for iparm,parm in enumerate(list(Parms)):
    3008                 parmVar = wx.CheckBox(G2frame.dataDisplay,label='Refine? Dist '+parm)
    3009                 parmVar.SetValue(Parms[parm][1])
    3010                 parmVar.Bind(wx.EVT_CHECKBOX, OnSelect)
    3011                 parmSizer.Add(parmVar,0,WACV)
    3012                 parmValue = wx.TextCtrl(G2frame.dataDisplay,value='%.3g'%(Parms[parm][0]),
    3013                     style=wx.TE_PROCESS_ENTER)
    3014                 parmValue.Bind(wx.EVT_TEXT_ENTER,OnValue)       
    3015                 parmValue.Bind(wx.EVT_KILL_FOCUS,OnValue)
    3016                 parmSizer.Add(parmValue,0,WACV)
    3017                 value = np.log10(Parms[parm][0])
    3018                 valMinMax = [value-1,value+1]
    3019                 parmSldr = wx.Slider(G2frame.dataDisplay,minValue=1000.*valMinMax[0],
    3020                     maxValue=1000.*valMinMax[1],value=1000.*value)
    3021                 Indx[parmVar.GetId()] = [Parms[parm],parm]
    3022                 Indx[parmValue.GetId()] = [Parms[parm],0,parmSldr]
    3023                 Indx[parmSldr.GetId()] = [Parms[parm],0,parmValue]
    3024                 parmSldr.Bind(wx.EVT_SLIDER,OnParmSlider)
    3025                 parmSizer.Add(parmSldr,1,wx.EXPAND)
    3026                 center = wx.CheckBox(G2frame.dataDisplay,label='Center? ')
    3027                 Indx[center.GetId()] = [parmSldr,value]
    3028                 center.Bind(wx.EVT_CHECKBOX,OnCenterSlider)
    3029                 parmSizer.Add(center,0,WACV)
    3030             for parm in list(FFargs):
    3031                 parmVar = wx.CheckBox(G2frame.dataDisplay,label='Refine? FF '+parm)
    3032                 parmVar.SetValue(FFargs[parm][1])
    3033                 Indx[parmVar.GetId()] = [FFargs[parm],1]
    3034                 parmVar.Bind(wx.EVT_CHECKBOX, OnSelect)
    3035                 parmSizer.Add(parmVar,0,WACV)
    3036                 parmValue = wx.TextCtrl(G2frame.dataDisplay,value='%.3g'%(FFargs[parm][0]),
    3037                     style=wx.TE_PROCESS_ENTER)
    3038                 parmValue.Bind(wx.EVT_TEXT_ENTER,OnValue)       
    3039                 parmValue.Bind(wx.EVT_KILL_FOCUS,OnValue)
    3040                 parmSizer.Add(parmValue,0,WACV)
    3041                 value = np.log10(FFargs[parm][0])
    3042                 valMinMax = [value-1,value+1]
    3043                 parmSldr = wx.Slider(G2frame.dataDisplay,minValue=1000.*valMinMax[0],
    3044                     maxValue=1000.*valMinMax[1],value=1000.*value)
    3045                 Indx[parmVar.GetId()] = [FFargs[parm],parm]
    3046                 Indx[parmValue.GetId()] = [FFargs[parm],0,parmSldr]
    3047                 Indx[parmSldr.GetId()] = [FFargs[parm],0,parmValue]
    3048                 parmSldr.Bind(wx.EVT_SLIDER,OnParmSlider)
    3049                 parmSizer.Add(parmSldr,1,wx.EXPAND)
    3050                 center = wx.CheckBox(G2frame.dataDisplay,label='Center? ')
    3051                 Indx[center.GetId()] = [parmSldr,value]
    3052                 center.Bind(wx.EVT_CHECKBOX,OnCenterSlider)
    3053                 parmSizer.Add(center,0,WACV)
     3048            parmOrder = ['Volume','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
     3049                'PkInt','PkPos','PkSig','PkGam',]
     3050            for parm in parmOrder:
     3051                if parm in Parms:
     3052                    parmVar = wx.CheckBox(G2frame.dataDisplay,label='Refine? Dist '+parm)
     3053                    parmVar.SetValue(Parms[parm][1])
     3054                    parmVar.Bind(wx.EVT_CHECKBOX, OnSelect)
     3055                    parmSizer.Add(parmVar,0,WACV)
     3056                    parmValue = wx.TextCtrl(G2frame.dataDisplay,value='%.3g'%(Parms[parm][0]),
     3057                        style=wx.TE_PROCESS_ENTER)
     3058                    parmValue.Bind(wx.EVT_TEXT_ENTER,OnValue)       
     3059                    parmValue.Bind(wx.EVT_KILL_FOCUS,OnValue)
     3060                    parmSizer.Add(parmValue,0,WACV)
     3061                    if parm == 'P':
     3062                        value = Parms[parm][0]
     3063                        valMinMax = [0.1,4.2]
     3064                    else:
     3065                        value = np.log10(Parms[parm][0])
     3066                        valMinMax = [value-1,value+1]
     3067                    parmSldr = wx.Slider(G2frame.dataDisplay,minValue=1000.*valMinMax[0],
     3068                        maxValue=1000.*valMinMax[1],value=1000.*value)
     3069                    Indx[parmVar.GetId()] = [Parms[parm],1]
     3070                    Indx[parmValue.GetId()] = [Parms,parm,0,parmSldr]
     3071                    Indx[parmSldr.GetId()] = [Parms,parm,0,parmValue]
     3072                    parmSldr.Bind(wx.EVT_SLIDER,OnParmSlider)
     3073                    parmSizer.Add(parmSldr,1,wx.EXPAND)
     3074                    if parm == 'P':
     3075                        parmSizer.Add((5,5),)
     3076                    else:
     3077                        center = wx.CheckBox(G2frame.dataDisplay,label='Center? ')
     3078                        Indx[center.GetId()] = [parmSldr,value]
     3079                        center.Bind(wx.EVT_CHECKBOX,OnCenterSlider)
     3080                        parmSizer.Add(center,0,WACV)
     3081            if level['Controls']['DistType'] not in ['Bragg']:
     3082                for parm in list(FFargs):
     3083                    parmVar = wx.CheckBox(G2frame.dataDisplay,label='Refine? FF '+parm)
     3084                    parmVar.SetValue(FFargs[parm][1])
     3085                    Indx[parmVar.GetId()] = [FFargs[parm],1]
     3086                    parmVar.Bind(wx.EVT_CHECKBOX, OnSelect)
     3087                    parmSizer.Add(parmVar,0,WACV)
     3088                    parmValue = wx.TextCtrl(G2frame.dataDisplay,value='%.3g'%(FFargs[parm][0]),
     3089                        style=wx.TE_PROCESS_ENTER)
     3090                    parmValue.Bind(wx.EVT_TEXT_ENTER,OnValue)       
     3091                    parmValue.Bind(wx.EVT_KILL_FOCUS,OnValue)
     3092                    parmSizer.Add(parmValue,0,WACV)
     3093                    value = np.log10(FFargs[parm][0])
     3094                    valMinMax = [value-1,value+1]
     3095                    parmSldr = wx.Slider(G2frame.dataDisplay,minValue=1000.*valMinMax[0],
     3096                        maxValue=1000.*valMinMax[1],value=1000.*value)
     3097                    Indx[parmVar.GetId()] = [FFargs[parm],1]
     3098                    Indx[parmValue.GetId()] = [FFargs[parm],0,parmSldr]
     3099                    Indx[parmSldr.GetId()] = [FFargs,parm,0,parmValue]
     3100                    parmSldr.Bind(wx.EVT_SLIDER,OnParmSlider)
     3101                    parmSizer.Add(parmSldr,1,wx.EXPAND)
     3102                    center = wx.CheckBox(G2frame.dataDisplay,label='Center? ')
     3103                    Indx[center.GetId()] = [parmSldr,value]
     3104                    center.Bind(wx.EVT_CHECKBOX,OnCenterSlider)
     3105                    parmSizer.Add(center,0,WACV)
    30543106            lvlSizer.Add(parmSizer,1,wx.EXPAND)
    30553107            partSizer.Add(lvlSizer,1,wx.EXPAND)
    … …  
    30823134    G2frame.dataFrame.Bind(wx.EVT_MENU, OnCopyModel, id=G2gd.wxID_MODELCOPY)
    30833135    G2frame.dataFrame.Bind(wx.EVT_MENU, OnFitModel, id=G2gd.wxID_MODELFIT)
     3136    G2frame.dataFrame.Bind(wx.EVT_MENU, OnUnDo, id=G2gd.wxID_MODELUNDO)
    30843137    G2frame.dataFrame.Bind(wx.EVT_MENU, OnAddModel, id=G2gd.wxID_MODELADD)
    30853138    Indx = {}

  • trunk/GSASIIsasd.py

    r1274 r1277  
    912912################################################################################
    913913
    914 def ModelFit(Profile,ProfDict,Limits,Substances,Sample,data):
    915     print 'do model fit'
     914def ModelFit(Profile,ProfDict,Limits,Substances,Sample,Model):
     915    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
     916        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
     917        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
     918        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
     919        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
     920
     921    def GetModelParms():
     922        parmDict = {}
     923        varyList = []
     924        values = []
     925        levelTypes = []
     926        Back = Model['Back']
     927        if Back[1]:
     928            varyList += ['Back',]
     929            values.append(Back[0])
     930        parmDict['Back'] = Back[0]
     931        partData = Model['Particle']
     932        parmDict['Matrix density'] = Substances['Substances'][partData['Matrix']['Name']].get('XAnom density',0.0)
     933        for i,level in enumerate(partData['Levels']):
     934            cid = str(i)+':'
     935            controls = level['Controls']
     936            Type = controls['DistType']
     937            levelTypes.append(Type)
     938            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
     939                if Type not in ['Monodosperse',]:
     940                    parmDict[cid+'NumPoints'] = controls['NumPoints']
     941                    parmDict[cid+'Cutoff'] = controls['Cutoff']
     942                parmDict[cid+'FormFact'] = shapes[controls['FormFact']][0]
     943                parmDict[cid+'FFVolume'] = shapes[controls['FormFact']][1]
     944                parmDict[cid+'XAnom density'] = Substances['Substances'][controls['Material']].get('XAnom density',0.0)
     945                for item in ['Aspect ratio','Length','Thickness','Diameter',]:
     946                    if item in controls['FFargs']:
     947                        parmDict[cid+item] = controls['FFargs'][item][0]
     948                        if controls['FFargs'][item][1]:
     949                            varyList.append(cid+item)
     950                            values.append(controls['FFargs'][item][0])
     951            distDict = controls['DistType']
     952            for item in level[distDict]:
     953                parmDict[cid+item] = level[distDict][item][0]
     954                if level[distDict][item][1]:
     955                    values.append(level[distDict][item][0])
     956                    varyList.append(cid+item)
     957        return levelTypes,parmDict,varyList,values
     958       
     959    def SetModelParms():
     960        print ' Refined parameters:'
     961        if 'Back' in varyList:
     962            Model['Back'][0] = parmDict['Back']
     963            print '  %15s %15.4f esd: %15.4g'%('Background:',parmDict['Back'],sigDict['Back'])
     964        partData = Model['Particle']
     965        for i,level in enumerate(partData['Levels']):
     966            Type = level['Controls']['DistType']
     967            print ' Component %d: Type: %s:'%(i,Type)
     968            cid = str(i)+':'
     969            controls = level['Controls']
     970            Type = controls['DistType']
     971            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
     972                for item in ['Aspect ratio','Length','Thickness','Diameter',]:
     973                    if cid+item in varyList:
     974                        controls['FFargs'][item][0] = parmDict[cid+item]
     975                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
     976            distDict = controls['DistType']
     977            for item in level[distDict]:
     978                if cid+item in varyList:
     979                    level[distDict][item][0] = parmDict[cid+item]
     980                    print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
     981                   
     982    def calcSASD(values,Q,Io,wt,levelTypes,parmDict,varyList):
     983        parmDict.update(zip(varyList,values))
     984        M = np.sqrt(wt)*(getSASD(Q,levelTypes,parmDict)-Io)
     985        return M
     986       
     987    def getSASD(Q,levelTypes,parmDict):
     988        Ic = np.ones_like(Q)
     989        Ic *= parmDict['Back']
     990        rhoMat = parmDict['Matrix density']
     991        for i,Type in enumerate(levelTypes):
     992            cid = str(i)+':'
     993            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
     994                FFfxn = parmDict[cid+'FormFact']
     995                Volfxn = parmDict[cid+'FFVolume']
     996                FFargs = []
     997                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
     998                    if item in parmDict:
     999                        FFargs.append(parmDict[item])
     1000                distDict = {}
     1001                for item in [cid+'Volume',cid+'Mean',cid+'StdDev',cid+'MinSize',]:
     1002                    if item in parmDict:
     1003                        distDict[item.split(':')[1]] = parmDict[item]
     1004                rho = parmDict[cid+'XAnom density']
     1005                contrast = rho**2-rhoMat**2
     1006                rBins,dBins,dist = MakeDiamDist(Type,parmDict[cid+'NumPoints'],parmDict[cid+'Cutoff'],distDict)
     1007                Gmat = G_matrix(Q,rBins,contrast,FFfxn,Volfxn,FFargs).T
     1008                dist *= parmDict[cid+'Volume']
     1009                Ic += np.dot(Gmat,dist)
     1010            elif 'Unified' in Type:
     1011                Rg,G,B,P = parmDict[cid+'Rg'],parmDict[cid+'G'],parmDict[cid+'B'],parmDict[cid+'P']
     1012                Qstar = Q/(scsp.erf(Q*Rg/np.sqrt(6)))**3
     1013                Guin = G*np.exp(-(Q*Rg)**2/3)
     1014                Porod = (B/Qstar**P)
     1015                Ic += Guin+Porod
     1016            elif 'Porod' in Type:
     1017                B,P,Rgco = parmDict[cid+'B'],parmDict[cid+'P'],parmDict[cid+'Cutoff']
     1018                Porod = (B/Q**P)*np.exp(-(Q*Rgco)**2/3)
     1019                Ic += Porod
     1020            elif 'Mono' in Type:
     1021                FFfxn = parmDict[cid+'FormFact']
     1022                Volfxn = parmDict[cid+'FormFact']
     1023                FFargs = []
     1024                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
     1025                    if item in parmDict:
     1026                        FFargs.append(parmDict[item])
     1027                rho = parmDict[cid+'XAnom density']
     1028                contrast = rho**2-rhoMat**2
     1029                R = parmDict[cid+'Radius']
     1030                Gmat = G_matrix(Q,R,contrast,FFfxn,Volfxn,FFargs)             
     1031                Ic += Gmat[0]*parmDict[cid+'Volume']
     1032            elif 'Bragg' in distFxn:
     1033                Ic += parmDict[cid+'PkInt']*G2pwd.getPsVoigt(parmDict[cid+'PkPos'],
     1034                    parmDict[cid+'PkSig'],parmDict[cid+'PkGam'],Q)
     1035        return Ic
     1036       
     1037    Q,Io,wt,Ic,Ib = Profile[:5]
     1038    Qmin = Limits[1][0]
     1039    Qmax = Limits[1][1]
     1040    wtFactor = ProfDict['wtFactor']
     1041    Ibeg = np.searchsorted(Q,Qmin)
     1042    Ifin = np.searchsorted(Q,Qmax)
     1043    Ic[:] = 0
     1044    levelTypes,parmDict,varyList,values = GetModelParms()
     1045    result = so.leastsq(calcSASD,values,full_output=True,   #ftol=Ftol,
     1046        args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],levelTypes,parmDict,varyList))
     1047    ncyc = int(result[2]['nfev']/2)
     1048    chisq = np.sum(result[2]['fvec']**2)
     1049    Rwp = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
     1050    GOF = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
     1051    parmDict.update(zip(varyList,result[0]))
     1052    Ic[Ibeg:Ifin] = getSASD(Q[Ibeg:Ifin],levelTypes,parmDict)
     1053    sigDict = dict(zip(varyList,np.sqrt(np.diag(result[1])*GOF)))
     1054    print ' Results of small angle data modelling fit:'
     1055    print 'Number of function calls:',result[2]['nfev'],' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
     1056    print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
     1057    SetModelParms()
    9161058   
    9171059def ModelFxn(Profile,ProfDict,Limits,Substances,Sample,sasdData):
    … …  
    9511093            contrast = rho**2-rhoMat**2
    9521094            parmDict = level[controls['DistType']]
    953             rBins,dBins,dist = MakeDiamDist(controls['DistType'],controls['NumPoints'],controls['Cutoff'],parmDict)
     1095            distDict = {}
     1096            for item in parmDict:
     1097                distDict[item] = parmDict[item][0]
     1098            rBins,dBins,dist = MakeDiamDist(controls['DistType'],controls['NumPoints'],controls['Cutoff'],distDict)
    9541099            Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs).T
    9551100            dist *= level[distFxn]['Volume'][0]
    … …  
    9581103            Dist.append(dist/(4.*dBins))
    9591104        elif 'Unified' in distFxn:
    960             rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
    9611105            parmDict = level[controls['DistType']]
    9621106            Rg,G,B,P = parmDict['Rg'][0],parmDict['G'][0],parmDict['B'][0],parmDict['P'][0]
    … …  
    9681112            Dist.append([])
    9691113        elif 'Porod' in distFxn:
    970             rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
    9711114            parmDict = level[controls['DistType']]
    9721115            B,P,Rgco = parmDict['B'][0],parmDict['P'][0],parmDict['Cutoff'][0]
    … …  
    9941137                parmDict['PkSig'][0],parmDict['PkGam'][0],Q[Ibeg:Ifin])
    9951138            Rbins.append([])
    996             Dist.append([])
    997            
     1139            Dist.append([])           
    9981140    sasdData['Size Calc'] = [Rbins,Dist]
    9991141   
    1000 def MakeDiamDist(DistName,nPoints,cutoff,parmDict):
     1142def MakeDiamDist(DistName,nPoints,cutoff,distDict):
    10011143   
    10021144    if 'LogNormal' in DistName:
    10031145        distFxn = 'LogNormalDist'
    10041146        cumeFxn = 'LogNormalCume'
    1005         pos = parmDict['MinSize'][0]
    1006         args = [parmDict['Mean'][0],parmDict['StdDev'][0]]
    1007         step = 4.*np.sqrt(np.exp(parmDict['StdDev'][0]**2)*(np.exp(parmDict['StdDev'][0]**2)-1.))
    1008         mode = parmDict['MinSize'][0]+parmDict['Mean'][0]/np.exp(parmDict['StdDev'][0]**2)
     1147        pos = distDict['MinSize']
     1148        args = [distDict['Mean'],distDict['StdDev']]
     1149        step = 4.*np.sqrt(np.exp(distDict['StdDev']**2)*(np.exp(distDict['StdDev']**2)-1.))
     1150        mode = distDict['MinSize']+distDict['Mean']/np.exp(distDict['StdDev']**2)
    10091151        minX = 1. #pos
    10101152        maxX = 1.e4 #np.min([mode+nPoints*step*40,1.e5])
    … …  
    10121154        distFxn = 'GaussDist'
    10131155        cumeFxn = 'GaussCume'
    1014         pos = parmDict['Mean'][0]
    1015         args = [parmDict['StdDev'][0]]
    1016         step = 0.02*parmDict['StdDev'][0]
    1017         mode = parmDict['Mean'][0]
    1018         minX = np.max([mode-4.*parmDict['StdDev'][0],1.])
    1019         maxX = np.min([mode+4.*parmDict['StdDev'][0],1.e5])
     1156        pos = distDict['Mean']
     1157        args = [distDict['StdDev']]
     1158        step = 0.02*distDict['StdDev']
     1159        mode = distDict['Mean']
     1160        minX = np.max([mode-4.*distDict['StdDev'],1.])
     1161        maxX = np.min([mode+4.*distDict['StdDev'],1.e5])
    10201162    elif 'LSW' in DistName:
    10211163        distFxn = 'LSWDist'
    10221164        cumeFxn = 'LSWCume'
    1023         pos = parmDict['Mean'][0]
     1165        pos = distDict['Mean']
    10241166        args = []
    10251167        minX,maxX = [0.,2.*pos]
    … …  
    10271169        distFxn = 'SchulzZimmDist'
    10281170        cumeFxn = 'SchulzZimmCume'
    1029         pos = parmDict['Mean'][0]
    1030         args = [parmDict['StdDev'][0]]
    1031         minX = np.max([1.,pos-4.*parmDict['StdDev'][0]])
    1032         maxX = np.min([pos+4.*parmDict['StdDev'][0],1.e5])
     1171        pos = distDict['Mean']
     1172        args = [distDict['StdDev']]
     1173        minX = np.max([1.,pos-4.*distDict['StdDev']])
     1174        maxX = np.min([pos+4.*distDict['StdDev'],1.e5])
    10331175    nP = 500
    10341176    Diam = np.logspace(0.,5.,nP,True)
Note: See TracChangeset for help on using the changeset viewer.