# -*- coding: utf-8 -*- """ Created on Thu Aug 16 00:07:08 2018 @author: USER """ # Codes are free to use. Do whatever you want from __future__ import absolute_import """CTR data fitting""" ####################### LIBRARY ############################# # exceptions library # Python stdlib imports import pickle import os # data processing library import numpy as np import pandas as pd import scipy.optimize as so import matplotlib.pyplot as plt import openpyxl as ox # pyrod library from read.read_parameters import initialization_parameters from read.read_raw_data import initialization_rhkl import tool.control as tc import tool.tools as tt from tool.properties import (roughness, vacancy, relax) import structure ####################### CONSTANT ############################ # constant data_path = os.path.abspath(os.path.dirname('ctr_optimised_result.pickle')) +\ '/data/ctr_optimised_result.pickle' ####################### FUNCTIONS ########################### ######################## CLASSS ############################# class var_refine(object): # class for modulate the parameters maunally or opimisly def __init__(self, parameters, experiment_data): # init var_list and table self.parameters = parameters self.para_path = os.path.abspath(os.path.dirname(parameters)) + parameters self.experiment_data = experiment_data self.data_path = os.path.abspath(os.path.dirname(experiment_data)) + experiment_data init_para = initialization_parameters(self.parameters) init_para._check_parameters() init_para._check_sheetname() init_data = initialization_rhkl(self.experiment_data) init_data._check_filename() # real time read xlsx, real time refresh self.var_list = init_para._var_list() self.var_table = init_para._var_table() # import factors self.properties = pd.read_excel(self.para_path,sheet_name = 'factors',index_col = 0) # init experiment data iq and intensity self.iq,self.i = init_data._read_rhkl(experiment_data) # whick rod self.h = self.properties.at['h','data'] self.k = self.properties.at['k','data'] # inti parameters self.absorption = self.properties.at['absorption','data'] self.roughness = self.properties.at['roughness','data'] self.scale = self.properties.at['scale','data'] self.intensity = self.properties.at['intensity','data'] # interface roughness self.iroughness = self.properties.at['interface roughness','data'] # contral parameters self.ubr_index = tc.contral_vars(self.var_list) def track_difference(self): var_list_rei = initialization_parameters._var_list(self.parameters) track_flag = 0 try: for slabi in var_list_rei: for keyi in var_list_rei[slabi]: # compare every key value in new excel and old excel if keyi in ['absorption','atoms_num','layers_num','roughness']: # refreshed data and old data value_rei = var_list_rei[slabi][keyi] value_old = self.var_list[slabi][keyi] # if same,continue. if change,flag = 1 if np.equal(value_rei,value_old): continue else: track_flag = 1 # these parameters are matrix elif keyi in ['dis','dw','ions','lattice','nel','ocu','pos']: # refreshed matrix and old matrix if False in np.equal(var_list_rei[slabi][keyi], self.var_list[slabi][keyi]).as_matrix(): track_flag = 1 break else: continue # beak from the inner loop if track_flag: break # break from the outer loop if track_flag: break except: track_flag = 1 properties = pd.read_excel(self.para_path,sheet_name = 'factors',index_col = 0) if False in np.equal(properties,self.properties).as_matrix(): track_flag = 1 # refresh if track_flag == 1: init_para = initialization_parameters(self.parameters) init_data = initialization_rhkl(self.experiment_data) self.var_list = init_para._var_list() self.var_table = init_para._var_table() # import factors self.properties = pd.read_excel(self.para_path,sheet_name = 'factors',index_col = 0) # init experiment data iq and intensity self.iq,self.i = init_data._read_rhkl(self.experiment_data) # whick rod self.h = self.properties.at['h','data'] self.k = self.properties.at['k','data'] # inti parameters self.absorption = self.properties.at['absorption','data'] self.roughness = self.properties.at['roughness','data'] self.scale = self.properties.at['scale','data'] self.intensity = self.properties.at['intensity','data'] # interface roughness self.iroughness = self.properties.at['interface roughness','data'] # contral parameters self.ubr_index = tc.contral_vars(self.var_list) return self,track_flag def re(self): init_para = initialization_parameters(self.parameters) init_data = initialization_rhkl(self.experiment_data) self.var_list = init_para._var_list() self.var_table = init_para._var_table() init_para = initialization_parameters(self.parameters) init_data = initialization_rhkl(self.experiment_data) # import factors self.properties = pd.read_excel(self.para_path,sheet_name = 'factors',index_col = 0) # init experiment data iq and intensity self.iq,self.i = init_data._read_rhkl(self.experiment_data) # whick rod self.h = self.properties.at['h','data'] self.k = self.properties.at['k','data'] # inti parameters self.absorption = self.properties.at['absorption','data'] self.roughness = self.properties.at['roughness','data'] self.scale = self.properties.at['scale','data'] self.intensity = self.properties.at['intensity','data'] # interface roughness self.iroughness = self.properties.at['interface roughness','data'] def disp(self): # list all the self parameters print("--------------------------------------------") print("parameters: %s\n" % self.parameters) print("experiment_data: %s\n" % self.experiment_data) print("h k qz: %s %s %s\n" % (int(self.h), int(self.k), [np.min(self.iq),np.max(self.iq)])) print("intensity: %s\n" % self.intensity) print("absorption: %s\n" % self.absorption) print("roughness: %s\n" % self.roughness) print("scale: %s\n" % self.scale) print("var_list: %s\n" %self.var_list.keys()) print('var_table: %s' %self.var_table.keys()) print("roughness interface: %s" %self.iroughness) print('--------------------------------------------') def sub_ctr(self): # the used keys key_list = ['absorption', 'roughness', 'scale'] # input substrate ubr and substrate index subr,sindex = tc.initialize_contral_vars(self.ubr_index, ['substrate'], key_list) substrate_ctr = structure.substrate_ctr(self.var_list, self.iq, subr, sindex, self.absorption, self.h, self.h, self.scale, key_list) # roughness cr = roughness(self.iq) r = cr.robinson_roughness(self.roughness) return np.multiply(self.intensity*substrate_ctr, np.mat(r).T) def slab_ctr(self): # used keys key_list = ['lattice_abc','dw'] # input slab ubr and slab index slab_ubr,slab_index = tc.initialize_contral_vars(self.ubr_index, np.unique(self.var_table['slab_list']).tolist(), key_list) fctr = structure.film_ctr(self.var_list, self.var_table, slab_ubr, slab_index, self.iq, self.h, self.k, key_list) # roughness cr = roughness(self.iq) r = cr.robinson_roughness(self.roughness) # interface roughness ir = cr.interface_roughness(self.iroughness) return np.multiply(np.multiply(self.intensity*fctr, np.mat(r).T),np.mat(ir).T) def lattice_refine(self, p_mask = 13, n_mask = 3, weight = 1): # Prepare the data mask # yin mask to mask the bragg peaks bn_mask = np.mat(tc.bragg_mask(self.iq, n_mask,1, mode = 'yin')) # yang mask to refine the signal around bragg peaks only bp_mask = np.mat(tc.bragg_mask(self.iq, p_mask,1, mode = 'yang')) bgm = np.asarray(np.multiply(bn_mask, bp_mask).T).reshape(-1)*weight # calculate ss tot.For only seveal signal points are calculate # mean value is np.sum(abs(i)*bgm)/np.sum(bgm) ss_tot = np.sum((abs(self.i)*bgm - \ np.sum(abs(self.i)*bgm)/np.sum(bgm))**2) # initialize the variable about the lattice slab_list = np.unique(self.var_table['slab_list']).tolist() subr,sindex = tc.initialize_contral_vars(self.ubr_index, slab_list, ['lattice_abc']) # substrate ctr # the used keys key_list = ['absorption', 'roughness', 'scale'] # input substrate ubr and substrate index sub_ubr,sub_index = tc.initialize_contral_vars(self.ubr_index, ['substrate'], key_list) substrate_ctr = structure.substrate_ctr(self.var_list, self.iq, sub_ubr, sub_index, self.absorption, self.h, self.h, self.scale, key_list) # roughness cr = roughness(self.iq) r = cr.robinson_roughness(self.roughness) def c_refine(c_var): # input lattice c variable into subr if len(slab_list) == 1: subr[2] = c_var elif len(slab_list) >= 2: for slabi in range(len(slab_list)): subr[slabi*3-1] = c_var[slabi] fctr = structure.film_ctr(self.var_list, self.var_table, subr, sindex, self.iq, self.h, self.k, ['lattice_abc']) ss = np.multiply(np.mat(fctr) + substrate_ctr,np.mat(r).T)*self.intensity # calculate r square sa = np.asarray(ss).reshape(-1) ss_res = np.sum((abs(self.i)*bgm - abs(sa)*bgm)**2) varience = ss_res/ss_tot print(int(varience*1e4)/1e4) return varience # select properity optimized method for lattice_c optimize # if there are only one or two variables, step_brute method is a direct method if len(slab_list) == 1: print('lattice c optimising....\nOPT_STEP_BRUTE method is used') rec = tt.opt_step_brute(c_refine,[[0.8,1.2]],grid_size = 20) elif len(slab_list) == 2: print('lattice c optimising....\nOPT_STEP_BRUTE method is used') rec = tt.opt_step_brute(c_refine,[[0.8,1.2],[0.8,1.2]],grid_size = 10) # for larger variable number, Nelder-Mead method is more comparable elif len(slab_list) >= 3: print('lattice c optimising....\nNELDER_MEAD method is used') c0 = np.ones(len(slab_list)).tolist() resc = so.minimize(c_refine,c0,method = 'Nelder-Mead',tol=10) rec = resc.x # update subr if len(slab_list) == 1: subr[2] = rec elif len(slab_list) >= 2: for slabi in range(len(slab_list)): subr[slabi*3-1] = rec[0][slabi] # update ubr_index self.ubr_index = tc.refresh_index(self.ubr_index,subr,sindex) # plot the refined result fctr = structure.film_ctr(self.var_list, self.var_table, subr, sindex, self.iq, self.h, self.k, ['lattice_abc']) ss = np.multiply(np.mat(fctr) + substrate_ctr,np.mat(r).T) plt.plot(self.iq,np.log(abs(ss*self.intensity))) plt.plot(self.iq,np.log(np.sqrt(self.i))) def dw_refine(self, p_mask = 13, n_mask = 3, weight = 1, key_list = ['dw']): # Prepare the data mask # yin mask to mask the bragg peaks bn_mask = np.mat(tc.bragg_mask(self.iq, n_mask,1, mode = 'yin')) # yang mask to refine the signal around bragg peaks only bp_mask = np.mat(tc.bragg_mask(self.iq, p_mask,1, mode = 'yang')) bgm = np.asarray(np.multiply(bn_mask, bp_mask).T).reshape(-1)*weight # calculate ss tot.For only seveal signal points are calculate # mean value is np.sum(abs(i)*bgm)/np.sum(bgm) ss_tot = np.sum((abs(self.i)*bgm - np.sum(abs(self.i)*bgm)/np.sum(bgm))**2) # initialize the variable about the lattice slab_list = np.unique(self.var_table['slab_list']).tolist() ubr_index = tc.contral_vars(self.var_list) subr,sindex = tc.initialize_contral_vars(ubr_index, slab_list, ['dw']) # roughness cr = roughness(self.iq) r = cr.robinson_roughness(self.roughness) # substrate ctr # the used keys key_list = ['absorption', 'roughness', 'scale'] # input substrate ubr and substrate index sub_ubr,sub_index = tc.initialize_contral_vars(self.ubr_index, ['substrate'], key_list) substrate_ctr = structure.substrate_ctr(self.var_list, self.iq, sub_ubr, sub_index, self.absorption, self.h, self.h, self.scale, key_list) def d_refine(d_var): fctr = structure.film_ctr(self.var_list, self.var_table, d_var, sindex, self.iq, self.h, self.k, ['dw']) ss = np.multiply(np.mat(fctr) + substrate_ctr,np.mat(r).T)*self.intensity # calculate r square sa = np.asarray(ss).reshape(-1) ss_res = np.sum((abs(self.i)*bgm - abs(sa)*bgm)**2) varience = ss_res/ss_tot print(int(varience*1e4)/1e4) return varience d0 = subr re = so.minimize(d_refine,d0,method = 'Nelder-Mead',tol=1) self.ubr_index = tc.refresh_index(self.ubr_index,re.x,sindex) # plot the refined result fctr = structure.film_ctr(self.var_list, self.var_table, re.x, sindex, self.iq, self.h, self.k, ['dw']) ss = np.multiply(np.mat(fctr) + substrate_ctr,np.mat(r).T) plt.plot(self.iq,np.log(abs(ss*self.intensity))) plt.plot(self.iq,np.log(np.sqrt(self.i))) # return the updated ubr_index, some parameters are list in ubr_index def update_var(self): # pl--parameters location in excel file pl = {'absorption':'B3', 'layers_num':'B4', 'lattice':'B6:G6', 'lattice_abc':'B6:D6', 'ocu':'B10:F10', 'dw':'B12:F14', 'dis':'B16:F18', 'pos':'B22:F24'} #p x--parameters xlsx workbook px = ox.load_workbook(self.para_path) # check the slab in ubr_index for slabi in self.ubr_index: # intensity is not a slab should be picked out if slabi == 'intensity': continue # start update the slab data else: # parameters loop for keyi in self.ubr_index[slabi]: # only several parameters are updated from ubr_index # scale roughness vacancy veta are not if keyi in ['scale','roughness','vacancy','beta']: pass else: # key is the parameters from ubr_index key = np.mat(self.ubr_index[slabi][keyi]) # psheet is the slabi sheet psheet = px[slabi] # cell is the cell of parameters in slabi sheet cell = psheet[pl[keyi]] # parameters absorptio is a singe value if keyi == 'absorption': # key_var is the value of parameters from var_list key_var = self.var_list[slabi][keyi] # key is the contral variable from ubr_index # key_var is the parameter value from excel file cell.value = key[0,0]*key_var elif keyi == 'lattice_abc': # print(key) # only apply lattice abc key_var = self.var_list[slabi]['lattice'] # print(key_var) ci = 0 # loop cell and update value # in openpyxl, cell can only be valued in this way for row in cell: for element in row: if ci <= 2: element.value = key[0, ci]*key_var.as_matrix()[0,ci] print(element.value) ci += 1 else: pass else: # for matrix parameters key_var = self.var_list[slabi][keyi] # the row index r = 0 # the column index c = 0 for row in cell: for element in row: # apply the contral variable to origin data element.value = key[r, c]*key_var.as_matrix()[r, c] # re init the column index if c == key.shape[1]-1: c = c - key.shape[1]+1 else: # loop column c += 1 # loop row r += 1 try: px.save(self.para_path) except: print('Error!!Please close the excel file %s and try again'%self.parameters) def pos_modulate(self, bc, sc, cl): subr0 = self.var_table['posz_list'] slab_index = 'posz' subr0 = relax.strain_relax(subr0, cl, bc, sc) key_list = ['pos'] print(subr0) fctr = structure.film_ctr(self.var_list, self.var_table, subr0, slab_index, self.iq, self.h, self.k, key_list) # roughness cr = roughness(self.iq) r = cr.robinson_roughness(self.roughness) # interface roughness ir = cr.interface_roughness(self.iroughness) return np.multiply(np.multiply(self.intensity*fctr, np.mat(r).T),np.mat(ir).T) def surface_modulate(self,thrface = 0.9,secface = 0.9,surface = 1.1): subr0 = self.var_table['posz_list'] slab_index = 'posz' subr0[-3] = thrface subr0[-1] = surface subr0[-2] = secface key_list = ['pos'] self.var_table['posz_list'][ 0] = thrface self.var_table['posz_list'][-1] = surface self.var_table['posz_list'][ 1] = secface # print(subr0) fctr = structure.film_ctr(self.var_list, self.var_table, subr0, slab_index, self.iq, self.h, self.k, key_list) # roughness cr = roughness(self.iq) r = cr.robinson_roughness(self.roughness) # interface roughness ir = cr.interface_roughness(self.iroughness) return np.multiply(np.multiply(self.intensity*fctr, np.mat(r).T),np.mat(ir).T) def pos_refine(self): # the used keys sub_key_list = ['absorption', 'roughness', 'scale'] # input substrate ubr and substrate index sub_ubr,sub_index = tc.initialize_contral_vars(self.ubr_index, ['substrate'], sub_key_list) substrate_ctr = structure.substrate_ctr(self.var_list, self.iq, sub_ubr, sub_index, self.absorption, self.h, self.h, self.scale, sub_key_list) # roughness cr = roughness(self.iq) r = cr.robinson_roughness(self.roughness) s = np.multiply(self.intensity*substrate_ctr, np.mat(r).T) subr0 = self.var_table['posz_list'] slab_index = 'posz' key_list = ['pos'] mask = tc.bragg_mask(self.iq,3,1,mode = 'yin',limit = 24) # fig, ax = plt.subplots() def varience(subr): fctr = structure.film_ctr(self.var_list, self.var_table, subr, slab_index, self.iq, self.h, self.k, key_list) # roughness cr = roughness(self.iq) # interface roughness ir = cr.interface_roughness(self.iroughness) f = np.multiply(np.multiply(self.intensity*fctr, np.mat(r).T),np.mat(ir).T) a = s + f e = np.sqrt(self.i) # plt.cla() # plt.plot(np.log10(e)) # plt.plot(np.log10(abs(a))) v = np.sum(abs(mask*(np.log10(abs(a)+1e-6) - np.log10(e+1e-6)))) print(v) return v re = so.minimize(varience,subr0, method = 'Nelder-Mead', options = {'maxiter':50}) subr1 = re.x fctr1 = structure.film_ctr(self.var_list, self.var_table, subr1, slab_index, self.iq, self.h, self.k, key_list) # roughness cr = roughness(self.iq) # interface roughness ir = cr.interface_roughness(self.iroughness) return np.multiply(np.multiply(self.intensity*fctr1, np.mat(r).T),np.mat(ir).T) # return the initlatised parameters def return_self(self): return self class vr(var_refine): # plot experiment data def p_rhkl(self): plt.plot(self.iq, np.log(abs(self.i) + 1e-5)) # plot sqrt experiment data def p_shkl(self): plt.plot(self.iq, np.log(np.sqrt(abs(self.i) + 1e-5)),'o') # plot substrate data def p_sctr(self): plt.plot(self.iq, np.log(abs(var_refine.sub_ctr(self)))) # plot slab data def p_fctr(self): plt.plot(self.iq, np.log(abs(var_refine.slab_ctr(self)))) # plot sumed data -- substrate data + slab data def p_actr(self): plt.plot(self.iq, np.log(abs(var_refine.sub_ctr(self) + var_refine.slab_ctr(self)))) # return the complex data of substrate, slab. The sqrt data of experiment data # compare the fitting data and experiment data def p_c(self): plt.plot(self.iq, np.log(abs(var_refine.sub_ctr(self) + var_refine.slab_ctr(self)))) plt.scatter(self.iq, np.log(np.sqrt(abs(self.i) + 1e-5)),s = 5, c = 'r') plt.plot(self.iq, np.log(np.sqrt(abs(self.i) + 1e-5)), c = 'r') # plot the pos_modulate result def p_p(self, bc, sc, cl): plt.plot(self.iq, np.log(abs(var_refine.sub_ctr(self) + var_refine.pos_modulate(self, bc, sc, cl)))) plt.scatter(self.iq, np.log(np.sqrt(abs(self.i) + 1e-5)),s = 5, c = 'r') # plot the surface modulate result def p_s(self,thrface,secface,surface): plt.plot(self.iq, np.log(abs(var_refine.sub_ctr(self) + var_refine.surface_modulate(self, thrface, secface, surface)))) plt.scatter(self.iq, np.log(np.sqrt(abs(self.i) + 1e-5)),s = 5, c = 'r') plt.plot(self.iq, np.log(np.sqrt(abs(self.i) + 1e-5)), c = 'r') # plot the pos_refine result def p_r(self): plt.plot(self.iq, np.log(abs(var_refine.sub_ctr(self) + var_refine.pos_refine(self)))) plt.scatter(self.iq, np.log(np.sqrt(abs(self.i) + 1e-5)),s = 5, c = 'r') plt.plot(self.iq, np.log(np.sqrt(abs(self.i) + 1e-5)), c = 'r') # return the fitting and sqrt experiment data def rdata(self): return {'iq':np.asarray(self.iq).reshape(-1), 'shkl':np.asarray(np.sqrt(self.i) + 1e-5).reshape(-1), 'substrate_ctr':np.asarray(var_refine.sub_ctr(self)).reshape(-1), 'slab_ctr':np.asarray(var_refine.slab_ctr(self)).reshape(-1)} # save modulated data to folder /data def save(self): ctr_optimised_result = {'q':np.asarray(self.iq).reshape(-1), 'shkl':np.asarray(np.sqrt(self.i) + 1e-5).reshape(-1), 'substrate_ctr':np.asarray(var_refine.sub_ctr(self)).reshape(-1), 'slab_ctr':np.asarray(var_refine.slab_ctr(self)).reshape(-1)} f = open(data_path, 'wb') pickle.dump(ctr_optimised_result, f) f.close()