Source code for haslib.convergence

from __future__ import division, print_function
from haslib.science_basics import A, m_He4
from multiprocessing import Pool
from tabulate import tabulate
import itertools
import numpy as np
from haslib.icc_2d_2 import iCC_2D_2 as icc2


[docs]def f(p): """ Calculate elastic intensity for a given peak Args: p: Parameters(a, h, D, xi, Ei, theta, g1, g2, gopen, gclose, sta_z, end_z, st_p_wl, cc, maxrad, m, peak) for CC simulation Returns: scalar: intensity for given peak """ # unpack parameters (a, h, D, xi, Ei, theta, g1, g2, gopen, gclose, sta_z, end_z, st_p_wl, cc, maxrad, m, peak) = p # use elastic simulation only floq = 0 AP = 0 phon = [0, 0, 0] # simulate iccResult = icc2(h, AP, a, D, xi, Ei, theta, g1, g2, gopen, gclose, phon, floq, sta_z, end_z, st_p_wl, CC_border=cc, maxrad=maxrad, verbose=False) # get result for specified peak spec_0 = iccResult[(iccResult[:, 0] == peak[0]) & (iccResult[:, 1] == peak[1])] # return peak intensity return spec_0[:, 6]
[docs]def auto_conv(D, xi, h, beam, theta, crystal, maxrad=20, cpus=1, m=m_He4, error=1e-7, gclosemax=70, st_p_wl_max=120, zrange=(20, 50), peak=(0, 0), verbose=False, rel=False): """ Find convergent parameter z, gclose and steps per wavelength Args: D: potentialdepth in meV xi: stiffness in 1/A h: corrugation beam: beam energy in meV theta: in degrees crystal: Crystal class see haslib.crystal_class maxrad: max radius(Default: 20) cpus: cpus to use in worker Pool(Default: 1) m: beam mass(Default: m_He4) error: minimum accuracy(Default: 1e-7) gclosemax: maximum closed channels to check(Default: 70) st_p_wl_max: maximum steps per wavelength to check(Default: 120) zrange: tuple of (start_z, end_z) to check in A(Default: (20, 50)) peak: tuple of two indices for g1 and g2(Default: (0, 0)) verbose: show results for each iteration(Default: False) rel: use error as relative error(Default: False) Returns: array: array(peak_intensity, gclose, end_z, start_z, st_p_wl) """ # get g vectors for given crystal g1 = crystal.b1[0:2] g2 = crystal.b2[0:2] # use all open channels gopen = -1 # define parameter range to check gclose = np.arange(0, gclosemax, 2).astype(int) end_z = np.arange(0, zrange[1], 1)*A sta_z = -np.arange(0, zrange[0], 1)*A st_p_wl = np.arange(10, st_p_wl_max, 5) # save all parameterranges as an array params = [gclose, end_z, sta_z, st_p_wl] # set the last values of the defines ranges as default/starting parameters default = [np.array([param[-1]]) for param in params] # loop over every type of parameter for ind in range(len(params)): # take the default parameters and replace the current # parameter type with its range p = default[:] p[ind] = params[ind] # create search space out of 3 default parameters # and one parameter range search = list(itertools.product(*p)) # create parameter list for each item in the search space # and store in payload # these parameters will be given to the cc-simulation cc = 10e-5 payload = [] for p in search: gclose, end_z, sta_z, st_p_wl = p payload.append([crystal.a1[0], h, D, xi, beam, theta, g1, g2, int(gopen), int(gclose), sta_z, end_z, int(st_p_wl), cc, maxrad, m, peak]) # check if multiple threads should be used and and run simulations if cpus >= 1: p = Pool(cpus) res = p.map(f, payload) p.close() else: res = map(f, payload) # create array from search space list search = np.vstack(search) # plot current results if verbose is defined if verbose: print(tabulate(np.hstack([res, search]), headers=["gclose", "end", "start", "step"])) # find the results which are converged if rel: # use a relative error good = abs(res/res[-1]-1) <= error else: # use an absolute error good = abs(res-res[-1]) <= error # search for position where the result converged ind2 = len(res)-1 for ind3 in range(len(res)): # check if all results are good after index ind3 if good[ind3:].all(): ind2 = ind3 break # if there is no convergence raise an exception if ind2 == len(res)-1: d = {0: "gclose", 1: "end_z", 2: "start_z", 3: "steps"} raise Warning("Boundries too low - {}({})".format(d[ind], ind)) # set the converged parameter as the new default parameter default[ind] = np.array([search[ind2, ind]]) # print default parameters if verbose=True if verbose: print(tabulate(np.array(default)[None, :], headers=["gclose", "end", "start", "step"])) # return default parameters after every parameter has converged return np.array(default)