Changeset 2774 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Apr 7, 2017 4:05:51 PM (5 years ago)
Author:
vondreele
Message:

implement MC/SA for reflectometry
more image importer speed ups

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r2772 r2774  
    19171917    print 'fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0])
    19181918   
     1919    class RandomDisplacementBounds(object):
     1920        """random displacement with bounds"""
     1921        def __init__(self, xmin, xmax, stepsize=0.5):
     1922            self.xmin = xmin
     1923            self.xmax = xmax
     1924            self.stepsize = stepsize
     1925   
     1926        def __call__(self, x):
     1927            """take a random step but ensure the new position is within the bounds"""
     1928            while True:
     1929                # this could be done in a much more clever way, but it will work for example purposes
     1930                steps = self.xmax-self.xmin
     1931                xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x))
     1932                if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
     1933                    break
     1934            return xnew
     1935   
    19191936    def GetModelParms():
    19201937        parmDict = {}
     
    19221939        values = []
    19231940        bounds = []
    1924         parmDict['Res'] = data['Resolution'][0]/100.       #%-->decimal
     1941        parmDict['Res'] = data['Resolution'][0]/(100.*np.sqrt(ateln2))     #% FWHM-->decimal sig
    19251942        for parm in ['Scale','FltBack']:
    19261943            parmDict[parm] = data[parm][0]
     
    19381955                if layer.get(parm,[0.,False])[1]:
    19391956                    varyList.append(cid+parm)
    1940                     values.append(layer[parm][0])
    1941                     bounds.append(Bounds[parm])
     1957                    value = layer[parm][0]
     1958                    values.append(value)
     1959                    if value:
     1960                        bound = [value*Bfac,value/Bfac]
     1961                    else:
     1962                        bound = [0.,10.]
     1963                    bounds.append(bound)
    19421964            if name not in ['vacuum','unit scatter']:
    19431965                parmDict[cid+'rho'] = Substances[name]['Scatt density']
     
    20072029        Ic += (A**2+B**2)*Scale     
    20082030        return Ic
     2031       
     2032    def estimateT0(takestep):
     2033        Mmax = -1.e-10
     2034        Mmin = 1.e10
     2035        for i in range(100):
     2036            x0 = takestep(values)
     2037            M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)
     2038            Mmin = min(M,Mmin)
     2039            MMax = max(M,Mmax)
     2040        return 1.5*(MMax-Mmin)
    20092041
    20102042    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
     
    20142046    Qmax = Limits[1][1]
    20152047    wtFactor = ProfDict['wtFactor']
     2048    Bfac = data['Toler']
    20162049    Ibeg = np.searchsorted(Q,Qmin)
    20172050    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
    20182051    Ic[:] = 0
    2019     Bounds = {'Scale':[data['Scale'][0]*.85,data['Scale'][0]/.85],'FltBack':[None,None],
    2020               'DenMul':[None,None],'Thick':[1.,None],'Rough':[0.,None],'Mag SLD':[-10.,10.],'iDenMul':[None,None]}
     2052    Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6],
     2053              'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]}
    20212054    parmDict,varyList,values,bounds = GetModelParms()
    20222055    Msg = 'Failed to converge'
     
    20302063            covM = result[1]
    20312064            newVals = result[0]
    2032         elif data['Minimizer'] == 'Global':
    2033             result = so.basinhopping(sumREFD,values,minimizer_kwargs={'method':'L-BFGS-B',
     2065        elif data['Minimizer'] == 'Basin Hopping':
     2066            xyrng = np.array(bounds).T
     2067            take_step = RandomDisplacementBounds(xyrng[0], xyrng[1])
     2068            T0 = estimateT0(take_step)
     2069            print ' Estimated temperature: %.3g'%(T0)
     2070            result = so.basinhopping(sumREFD,values,take_step=take_step,disp=True,T=T0,stepsize=Bfac,
     2071                interval=20,niter=200,minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
    20342072                'args':(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList)})
    20352073            chisq = result.fun
     
    20372075            newVals = result.x
    20382076            covM = []
     2077        elif data['Minimizer'] == 'MC/SA Anneal':
     2078            xyrng = np.array(bounds).T
     2079            result = G2mth.anneal(sumREFD, values, args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],parmDict,varyList),
     2080                schedule='log', full_output=True,maxeval=None, maxaccept=None, maxiter=10,dwell=1000,
     2081                boltzmann=10.0, feps=1e-6,lower=xyrng[0], upper=xyrng[1], slope=0.9,ranStart=True,
     2082                ranRange=0.20,autoRan=False,dlg=None)
     2083            newVals = result[0]
     2084            parmDict.update(zip(varyList,newVals))
     2085            chisq = result[1]
     2086            ncalc = result[3]
     2087            covM = []
     2088            print ' MC/SA final temperature: %.4g'%(result[2])
    20392089        elif data['Minimizer'] == 'L-BFGS-B':
    20402090            result = so.minimize(sumREFD,values,method='L-BFGS-B',bounds=bounds,   #ftol=Ftol,
Note: See TracChangeset for help on using the changeset viewer.