Changeset 2757 for trunk/GSASIIpwd.py


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

workable reflectometry fitting by Levenberg-Marquart least squares.

File:
1 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]
Note: See TracChangeset for help on using the changeset viewer.