Changeset 2758 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Mar 21, 2017 3:34:22 PM (5 years ago)
Author:
vondreele
Message:

implement L-BFGS-B and basinhopping (preliminary) minimizers for reflectometry

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r2757 r2758  
    19141914
    19151915def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data):
    1916     print 'fit REFD data'
     1916    print 'fit REFD data by '+data['Minimizer']
    19171917   
    19181918    def GetModelParms():
     
    19201920        varyList = []
    19211921        values = []
     1922        bounds = []
    19221923        for parm in ['Scale','FltBack']:
    19231924            parmDict[parm] = data[parm][0]
     
    19251926                varyList.append(parm)
    19261927                values.append(data[parm][0])
     1928                bounds.append(Bounds[parm])
    19271929        parmDict['nLayers'] = len(data['Layers'])
    19281930        for ilay,layer in enumerate(data['Layers']):
     
    19341936                    varyList.append(cid+parm)
    19351937                    values.append(layer[parm][0])
     1938                    bounds.append(Bounds[parm])
    19361939            parmDict[cid+'rho'] = Substances[name]['Scatt density']
    19371940            parmDict[cid+'irho'] = Substances[name].get('XImag density',0.)
    1938         return parmDict,varyList,values
     1941        return parmDict,varyList,values,bounds
    19391942   
    19401943    def SetModelParms():
     
    19701973        return M
    19711974   
     1975    def sumREFD(values,Q,Io,wt,parmDict,varyList):
     1976        parmDict.update(zip(varyList,values))
     1977        M = np.sqrt(wt)*(getREFD(Q,parmDict)-Io)
     1978        return np.sum(M**2)
     1979   
    19721980    def getREFD(Q,parmDict):
    19731981        Ic = np.ones_like(Q)*parmDict['FltBack']
     
    19952003    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
    19962004    Ic[:] = 0
    1997     parmDict,varyList,values = GetModelParms()
     2005    Bounds = {'Scale':[data['Scale'][0]*.85,data['Scale'][0]/.85],'FltBack':[None,None],'DenMul':[0.,None],'Thick':[1.,None],'Rough':[0.,None]}
     2006    parmDict,varyList,values,bounds = GetModelParms()
     2007    Msg = 'Failed to converge'
    19982008    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]
     2009        if data['Minimizer'] == 'LMLS':
     2010            result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
     2011                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList))
     2012            parmDict.update(zip(varyList,result[0]))
     2013            chisq = np.sum(result[2]['fvec']**2)
     2014            ncalc = result[2]['nfev']
     2015            covM = result[1]
     2016            newVals = result[0]
     2017        elif data['Minimizer'] == 'Global':
     2018            result = so.basinhopping(sumREFD,values,minimizer_kwargs={'method':'L-BFGS-B',
     2019                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)})
     2020            chisq = result.fun
     2021            ncalc = result.nfev
     2022            newVals = result.x
     2023            covM = []
     2024        elif data['Minimizer'] == 'L-BFGS-B':
     2025            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
     2026                args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList))
     2027            parmDict.update(zip(varyList,result['x']))
     2028            chisq = result.fun
     2029            ncalc = result.nfev
     2030            newVals = result.x
     2031            covM = []
    20052032    else:   #nothing varied
    20062033        M = calcREFD(values,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)
     
    20162043    Ic[Ibeg:Ifin] = getREFD(Q[Ibeg:Ifin],parmDict)
    20172044    Ib[Ibeg:Ifin] = parmDict['FltBack']
    2018     Msg = 'Failed to converge'
    20192045    try:
    2020         Nans = np.isnan(result[0])
     2046        if not len(varyList):
     2047            Msg += ' - nothing refined'
     2048            raise ValueError
     2049        Nans = np.isnan(newVals)
    20212050        if np.any(Nans):
    20222051            name = varyList[Nans.nonzero(True)[0]]
    20232052            Msg += ' Nan result for '+name+'!'
    20242053            raise ValueError
    2025         Negs = np.less_equal(result[0],0.)
     2054        Negs = np.less_equal(newVals,0.)
    20262055        if np.any(Negs):
    20272056            indx = Negs.nonzero()
     
    20322061        if len(covM):
    20332062            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
    2034             sigDict = dict(zip(varyList,sig))
     2063            covMatrix = covM*Rvals['GOF']
     2064        else:
     2065            sig = np.zeros(len(varyList))
     2066            covMatrix = []
     2067        sigDict = dict(zip(varyList,sig))
    20352068        print ' Results of reflectometry data modelling fit:'
    20362069        print 'Number of function calls:',ncalc,' Number of observations: ',Ifin-Ibeg,' Number of parameters: ',len(varyList)
    20372070        print 'Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])
    20382071        SetModelParms()
    2039         covMatrix = covM*Rvals['GOF']
    20402072        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
    20412073    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
Note: See TracChangeset for help on using the changeset viewer.