Changeset 1277 for trunk/GSASIIsasd.py


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

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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.