Changeset 2757


Ignore:
Timestamp:
Mar 21, 2017 11:01:51 AM (6 years ago)
Author:
vondreele
Message:

workable reflectometry fitting by Levenberg-Marquart least squares.

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIpwd.py

    r2755 r2757  
    19131913################################################################################
    19141914
    1915 def REFDModelFxn(Profile,ProfDict,Inst,Limits,Substances,data):
     1915def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
     1916    print 'fit REFD data'
     1917   
     1918    def GetModelParms():
     1919        parmDict = {}
     1920        varyList = []
     1921        values = []
     1922        for parm in ['Scale','FltBack']:
     1923            parmDict[parm] = data[parm][0]
     1924            if data[parm][1]:
     1925                varyList.append(parm)
     1926                values.append(data[parm][0])
     1927        parmDict['nLayers'] = len(data['Layers'])
     1928        for ilay,layer in enumerate(data['Layers']):
     1929            name = layer['Name']
     1930            cid = str(ilay)+';'
     1931            for parm in ['Thick','Rough','DenMul']:
     1932                parmDict[cid+parm] = layer.get(parm,[0.,False])[0]
     1933                if layer.get(parm,[0.,False])[1]:
     1934                    varyList.append(cid+parm)
     1935                    values.append(layer[parm][0])
     1936            parmDict[cid+'rho'] = Substances[name]['Scatt density']
     1937            parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
     1938        return parmDict,varyList,values
     1939   
     1940    def SetModelParms():
     1941        line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale'])
     1942        if 'Scale' in varyList:
     1943            data['Scale'] = parmDict['Scale']
     1944            line += ' esd: %.4g'%(sigDict['Scale'])                                                             
     1945        print line
     1946        line = ' Flat background: %15.4g'%(parmDict['FltBack'])
     1947        if 'FltBack' in varyList:
     1948            data['FltBack'] = parmDict['FltBack']
     1949            line += ' esd: %15.3g'%(sigDict['FltBack'])
     1950        print line
     1951        for ilay,layer in enumerate(data['Layers']):
     1952            name = layer['Name']
     1953            print ' Parameters for layer: %d %s'%(ilay,name)
     1954            cid = str(ilay)+';'
     1955            line = ' '
     1956            line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul'])
     1957            line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)**parmDict[cid+'DenMul'])
     1958            for parm in ['Thick','Rough','DenMul']:
     1959                if parm in layer:
     1960                    layer[parm][0] = parmDict[cid+parm]
     1961                    line += ' %s: %.3f'%(parm,layer[parm][0])
     1962                    if cid+parm in varyList:
     1963                        line += ' esd: %.3g'%(sigDict[cid+parm])
     1964            print line
     1965            print line2
     1966   
     1967    def calcREFD(values,Q,Io,wt,parmDict,varyList):
     1968        parmDict.update(zip(varyList,values))
     1969        M = np.sqrt(wt)*(getREFD(Q,parmDict)-Io)
     1970        return M
     1971   
     1972    def getREFD(Q,parmDict):
     1973        Ic = np.ones_like(Q)*parmDict['FltBack']
     1974        Scale = parmDict['Scale']
     1975        Nlayers = parmDict['nLayers']
     1976        depth = np.zeros(Nlayers)
     1977        rho = np.zeros(Nlayers)
     1978        irho = np.zeros(Nlayers)
     1979        sigma = np.zeros(Nlayers)
     1980        for ilay in range(Nlayers):
     1981            cid = str(ilay)+';'
     1982            depth[ilay] = parmDict[cid+'Thick']
     1983            sigma[ilay] = parmDict[cid+'Rough']
     1984            rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul']
     1985            irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul']
     1986            A,B = abeles(0.5*Q,depth,rho,irho,sigma[1:])     #Q --> k, offset roughness for abeles
     1987        Ic += (A**2+B**2)*Scale     
     1988        return Ic
     1989
     1990    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
     1991    Qmin = Limits[1][0]
     1992    Qmax = Limits[1][1]
     1993    wtFactor = ProfDict['wtFactor']
     1994    Ibeg = np.searchsorted(Q,Qmin)
     1995    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
     1996    Ic[:] = 0
     1997    parmDict,varyList,values = GetModelParms()
     1998    if varyList:
     1999        result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
     2000            args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList))
     2001        parmDict.update(zip(varyList,result[0]))
     2002        chisq = np.sum(result[2]['fvec']**2)
     2003        ncalc = result[2]['nfev']
     2004        covM = result[1]
     2005    else:   #nothing varied
     2006        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)
     2007        chisq = np.sum(M**2)
     2008        ncalc = 0
     2009        covM = []
     2010        sig = []
     2011        sigDict = {}
     2012        result = []
     2013    Rvals = {}
     2014    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
     2015    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
     2016    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],parmDict)
     2017    Ib[Ibeg:Ifin] = parmDict['FltBack']
     2018    Msg = 'Failed to converge'
     2019    try:
     2020        Nans = np.isnan(result[0])
     2021        if np.any(Nans):
     2022            name = varyList[Nans.nonzero(True)[0]]
     2023            Msg += ' Nan result for '+name+'!'
     2024            raise ValueError
     2025        Negs = np.less_equal(result[0],0.)
     2026        if np.any(Negs):
     2027            indx = Negs.nonzero()
     2028            name = varyList[indx[0][0]]
     2029            if name != 'FltBack':
     2030                Msg += ' negative coefficient for '+name+'!'
     2031                raise ValueError
     2032        if len(covM):
     2033            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
     2034            sigDict = dict(zip(varyList,sig))
     2035        print ' Results of reflectometry data modelling fit:'
     2036        print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
     2037        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
     2038        SetModelParms()
     2039        covMatrix = covM*Rvals['GOF']
     2040        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
     2041    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
     2042        print Msg
     2043        return False,0,0,0,0,0,0,Msg
     2044
     2045def REFDModelFxn(Profile,Inst,Limits,Substances,data):
    19162046   
    19172047    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]

  • trunk/GSASIIpwdGUI.py

    r2755 r2757  
    47064706    def OnFitModel(event):
    47074707       
    4708         print 'fit model'
    4709 #        SaveState()
    4710         G2pwd.REFDModelFxn(Profile,ProfDict,Inst,Limits,Substances,data)
    4711         G2plt.PlotPatterns(G2frame,plotType='REFD')
    4712         event.Skip()       
     4708        SaveState()
     4709        G2pwd.REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data)
     4710        wx.CallAfter(G2plt.PlotPatterns,G2frame,plotType='REFD')
     4711        wx.CallAfter(UpdateREFDModelsGrid,G2frame,data)
    47134712       
    47144713    def OnFitModelAll(event):
    … …  
    47204719        data = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,
    47214720            G2frame.PatternId,'Models'))
    4722         G2frame.dataFrame.RefdUndo.Enable(False)
     4721        G2frame.dataFrame.REFDUndo.Enable(False)
    47234722        UpdateREFDModelsGrid(G2frame,data)
    4724         G2pwd.REFDModelFxn(Profile,ProfDict,Inst,Limits,Substances,data)
     4723        G2pwd.REFDModelFxn(Profile,Inst,Limits,Substances,data)
    47254724
    47264725    def DoUnDo():
    47274726        print 'Undo last refinement'
    4728         file = open(G2frame.undosasd,'rb')
     4727        file = open(G2frame.undorefd,'rb')
    47294728        PatternId = G2frame.PatternId
    47304729        G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Models'),cPickle.load(file))
    4731         print ' Models recovered'
     4730        print ' Model recovered'
    47324731        file.close()
    47334732       
    … …  
    47364735        file = open(G2frame.undorefd,'wb')
    47374736        PatternId = G2frame.PatternId
    4738         for item in ['Models']:
    4739             cPickle.dump(G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,item)),file,1)
     4737        cPickle.dump(G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Models')),file,1)
    47404738        file.close()
    4741         G2frame.dataFrame.RefdUndo.Enable(True)
     4739        G2frame.dataFrame.REFDUndo.Enable(True)
    47424740   
    47434741    def ControlSizer():
    … …  
    47854783            if invalid:
    47864784                return
    4787             G2pwd.REFDModelFxn(Profile,ProfDict,Inst,Limits,Substances,data)
     4785            G2pwd.REFDModelFxn(Profile,Inst,Limits,Substances,data)
    47884786            G2plt.PlotPatterns(G2frame,plotType='REFD')
    47894787
    … …  
    48154813            data['Layers'][item]['Rough'] = [0.,False]
    48164814            data['Layers'][item]['Thick'] = [1.,False]
    4817             G2pwd.REFDModelFxn(Profile,ProfDict,Inst,Limits,Substances,data)
     4815            G2pwd.REFDModelFxn(Profile,Inst,Limits,Substances,data)
    48184816            G2plt.PlotPatterns(G2frame,plotType='REFD')
    48194817            wx.CallAfter(UpdateREFDModelsGrid,G2frame,data)
    … …  
    48284826            ind = Indx[Obj.GetId()]
    48294827            data['Layers'].insert(ind+1,{'Name':'vacuum','DenMul':[1.0,False],})
    4830             G2pwd.REFDModelFxn(Profile,ProfDict,Inst,Limits,Substances,data)
     4828            G2pwd.REFDModelFxn(Profile,Inst,Limits,Substances,data)
    48314829            G2plt.PlotPatterns(G2frame,plotType='REFD')
    48324830            wx.CallAfter(UpdateREFDModelsGrid,G2frame,data)
    … …  
    48364834            ind = Indx[Obj.GetId()]
    48374835            del data['Layers'][ind]
    4838             G2pwd.REFDModelFxn(Profile,ProfDict,Inst,Limits,Substances,data)
     4836            G2pwd.REFDModelFxn(Profile,Inst,Limits,Substances,data)
    48394837            G2plt.PlotPatterns(G2frame,plotType='REFD')
    48404838            wx.CallAfter(UpdateREFDModelsGrid,G2frame,data)
    … …  
    48434841            if invalid:
    48444842                return
    4845             G2pwd.REFDModelFxn(Profile,ProfDict,Inst,Limits,Substances,data)
     4843            G2pwd.REFDModelFxn(Profile,Inst,Limits,Substances,data)
    48464844            G2plt.PlotPatterns(G2frame,plotType='REFD')
    48474845            wx.CallAfter(UpdateREFDModelsGrid,G2frame,data)
    … …  
    55085506                data['delt-G(R)'][2] += ('-\n'+subData[2])
    55095507                G2plt.PlotISFG(G2frame,data,newPlot=True,plotType='delt-G(R)')
     5508                wx.CallAfter(UpdatePDFGrid,G2frame,data)
     5509               
     5510        def OnMult(invalid,value,tc):
     5511            if invalid: return
     5512            id = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,data['diffGRname'])
     5513            pId = G2gd.GetPatternTreeItemId(G2frame,id,'PDF Controls')
     5514            subData = G2frame.PatternTree.GetItemPyData(pId)['G(R)']
     5515            data['delt-G(R)'][1] = np.array([subData[1][0],data['G(R)'][1][1]-data['diffMult']*subData[1][1]])
     5516            G2plt.PlotISFG(G2frame,data,newPlot=True,plotType='delt-G(R)')
    55105517       
    55115518        diffSizer = wx.BoxSizer(wx.HORIZONTAL)
    … …  
    55155522            style=wx.CB_READONLY|wx.CB_DROPDOWN)
    55165523        grName.Bind(wx.EVT_COMBOBOX,OnSelectGR)       
    5517         diffSizer.Add(grName,0,)
     5524        diffSizer.Add(grName,0,WACV)
     5525        if data['diffGRname']:
     5526            diffSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Mult: '),0,WACV)
     5527            mult = G2G.ValidatedTxtCtrl(G2frame.dataDisplay,data,'diffMult',nDig=(10,3),
     5528                    typeHint=float,OnLeave=OnMult)
     5529            diffSizer.Add(mult,0,WACV)
     5530            OnMult(False,None,None)
    55185531        return diffSizer
    55195532           
    … …  
    57605773    if 'diffGRname' not in data:
    57615774        data['diffGRname'] = ''
     5775    if 'diffMult' not in data:
     5776        data['diffMult'] = 1.0
    57625777    if G2frame.dataDisplay:
    57635778        G2frame.dataFrame.Clear()
Note: See TracChangeset for help on using the changeset viewer.