Ignore:
Timestamp:
Jul 19, 2017 4:14:08 PM (4 years ago)
Author:
vondreele
Message:

put a self.SetAutoLayout?(True) in SetDataSize?
cleanup anneal - remove boltzmann & cauchy option & parameters - fundamentlly flawed
add new function (makeTsched) to generate anneal T schedule - show plot if coeff changed
add a final 'L-BFGS-B' refinement of each anneal result - gives much improved results
show plot of reflection covariances for MC/SA runs
cleanup MC/SA GUI stuff
cleanup Index list stuff - fix bug as well

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branch/2frame/GSASIImath.py

    r2930 r2937  
    38833883
    38843884import numpy
    3885 from numpy import asarray, tan, exp, squeeze, sign, \
    3886      all, log, pi, shape, array, where
     3885from numpy import asarray, exp, squeeze, sign, \
     3886     all, shape, array, where
    38873887from numpy import random
    38883888
     
    38943894    def __init__(self):
    38953895        self.dwell = 20
    3896         self.learn_rate = 0.5
    38973896        self.lower = -10
    38983897        self.upper = 10
     
    39603959            self.accepted += 1
    39613960            return 1
    3962         p = exp(-dE*1.0/self.boltzmann/T)
     3961        p = exp(-dE*1.0/T)
    39633962        if (p > random.uniform(0.0, 1.0)):
    39643963            self.accepted += 1
     
    39723971        pass
    39733972
    3974 
    39753973#  A schedule due to Lester Ingber modified to use bounds - OK
    39763974class fast_sa(base_schedule):
    39773975    def init(self, **options):
    39783976        self.__dict__.update(options)
    3979         if self.m is None:
    3980             self.m = 1.0
    3981         if self.n is None:
    3982             self.n = 1.0
    3983         self.c = self.m * exp(-self.n * self.quench)
    39843977
    39853978    def update_guess(self, x0):
     
    39903983        xnew = xc*(self.upper - self.lower)+self.lower
    39913984        return xnew
    3992 #        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
    3993 #        xc = y*(self.upper - self.lower)
    3994 #        xnew = x0 + xc
    3995 #        return xnew
    39963985
    39973986    def update_temp(self):
     
    40003989        return
    40013990
    4002 class cauchy_sa(base_schedule):     #modified to use bounds - not good
    4003     def update_guess(self, x0):
    4004         x0 = asarray(x0)
    4005         numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
    4006         xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
    4007         xnew = xc*(self.upper - self.lower)+self.lower
    4008         return xnew
    4009 #        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
    4010 #        xc = self.learn_rate * self.T * tan(numbers)
    4011 #        xnew = x0 + xc
    4012 #        return xnew
    4013 
    4014     def update_temp(self):
    4015         self.T = self.T0/(1+self.k)
    4016         self.k += 1
    4017         return
    4018 
    4019 class boltzmann_sa(base_schedule):
    4020 #    def update_guess(self, x0):
    4021 #        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
    4022 #        x0 = asarray(x0)
    4023 #        xc = squeeze(random.normal(0, 1.0, size=self.dims))
    4024 #
    4025 #        xnew = x0 + xc*std*self.learn_rate
    4026 #        return xnew
    4027 
    4028     def update_temp(self):
    4029         self.k += 1
    4030         self.T = self.T0 / log(self.k+1.0)
    4031         return
    4032 
    40333991class log_sa(base_schedule):        #OK
    40343992
     
    40373995       
    40383996    def update_guess(self,x0):     #same as default #TODO - is this a reasonable update procedure?
    4039         return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
     3997        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
     3998        T = self.T
     3999        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
     4000        xnew = xc*(self.upper - self.lower)+self.lower
     4001        return xnew
    40404002       
    40414003    def update_temp(self):
     
    40484010        self.cost = None
    40494011
    4050 # TODO:
    4051 #     allow for general annealing temperature profile
    4052 #     in that case use update given by alpha and omega and
    4053 #     variation of all previous updates and temperature?
    4054 
    4055 # Simulated annealing   #TODO - should we switch to scipy basinhopping?
    4056 
     4012def makeTsched(data):
     4013    if data['Algorithm'] == 'fast':
     4014        sched = fast_sa()
     4015        sched.quench = data['fast parms'][0]
     4016        sched.c = data['fast parms'][1]
     4017    elif data['Algorithm'] == 'log':
     4018        sched = log_sa()
     4019        sched.slope = data['log slope']
     4020    sched.T0 = data['Annealing'][0]
     4021    Tf = data['Annealing'][1]
     4022    Tsched = [sched.T0,]
     4023    while Tsched[-1] > Tf:
     4024        sched.update_temp()
     4025        Tsched.append(sched.T)
     4026    return Tsched[1:]
     4027   
    40574028def anneal(func, x0, args=(), schedule='fast', full_output=0,
    40584029           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
    4059            boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
     4030           feps=1e-6, quench=1.0, c=1.0,
    40604031           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
    40614032           ranRange=0.10,autoRan=False,dlg=None):
     
    40884059    :param float learn_rate:
    40894060        Scale constant for adjusting guesses.
    4090     :param float boltzmann:
    4091         Boltzmann constant in acceptance test
    4092         (increase for less stringent test at each temperature).
    40934061    :param float feps:
    40944062        Stopping relative error tolerance for the function value in
    40954063        last four coolings.
    4096     :param float quench,m,n:
     4064    :param float quench,c:
    40974065        Parameters to alter fast_sa schedule.
    40984066    :param float/ndarray lower,upper:
     
    41564124        x_new = x_old + xc
    41574125
    4158         c = n * exp(-n * quench)
    41594126        T_new = T0 * exp(-c * k**quench)
    41604127
    4161 
    4162     In the 'cauchy' schedule the updates are ::
    4163 
    4164         u ~ Uniform(-pi/2, pi/2, size=d)
    4165         xc = learn_rate * T * tan(u)
    4166         x_new = x_old + xc
    4167 
    4168         T_new = T0 / (1+k)
    4169 
    4170     In the 'boltzmann' schedule the updates are ::
    4171 
    4172         std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
    4173         y ~ Normal(0, std, size=d)
    4174         x_new = x_old + learn_rate * y
    4175 
    4176         T_new = T0 / log(1+k)
    41774128
    41784129    """
     
    41834134    schedule = eval(schedule+'_sa()')
    41844135    #   initialize the schedule
    4185     schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
    4186                   learn_rate=learn_rate, lower=lower, upper=upper,
    4187                   m=m, n=n, quench=quench, dwell=dwell, slope=slope)
     4136    schedule.init(dims=shape(x0),func=func,args=args,T0=T0,lower=lower, upper=upper,
     4137        c=c, quench=quench, dwell=dwell, slope=slope)
    41884138
    41894139    current_state, last_state, best_state = _state(), _state(), _state()
     
    42584208            if abs(af[-1]-best_state.cost) > feps*10:
    42594209                retval = 5
    4260 #                print "Warning: Cooled to %f at %s but this is not" \
    4261 #                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
    4262 #                      + " the smallest point found."
     4210                print " Warning: Cooled to %.4f > selected Tmin %.4f"%(squeeze(last_state.cost),Tf)
    42634211            break
    42644212        if (Tf is not None) and (schedule.T < Tf):
     4213            print ' Minimum T reached'
    42654214            retval = 1
    42664215            break
     
    42694218            break
    42704219        if (iters > maxiter):
    4271             print "Warning: Maximum number of iterations exceeded."
     4220            print  " Warning: Maximum number of iterations exceeded."
    42724221            retval = 3
    42734222            break
     
    45524501        M = np.inner(refList[6],np.inner(rcov,refList[6]))
    45534502        tsum += (time.time()-t0)
    4554         return M/np.sum(refList[4]**2)
     4503        return np.sqrt(M/np.sum(refList[4]**2))
    45554504   
    45564505    def MCSAcallback(x, f, fmin,accept):
     
    45584507            newmsg='%s%8.4f%s'%('MC/SA Residual:',fmin*100,'%'))[0]
    45594508
    4560     sq8ln2 = np.sqrt(8*np.log(2))
    45614509    sq2pi = np.sqrt(2*np.pi)
    45624510    sq4pi = np.sqrt(4*np.pi)
     
    46204568            if d >= MCSA['dmin']:
    46214569                sig = np.sqrt(sig)      #var -> sig in centideg
    4622                 sig = G2pwd.getgamFW(gam,sig)/sq8ln2        #FWHM -> sig in centideg
     4570                sig = .01*G2pwd.getgamFW(gam,sig)        #sig,gam -> FWHM in deg
    46234571                SQ = 0.25/d**2
    46244572                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
     
    46904638    x0 = [parmDict[val] for val in varyList]
    46914639    ifInv = SGData['SGInv']
     4640    bounds = np.array(zip(lower,upper))
    46924641    if MCSA['Algorithm'] == 'Basin Hopping':
    46934642        import basinhopping as bs
    4694         bounds = np.array(zip(lower,upper))
    46954643        take_step = RandomDisplacementBounds(np.array(lower), np.array(upper))
    46964644        results = bs.basinhopping(mcsaCalc,x0,take_step=take_step,disp=True,T=MCSA['Annealing'][0],
    4697                 interval=MCSA['Annealing'][2]/5,niter=MCSA['Annealing'][2],minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
     4645                interval=MCSA['Annealing'][2]/10,niter=MCSA['Annealing'][2],minimizer_kwargs={'method':'L-BFGS-B','bounds':bounds,
    46984646                'args':(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)},callback=MCSAcallback)
    4699         mcsaCalc(results['x'],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
    4700         Result = [False,False,results['fun'],0.0,]+list(results['x'])
    47014647    else:
    47024648        results = anneal(mcsaCalc,x0,args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict),
    4703             schedule=MCSA['Algorithm'], full_output=True,
    4704             T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
    4705             boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
    4706             quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
     4649            schedule=MCSA['Algorithm'], full_output=True,dwell=MCSA['Annealing'][2],maxiter=10000,
     4650            T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],
     4651            quench=MCSA['fast parms'][0], c=MCSA['fast parms'][1],
    47074652            lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False),
    47084653            ranRange=MCSA.get('ranRange',10.)/100.,autoRan=MCSA.get('autoRan',False),dlg=pgbar)
    4709         mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
    4710         Result = [False,False,results[1],results[2],]+list(results[0])
    4711 #    for ref in refs.T:
    4712 #        print ' %4d %4d %4d %10.3f %10.3f %10.3f'%(int(ref[0]),int(ref[1]),int(ref[2]),ref[4],ref[5],ref[6])
    4713 #    print np.sqrt((np.sum(refs[6]**2)/np.sum(refs[4]**2)))
     4654        results = so.minimize(mcsaCalc,results[0],method='L-BFGS-B',args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict),
     4655            bounds=bounds,)
     4656    mcsaCalc(results['x'],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
     4657    Result = [False,False,results['fun'],0.0,]+list(results['x'])
     4658#        Result = [False,False,results[1],results[2],]+list(results[0])
    47144659    Result.append(varyList)
    4715     return Result,tsum
     4660    return Result,tsum,rcov
    47164661
    47174662       
Note: See TracChangeset for help on using the changeset viewer.