Changeset 2937 for branch/2frame/GSASIImath.py
 Timestamp:
 Jul 19, 2017 4:14:08 PM (5 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branch/2frame/GSASIImath.py
r2930 r2937 3883 3883 3884 3884 import numpy 3885 from numpy import asarray, tan,exp, squeeze, sign, \3886 all, log, pi,shape, array, where3885 from numpy import asarray, exp, squeeze, sign, \ 3886 all, shape, array, where 3887 3887 from numpy import random 3888 3888 … … 3894 3894 def __init__(self): 3895 3895 self.dwell = 20 3896 self.learn_rate = 0.53897 3896 self.lower = 10 3898 3897 self.upper = 10 … … 3960 3959 self.accepted += 1 3961 3960 return 1 3962 p = exp(dE*1.0/ self.boltzmann/T)3961 p = exp(dE*1.0/T) 3963 3962 if (p > random.uniform(0.0, 1.0)): 3964 3963 self.accepted += 1 … … 3972 3971 pass 3973 3972 3974 3975 3973 # A schedule due to Lester Ingber modified to use bounds  OK 3976 3974 class fast_sa(base_schedule): 3977 3975 def init(self, **options): 3978 3976 self.__dict__.update(options) 3979 if self.m is None:3980 self.m = 1.03981 if self.n is None:3982 self.n = 1.03983 self.c = self.m * exp(self.n * self.quench)3984 3977 3985 3978 def update_guess(self, x0): … … 3990 3983 xnew = xc*(self.upper  self.lower)+self.lower 3991 3984 return xnew 3992 # y = sign(u0.5)*T*((1+1.0/T)**abs(2*u1)1.0)3993 # xc = y*(self.upper  self.lower)3994 # xnew = x0 + xc3995 # return xnew3996 3985 3997 3986 def update_temp(self): … … 4000 3989 return 4001 3990 4002 class cauchy_sa(base_schedule): #modified to use bounds  not good4003 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.lower4008 return xnew4009 # numbers = squeeze(random.uniform(pi/2, pi/2, size=self.dims))4010 # xc = self.learn_rate * self.T * tan(numbers)4011 # xnew = x0 + xc4012 # return xnew4013 4014 def update_temp(self):4015 self.T = self.T0/(1+self.k)4016 self.k += 14017 return4018 4019 class boltzmann_sa(base_schedule):4020 # def update_guess(self, x0):4021 # std = minimum(sqrt(self.T)*ones(self.dims), (self.upperself.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_rate4026 # return xnew4027 4028 def update_temp(self):4029 self.k += 14030 self.T = self.T0 / log(self.k+1.0)4031 return4032 4033 3991 class log_sa(base_schedule): #OK 4034 3992 … … 4037 3995 4038 3996 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.upperself.lower)+self.lower 3997 u = squeeze(random.uniform(0.0, 1.0, size=self.dims)) 3998 T = self.T 3999 xc = (sign(u0.5)*T*((1+1.0/T)**abs(2*u1)1.0)+1.0)/2.0 4000 xnew = xc*(self.upper  self.lower)+self.lower 4001 return xnew 4040 4002 4041 4003 def update_temp(self): … … 4048 4010 self.cost = None 4049 4011 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 4012 def 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 4057 4028 def anneal(func, x0, args=(), schedule='fast', full_output=0, 4058 4029 T0=None, Tf=1e12, maxeval=None, maxaccept=None, maxiter=400, 4059 boltzmann=1.0, learn_rate=0.5, feps=1e6, quench=1.0, m=1.0, n=1.0,4030 feps=1e6, quench=1.0, c=1.0, 4060 4031 lower=100, upper=100, dwell=50, slope=0.9,ranStart=False, 4061 4032 ranRange=0.10,autoRan=False,dlg=None): … … 4088 4059 :param float learn_rate: 4089 4060 Scale constant for adjusting guesses. 4090 :param float boltzmann:4091 Boltzmann constant in acceptance test4092 (increase for less stringent test at each temperature).4093 4061 :param float feps: 4094 4062 Stopping relative error tolerance for the function value in 4095 4063 last four coolings. 4096 :param float quench, m,n:4064 :param float quench,c: 4097 4065 Parameters to alter fast_sa schedule. 4098 4066 :param float/ndarray lower,upper: … … 4156 4124 x_new = x_old + xc 4157 4125 4158 c = n * exp(n * quench)4159 4126 T_new = T0 * exp(c * k**quench) 4160 4127 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 + xc4167 4168 T_new = T0 / (1+k)4169 4170 In the 'boltzmann' schedule the updates are ::4171 4172 std = minimum( sqrt(T) * ones(d), (upperlower) / (3*learn_rate) )4173 y ~ Normal(0, std, size=d)4174 x_new = x_old + learn_rate * y4175 4176 T_new = T0 / log(1+k)4177 4128 4178 4129 """ … … 4183 4134 schedule = eval(schedule+'_sa()') 4184 4135 # 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) 4188 4138 4189 4139 current_state, last_state, best_state = _state(), _state(), _state() … … 4258 4208 if abs(af[1]best_state.cost) > feps*10: 4259 4209 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) 4263 4211 break 4264 4212 if (Tf is not None) and (schedule.T < Tf): 4213 print ' Minimum T reached' 4265 4214 retval = 1 4266 4215 break … … 4269 4218 break 4270 4219 if (iters > maxiter): 4271 print "Warning: Maximum number of iterations exceeded."4220 print " Warning: Maximum number of iterations exceeded." 4272 4221 retval = 3 4273 4222 break … … 4552 4501 M = np.inner(refList[6],np.inner(rcov,refList[6])) 4553 4502 tsum += (time.time()t0) 4554 return M/np.sum(refList[4]**2)4503 return np.sqrt(M/np.sum(refList[4]**2)) 4555 4504 4556 4505 def MCSAcallback(x, f, fmin,accept): … … 4558 4507 newmsg='%s%8.4f%s'%('MC/SA Residual:',fmin*100,'%'))[0] 4559 4508 4560 sq8ln2 = np.sqrt(8*np.log(2))4561 4509 sq2pi = np.sqrt(2*np.pi) 4562 4510 sq4pi = np.sqrt(4*np.pi) … … 4620 4568 if d >= MCSA['dmin']: 4621 4569 sig = np.sqrt(sig) #var > sig in centideg 4622 sig = G2pwd.getgamFW(gam,sig)/sq8ln2 #FWHM > sig in centideg4570 sig = .01*G2pwd.getgamFW(gam,sig) #sig,gam > FWHM in deg 4623 4571 SQ = 0.25/d**2 4624 4572 allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM)) … … 4690 4638 x0 = [parmDict[val] for val in varyList] 4691 4639 ifInv = SGData['SGInv'] 4640 bounds = np.array(zip(lower,upper)) 4692 4641 if MCSA['Algorithm'] == 'Basin Hopping': 4693 4642 import basinhopping as bs 4694 bounds = np.array(zip(lower,upper))4695 4643 take_step = RandomDisplacementBounds(np.array(lower), np.array(upper)) 4696 4644 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':'LBFGSB','bounds':bounds,4645 interval=MCSA['Annealing'][2]/10,niter=MCSA['Annealing'][2],minimizer_kwargs={'method':'LBFGSB','bounds':bounds, 4698 4646 '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'])4701 4647 else: 4702 4648 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], 4707 4652 lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False), 4708 4653 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='LBFGSB',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]) 4714 4659 Result.append(varyList) 4715 return Result,tsum 4660 return Result,tsum,rcov 4716 4661 4717 4662
Note: See TracChangeset
for help on using the changeset viewer.