diff --git a/src/py2dmat/solver/lib_make_convolution.py b/src/py2dmat/solver/lib_make_convolution.py index 8d591013..a4cc12c4 100644 --- a/src/py2dmat/solver/lib_make_convolution.py +++ b/src/py2dmat/solver/lib_make_convolution.py @@ -1,40 +1,50 @@ import numpy as np -import sys import copy + def calc(RC_data_org, number_of_beams, number_of_glancing_angles, omega, verbose_mode): - sigma = 0.5 * omega / (np.sqrt(2.0*np.log(2.0))) + sigma = 0.5 * omega / (np.sqrt(2.0 * np.log(2.0))) def g(x): - g = (1.0 / (sigma*np.sqrt(2.0*np.pi))) * np.exp(-0.5 * x ** 2.0 / sigma ** 2.0) + g = (1.0 / (sigma * np.sqrt(2.0 * np.pi))) * np.exp( + -0.5 * x**2.0 / sigma**2.0 + ) return g - RC_data_cnv = np.zeros((number_of_glancing_angles, number_of_beams+1)) - #copy glancing angle - RC_data_cnv[:,0] = copy.deepcopy(RC_data_org[:,0]) + RC_data_cnv = np.zeros((number_of_glancing_angles, number_of_beams + 1)) + # copy glancing angle + RC_data_cnv[:, 0] = copy.deepcopy(RC_data_org[:, 0]) if verbose_mode: - print("RC_data_org =\n",RC_data_org) - print("RC_data_cnv =\n",RC_data_cnv) - print('angle_interval=', angle_interval) + print("RC_data_org =\n", RC_data_org) + print("RC_data_cnv =\n", RC_data_cnv) for beam_index in range(number_of_beams): for index in range(number_of_glancing_angles): integral = 0.0 angle_interval = 0.0 for index2 in range(number_of_glancing_angles): - if index2 == number_of_glancing_angles - 1 : - pass - else : - angle_interval = RC_data_org[index2+1,0] - RC_data_org[index2,0] - integral += RC_data_org[index2,beam_index+1] * g(RC_data_org[index,0] - RC_data_org[index2,0]) * angle_interval + if index2 < number_of_glancing_angles - 1: + angle_interval = RC_data_org[index2 + 1, 0] - RC_data_org[index2, 0] + integral += ( + RC_data_org[index2, beam_index + 1] + * g(RC_data_org[index, 0] - RC_data_org[index2, 0]) + * angle_interval + ) if verbose_mode: - print("beam_index, index, index2, g(RC_data_org[index,0] - RC_data_org[index2,0]) =",beam_index, index, index2, g(RC_data_org[index,0] - RC_data_org[index2,0])) - RC_data_cnv[index, beam_index+1]=integral + print( + "beam_index, index, index2, g(RC_data_org[index,0] - RC_data_org[index2,0]) =", + beam_index, + index, + index2, + g(RC_data_org[index, 0] - RC_data_org[index2, 0]), + ) + RC_data_cnv[index, beam_index + 1] = integral if verbose_mode: - print("RC_data_cnv =\n",RC_data_cnv) + print("RC_data_cnv =\n", RC_data_cnv) return RC_data_cnv + diff --git a/src/py2dmat/solver/sim_trhepd_rheed.py b/src/py2dmat/solver/sim_trhepd_rheed.py index e3bc5f74..e27838bc 100644 --- a/src/py2dmat/solver/sim_trhepd_rheed.py +++ b/src/py2dmat/solver/sim_trhepd_rheed.py @@ -2,10 +2,8 @@ import os import os.path import shutil -import subprocess import time import ctypes -import copy import numpy as np @@ -13,13 +11,14 @@ from py2dmat import exception, mpi from . import lib_make_convolution -#for type hints +# for type hints from pathlib import Path from typing import List, Dict, Optional, TYPE_CHECKING if TYPE_CHECKING: from mpi4py import MPI + class Solver(py2dmat.solver.SolverBase): mpicomm: Optional["MPI.Comm"] mpisize: int @@ -28,47 +27,47 @@ class Solver(py2dmat.solver.SolverBase): isLogmode: bool detail_timer: Dict path_to_solver: Path - + def __init__(self, info: py2dmat.Info): super().__init__(info) self.mpicomm = mpi.comm() self.mpisize = mpi.size() self.mpirank = mpi.rank() - + self._name = "sim_trhepd_rheed_mb_connect" - self.run_scheme = info.solver.get("run_scheme","") - scheme_list = ["subprocess","connect_so"] + self.run_scheme = info.solver.get("run_scheme", "") + scheme_list = ["subprocess", "connect_so"] scheme_judge = [i == self.run_scheme for i in scheme_list] if not any(scheme_judge): - raise exception.InputError( - "ERROR : Input scheme is incorrect." - ) + raise exception.InputError("ERROR : Input scheme is incorrect.") if self.run_scheme == "connect_so": self.load_so() - + elif self.run_scheme == "subprocess": - #path to surf.exe + # path to surf.exe p2solver = info.solver["config"].get("surface_exec_file", "surf.exe") if os.path.dirname(p2solver) != "": # ignore ENV[PATH] self.path_to_solver = self.root_dir / Path(p2solver).expanduser() else: - for P in itertools.chain([self.root_dir], os.environ["PATH"].split(":")): + for P in itertools.chain( + [self.root_dir], os.environ["PATH"].split(":") + ): self.path_to_solver = Path(P) / p2solver if os.access(self.path_to_solver, mode=os.X_OK): break if not os.access(self.path_to_solver, mode=os.X_OK): raise exception.InputError(f"ERROR: solver ({p2solver}) is not found") - + self.isLogmode = False self.set_detail_timer() - self.input = Solver.Input(info,self.isLogmode,self.detail_timer) - self.output = Solver.Output(info,self.isLogmode,self.detail_timer) - + self.input = Solver.Input(info, self.isLogmode, self.detail_timer) + self.output = Solver.Output(info, self.isLogmode, self.detail_timer) + def set_detail_timer(self) -> None: # TODO: Operate log_mode with toml file. Generate txt of detail_timer. if self.isLogmode: @@ -106,55 +105,55 @@ def prepare(self, message: py2dmat.Message) -> None: def get_results(self) -> float: return self.output.get_results(self.work_dir) - + def load_so(self) -> None: - self.lib = np.ctypeslib.load_library("surf.so", - os.path.dirname(__file__)) + self.lib = np.ctypeslib.load_library("surf.so", os.path.dirname(__file__)) self.lib.surf_so.argtypes = ( - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_int), - np.ctypeslib.ndpointer(), - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_int), - np.ctypeslib.ndpointer(), - np.ctypeslib.ndpointer() - ) + ctypes.POINTER(ctypes.c_int), + ctypes.POINTER(ctypes.c_int), + np.ctypeslib.ndpointer(), + ctypes.POINTER(ctypes.c_int), + ctypes.POINTER(ctypes.c_int), + np.ctypeslib.ndpointer(), + np.ctypeslib.ndpointer(), + ) self.lib.surf_so.restype = ctypes.c_void_p - + def launch_so(self) -> None: - + n_template_file = len(self.input.template_file) - m_template_file = self.input.surf_tempalte_width_for_fortran + m_template_file = self.input.surf_template_width_for_fortran n_bulk_file = len(self.input.bulk_file) m_bulk_file = self.input.bulk_out_width_for_fortran - + # NOTE: The "20480" is related to the following directive in surf_so.f90. # character(c_char), intent(inout) :: surf_out(20480) - emp_str = ' '*20480 + emp_str = " " * 20480 self.output.surf_output = np.array([emp_str.encode()]) self.lib.surf_so( - ctypes.byref(ctypes.c_int(n_template_file)), - ctypes.byref(ctypes.c_int(m_template_file)), - self.input.template_file, - ctypes.byref(ctypes.c_int(n_bulk_file)), - ctypes.byref(ctypes.c_int(m_bulk_file)), - self.input.bulk_file, - self.output.surf_output - ) + ctypes.byref(ctypes.c_int(n_template_file)), + ctypes.byref(ctypes.c_int(m_template_file)), + self.input.template_file, + ctypes.byref(ctypes.c_int(n_bulk_file)), + ctypes.byref(ctypes.c_int(m_bulk_file)), + self.input.bulk_file, + self.output.surf_output, + ) self.output.surf_output = self.output.surf_output[0].decode().splitlines() - + def run(self, nprocs: int = 1, nthreads: int = 1) -> None: - if self.isLogmode : time_sta = time.perf_counter() - + if self.isLogmode: + time_sta = time.perf_counter() + if self.run_scheme == "connect_so": self.launch_so() elif self.run_scheme == "subprocess": self._run_by_subprocess([str(self.path_to_solver)]) - + if self.isLogmode: time_end = time.perf_counter() self.detail_timer["launch_STR"] += time_end - time_sta - + class Input(object): mpicomm: Optional["MPI.Comm"] mpisize: int @@ -175,7 +174,7 @@ def __init__(self, info, isLogmode, detail_timer): self.mpicomm = mpi.comm() self.mpisize = mpi.size() self.mpirank = mpi.rank() - + self.isLogmode = isLogmode self.detail_timer = detail_timer @@ -186,19 +185,19 @@ def __init__(self, info, isLogmode, detail_timer): self.dimension = info.solver["dimension"] else: self.dimension = info.base["dimension"] - - #read info + + # read info info_s = info.solver self.run_scheme = info_s["run_scheme"] self.generate_rocking_curve = info_s.get("generate_rocking_curve", False) - + # NOTE: - # surf_tempalte_width_for_fortran: Number of strings per line of template.txt data for surf.so. + # surf_template_width_for_fortran: Number of strings per line of template.txt data for surf.so. # bulk_out_width_for_fortran: Number of strings per line of bulkP.txt data for surf.so. - if self.run_scheme=="connect_so": - self.surf_tempalte_width_for_fortran = 128 + if self.run_scheme == "connect_so": + self.surf_template_width_for_fortran = 128 self.bulk_out_width_for_fortran = 1024 - + info_param = info_s.get("param", {}) v = info_param.setdefault("string_list", ["value_01", "value_02"]) if len(v) != self.dimension: @@ -219,15 +218,15 @@ def __init__(self, info, isLogmode, detail_timer): raise exception.InputError( f"ERROR: surface_template_file ({self.surface_template_file}) does not exist" ) - + if self.mpirank == 0: self._check_template() - temp_origin = self.load_surface_template_file(filename) + temp_origin = self.load_surface_template_file(filename) else: temp_origin = None - self.template_file_origin = self.mpicomm.bcast(temp_origin,root=0) + self.template_file_origin = self.mpicomm.bcast(temp_origin, root=0) - if self.run_scheme == "connect_so": + if self.run_scheme == "connect_so": filename = info_config.get("bulk_output_file", "bulkP.txt") filename = Path(filename).expanduser().resolve() self.bulk_output_file = self.root_dir / filename @@ -235,12 +234,12 @@ def __init__(self, info, isLogmode, detail_timer): raise exception.InputError( f"ERROR: bulk_output_file ({self.bulk_output_file}) does not exist" ) - + if self.mpirank == 0: bulk_f = self.load_bulk_output_file(filename) else: bulk_f = None - self.bulk_file = self.mpicomm.bcast(bulk_f,root=0) + self.bulk_file = self.mpicomm.bcast(bulk_f, root=0) else: filename = info_config.get("bulk_output_file", "bulkP.b") @@ -251,14 +250,14 @@ def __init__(self, info, isLogmode, detail_timer): f"ERROR: bulk_output_file ({self.bulk_output_file}) does not exist" ) - def load_surface_template_file(self, filename) : + def load_surface_template_file(self, filename): template_file = [] with open(self.surface_template_file) as f: for line in f: template_file.append(line) return template_file - def load_bulk_output_file(self, filename) : + def load_bulk_output_file(self, filename): bulk_file = [] with open(self.bulk_output_file) as f: for line in f: @@ -269,12 +268,13 @@ def load_bulk_output_file(self, filename) : return bulk_f def prepare(self, message: py2dmat.Message): - if self.isLogmode : time_sta = time.perf_counter() - + if self.isLogmode: + time_sta = time.perf_counter() + x_list = message.x step = message.step iset = message.set - + dimension = self.dimension string_list = self.string_list bulk_output_file = self.bulk_output_file @@ -285,33 +285,36 @@ def prepare(self, message: py2dmat.Message): fitted_value += format(x_list[index], ".8f") fitted_value = fitted_value[: len(string_list[index])] fitted_x_list.append(fitted_value) - - if self.isLogmode : - time_end = time.perf_counter() + + if self.isLogmode: + time_end = time.perf_counter() self.detail_timer["make_surf_input"] += time_end - time_sta - - if self.isLogmode : time_sta = time.perf_counter() - + + if self.isLogmode: + time_sta = time.perf_counter() + + folder_name = "." if self.generate_rocking_curve: folder_name = self._pre_bulk(step, bulk_output_file, iset) else: if self.run_scheme == "connect_so": folder_name = "." - + elif self.run_scheme == "subprocess": - #make workdir and copy bulk output file + # make workdir and copy bulk output file folder_name = self._pre_bulk(step, bulk_output_file, iset) - - if self.isLogmode : - time_end = time.perf_counter() + + if self.isLogmode: + time_end = time.perf_counter() self.detail_timer["prepare_Log-directory"] += time_end - time_sta - if self.isLogmode : time_sta = time.perf_counter() - + if self.isLogmode: + time_sta = time.perf_counter() + self._replace(fitted_x_list, folder_name) - - if self.isLogmode : - time_end = time.perf_counter() + + if self.isLogmode: + time_end = time.perf_counter() self.detail_timer["make_surf_input"] += time_end - time_sta return fitted_x_list, folder_name @@ -321,7 +324,7 @@ def _pre_bulk(self, Log_number, bulk_output_file, iset): os.makedirs(folder_name, exist_ok=True) if self.run_scheme == "connect_so": pass - else: #self.run_scheme == "subprocess": + else: # self.run_scheme == "subprocess": shutil.copy( bulk_output_file, os.path.join(folder_name, bulk_output_file.name) ) @@ -330,24 +333,24 @@ def _pre_bulk(self, Log_number, bulk_output_file, iset): def _replace(self, fitted_x_list, folder_name): template_file = [] if self.run_scheme == "subprocess": - file_output = open(os.path.join(folder_name, self.surface_input_file), "w") + file_output = open( + os.path.join(folder_name, self.surface_input_file), "w" + ) for line in self.template_file_origin: for index in range(self.dimension): if line.find(self.string_list[index]) != -1: line = line.replace( - self.string_list[index], - fitted_x_list[index] - ) + self.string_list[index], fitted_x_list[index] + ) if self.run_scheme == "connect_so": line = line.replace("\t", " ").replace("\n", " ") - line = line.encode().ljust(self.surf_tempalte_width_for_fortran) + line = line.encode().ljust(self.surf_template_width_for_fortran) template_file.append(line) - + elif self.run_scheme == "subprocess": file_output.write(line) - if self.run_scheme == "connect_so": self.template_file = np.array(template_file) elif self.run_scheme == "subprocess": @@ -372,6 +375,7 @@ class Output(object): """ Output manager. """ + mpicomm: Optional["MPI.Comm"] mpisize: int mpirank: int @@ -412,9 +416,9 @@ def __init__(self, info, isLogmode, detail_timer): info_s = info.solver self.run_scheme = info_s["run_scheme"] - - #If self.run_scheme == "connect_so", - #the contnts of surface_output_file are retailned in self.surf_output. + + # If self.run_scheme == "connect_so", + # the contents of surface_output_file are retailned in self.surf_output. self.surf_output = np.array([]) self.generate_rocking_curve = info_s.get("generate_rocking_curve", False) @@ -423,58 +427,60 @@ def __init__(self, info, isLogmode, detail_timer): v = info_post.get("normalization", "") available_normalization = ["TOTAL", "MANY_BEAM"] if v == "MAX": - raise exception.InputError('ERROR: normalization == "MAX" is not available') + raise exception.InputError( + 'ERROR: normalization == "MAX" is not available' + ) if v not in available_normalization: - msg ="ERROR: normalization must be " - msg+="MANY_BEAM or TOTAL" + msg = "ERROR: normalization must be " + msg += "MANY_BEAM or TOTAL" raise exception.InputError(msg) self.normalization = v v = info_post.get("weight_type", None) available_weight_type = ["calc", "manual"] if self.normalization == "MANY_BEAM": - if v is None : - msg ='ERROR: If normalization = "MANY_BEAM", ' - msg+='"weight_type" must be set in [solver.post].' + if v is None: + msg = 'ERROR: If normalization = "MANY_BEAM", ' + msg += '"weight_type" must be set in [solver.post].' raise exception.InputError(msg) elif v not in available_weight_type: - msg ="ERROR: weight_type must be " - msg+="calc or manual" + msg = "ERROR: weight_type must be " + msg += "calc or manual" raise exception.InputError(msg) else: - if v is not None : + if v is not None: if self.mpirank == 0: - msg ='NOTICE: If normalization = "MANY_BEAM" is not set, ' - msg+='"weight_type" is NOT considered in the calculation.' - print(msg) + msg = 'NOTICE: If normalization = "MANY_BEAM" is not set, ' + msg += '"weight_type" is NOT considered in the calculation.' + print(msg) self.weight_type = None self.weight_type = v v = info_post.get("spot_weight", []) - if self.normalization=="MANY_BEAM" and self.weight_type=="manual": - if v == [] : - msg ='ERROR: With normalization="MANY_BEAM" and ' - msg+='weight_type=="manual", the "spot_weight" option ' - msg+='must be set in [solver.post].' + if self.normalization == "MANY_BEAM" and self.weight_type == "manual": + if v == []: + msg = 'ERROR: With normalization="MANY_BEAM" and ' + msg += 'weight_type=="manual", the "spot_weight" option ' + msg += "must be set in [solver.post]." raise exception.InputError(msg) - self.spot_weight = np.array(v)/sum(v) + self.spot_weight = np.array(v) / sum(v) else: if v != []: if self.mpirank == 0: - msg ='NOTICE: With the specified "normalization" option, ' - msg+='the "spot_weight" you specify in the toml file is ignored, ' - msg+='since the spot_weight is automatically calculated within the program.' + msg = 'NOTICE: With the specified "normalization" option, ' + msg += 'the "spot_weight" you specify in the toml file is ignored, ' + msg += "since the spot_weight is automatically calculated within the program." print(msg) - if self.normalization== "TOTAL": + if self.normalization == "TOTAL": self.spot_weight = np.array([1]) v = info_post.get("Rfactor_type", "A") if v not in ["A", "B", "A2"]: raise exception.InputError("ERROR: Rfactor_type must be A, A2 or B") - if self.normalization=="MANY_BEAM": - if (v!="A") and (v!="A2") : - msg ='With normalization="MANY_BEAM", ' - msg+='only Rfactor_type="A" or Rfactor_type="A2" is valid.' + if self.normalization == "MANY_BEAM": + if (v != "A") and (v != "A2"): + msg = 'With normalization="MANY_BEAM", ' + msg += 'only Rfactor_type="A" or Rfactor_type="A2" is valid.' raise exception.InputError(msg) self.Rfactor_type = v @@ -482,10 +488,9 @@ def __init__(self, info, isLogmode, detail_timer): if v <= 0.0: raise exception.InputError("ERROR: omega should be positive") self.omega = v - + self.remove_work_dir = info_post.get("remove_work_dir", False) - - + # solver.param info_param = info_s.get("param", {}) v = info_param.setdefault("string_list", ["value_01", "value_02"]) @@ -505,7 +510,7 @@ def __init__(self, info, isLogmode, detail_timer): firstline = v # None is dummy value - # If "reference_last_line" is not specified in the toml file, + # If "reference_last_line" is not specified in the toml file, # the last line of the reference file is used for the R-factor calculation. v = info_ref.setdefault("reference_last_line", None) if v is None: @@ -520,83 +525,79 @@ def __init__(self, info, isLogmode, detail_timer): reference_path = info_ref.get("path", "experiment.txt") data_experiment = self.read_experiment( - reference_path, firstline, lastline, - reference_are_read_to_final_line ) + reference_path, firstline, lastline, reference_are_read_to_final_line + ) self.angle_number_experiment = data_experiment.shape[0] self.beam_number_exp_raw = data_experiment.shape[1] - + v = info_ref.get("exp_number", []) - - if len(v) == 0 : - raise exception.InputError( - "ERROR: You have to set the 'exp_number'." - ) + + if len(v) == 0: + raise exception.InputError("ERROR: You have to set the 'exp_number'.") if not isinstance(v, list): - raise exception.InputError( - "ERROR: 'exp_number' must be a list type." - ) + raise exception.InputError("ERROR: 'exp_number' must be a list type.") - if max(v) > self.beam_number_exp_raw : - raise exception.InputError( - "ERROR: The 'exp_number' setting is wrong." - ) + if max(v) > self.beam_number_exp_raw: + raise exception.InputError("ERROR: The 'exp_number' setting is wrong.") - if self.normalization=="MANY_BEAM" and self.weight_type=="manual": + if self.normalization == "MANY_BEAM" and self.weight_type == "manual": if len(v) != len(self.spot_weight): raise exception.InputError( - "ERROR:len('exp_number') and len('spot_weight') do not match." + "ERROR:len('exp_number') and len('spot_weight') do not match." ) - if self.normalization=="TOTAL" and len(v)!=1: - msg ='When normalization=="TOTAL" is specified, ' - msg+='only one beam data can be specified with ' - msg+='"exp_number" option.' + if self.normalization == "TOTAL" and len(v) != 1: + msg = 'When normalization=="TOTAL" is specified, ' + msg += "only one beam data can be specified with " + msg += '"exp_number" option.' raise exception.InputError(msg) self.exp_number = v - - #Normalization of reference data + + # Normalization of reference data self.beam_number_experiment = len(self.exp_number) for loop_index in range(self.beam_number_experiment): exp_index = self.exp_number[loop_index] - I_reference = data_experiment[:,exp_index] + I_reference = data_experiment[:, exp_index] if self.normalization == "TOTAL": I_reference_norm = np.sum(I_reference) - I_reference_normalized = I_reference/I_reference_norm + I_reference_normalized = I_reference / I_reference_norm I_reference_norm_l = np.array([I_reference_norm]) self.I_reference_normalized_l = np.array([I_reference_normalized]) - elif self.normalization=="MANY_BEAM" and self.weight_type=="calc": + elif self.normalization == "MANY_BEAM" and self.weight_type == "calc": I_reference_norm = np.sum(I_reference) - I_reference_normalized = I_reference/I_reference_norm - if loop_index == 0: #first loop + I_reference_normalized = I_reference / I_reference_norm + if loop_index == 0: # first loop I_reference_norm_l = np.array([I_reference_norm]) - self.I_reference_normalized_l = np.array([I_reference_normalized]) - else : # N-th loop + self.I_reference_normalized_l = np.array( + [I_reference_normalized] + ) + else: # N-th loop I_reference_norm_l = np.block( - [I_reference_norm_l, I_reference_norm] - ) + [I_reference_norm_l, I_reference_norm] + ) self.I_reference_normalized_l = np.block( - [[self.I_reference_normalized_l], - [I_reference_normalized]] - ) - elif self.normalization=="MANY_BEAM" and self.weight_type=="manual": + [[self.I_reference_normalized_l], [I_reference_normalized]] + ) + elif self.normalization == "MANY_BEAM" and self.weight_type == "manual": I_reference_norm = np.sum(I_reference) - I_reference_normalized = I_reference/I_reference_norm - if loop_index == 0: #first loop + I_reference_normalized = I_reference / I_reference_norm + if loop_index == 0: # first loop I_reference_norm_l = np.array([I_reference_norm]) - self.I_reference_normalized_l = np.array([I_reference_normalized]) - else : # N-th loop + self.I_reference_normalized_l = np.array( + [I_reference_normalized] + ) + else: # N-th loop I_reference_norm_l = np.block( [I_reference_norm_l, I_reference_norm] - ) + ) self.I_reference_normalized_l = np.block( - [[self.I_reference_normalized_l], - [I_reference_normalized]] - ) + [[self.I_reference_normalized_l], [I_reference_normalized]] + ) else: - msg ="ERROR: normalization must be " - msg+="MANY_BEAM or TOTAL" + msg = "ERROR: normalization must be " + msg += "MANY_BEAM or TOTAL" raise exception.InputError(msg) # solver.config @@ -613,31 +614,32 @@ def __init__(self, info, isLogmode, detail_timer): self.calculated_first_line = v v = info_config.get( - "calculated_last_line", - self.calculated_first_line + self.angle_number_experiment-1 - ) + "calculated_last_line", + self.calculated_first_line + self.angle_number_experiment - 1, + ) if not (isinstance(v, int) and v >= 0): raise exception.InputError( "ERROR: calculated_last_line should be non-negative integer" ) self.calculated_last_line = v - - # Number of lines in the computation file + + # Number of lines in the computation file # used for R-factor calculations. - self.calculated_nlines = (self.calculated_last_line - - self.calculated_first_line + 1 ) - - if self.angle_number_experiment != self.calculated_nlines : + self.calculated_nlines = ( + self.calculated_last_line - self.calculated_first_line + 1 + ) + + if self.angle_number_experiment != self.calculated_nlines: raise exception.InputError( "ERROR: The number of glancing angles in the calculation data does not match the number of glancing angles in the experimental data." ) - - # Variable indicating whether the warning - # "self.calculated_nlines does not match - # the number of glancing angles in the calculated file" + + # Variable indicating whether the warning + # "self.calculated_nlines does not match + # the number of glancing angles in the calculated file" # is displayed. self.isWarning_calcnline = False - + v = info_config.get("calculated_info_line", 2) if not (isinstance(v, int) and v >= 0): raise exception.InputError( @@ -645,26 +647,22 @@ def __init__(self, info, isLogmode, detail_timer): ) self.calculated_info_line = v - v = info_config.get("cal_number",[]) - if len(v) == 0 : - raise exception.InputError( - "ERROR: You have to set the 'cal_number'." - ) - + v = info_config.get("cal_number", []) + if len(v) == 0: + raise exception.InputError("ERROR: You have to set the 'cal_number'.") + if not isinstance(v, list): - raise exception.InputError( - "ERROR: 'cal_number' must be a list type." - ) - - if self.normalization=="MANY_BEAM" and self.weight_type=="manual": + raise exception.InputError("ERROR: 'cal_number' must be a list type.") + + if self.normalization == "MANY_BEAM" and self.weight_type == "manual": if len(self.spot_weight) != len(v): raise exception.InputError( "len('cal_number') and len('spot_weight') do not match." ) - if self.normalization=="TOTAL" and len(v)!=1: - msg ='When normalization=="TOTAL" is specified, ' - msg+='only one beam data can be specified with ' - msg+='"cal_number" option.' + if self.normalization == "TOTAL" and len(v) != 1: + msg = 'When normalization=="TOTAL" is specified, ' + msg += "only one beam data can be specified with " + msg += '"cal_number" option.' raise exception.InputError(msg) self.cal_number = v @@ -675,28 +673,28 @@ def read_experiment(self, ref_path, first, last, read_to_final_line): Elines = file_input.readlines() firstline = first - if read_to_final_line : + if read_to_final_line: lastline = len(Elines) else: lastline = last - - n_exp_row = lastline-firstline+1 - + + n_exp_row = lastline - firstline + 1 + # get value from data - for row_index in range(n_exp_row) : + for row_index in range(n_exp_row): line_index = firstline + row_index - 1 line = Elines[line_index] data = line.split() # first loop if row_index == 0: n_exp_column = len(data) - data_e = np.zeros((n_exp_row, n_exp_column)) #init value - + data_e = np.zeros((n_exp_row, n_exp_column)) # init value + for column_index in range(n_exp_column): data_e[row_index, column_index] = float(data[column_index]) else: data_e = None - data_exp = self.mpicomm.bcast(data_e,root=0) + data_exp = self.mpicomm.bcast(data_e, root=0) return data_exp def prepare(self, fitted_x_list): @@ -707,57 +705,62 @@ def get_results(self, work_dir) -> float: Get Rfactor obtained by the solver program. Returns ------- - """ + """ # Calculate Rfactor and Output numerical results cwd = os.getcwd() os.chdir(work_dir) Rfactor = self._post(self.fitted_x_list) os.chdir(cwd) - #delete Log-directory - if self.isLogmode : time_sta = time.perf_counter() - + # delete Log-directory + if self.isLogmode: + time_sta = time.perf_counter() + if self.run_scheme == "subprocess": if self.remove_work_dir: + def rmtree_error_handler(function, path, excinfo): print(f"WARNING: Failed to remove a working directory, {path}") + shutil.rmtree(work_dir, onerror=rmtree_error_handler) - - if self.isLogmode : + + if self.isLogmode: time_end = time.perf_counter() - self.detail_timer["delete_Log-directory"] += time_end - time_sta + self.detail_timer["delete_Log-directory"] += time_end - time_sta return Rfactor def _post(self, fitted_x_list): I_experiment_normalized_l = self.I_reference_normalized_l - + ( glancing_angle, conv_number_of_g_angle, conv_I_calculated_norm_l, conv_I_calculated_normalized_l, ) = self._calc_I_from_file() - - if self.isLogmode : time_sta = time.perf_counter() + + if self.isLogmode: + time_sta = time.perf_counter() Rfactor = self._calc_Rfactor( - conv_number_of_g_angle, - conv_I_calculated_normalized_l, - I_experiment_normalized_l, - ) - if self.isLogmode : + conv_number_of_g_angle, + conv_I_calculated_normalized_l, + I_experiment_normalized_l, + ) + if self.isLogmode: time_end = time.perf_counter() self.detail_timer["calculate_R-factor"] += time_end - time_sta - - #generate RockingCurve_calculated.txt + + # generate RockingCurve_calculated.txt dimension = self.dimension string_list = self.string_list cal_number = self.cal_number spot_weight = self.spot_weight - weight_type = self.weight_type + weight_type = self.weight_type Rfactor_type = self.Rfactor_type normalization = self.normalization - if self.generate_rocking_curve : - if self.isLogmode : time_sta = time.perf_counter() + if self.generate_rocking_curve: + if self.isLogmode: + time_sta = time.perf_counter() with open("RockingCurve_calculated.txt", "w") as file_RC: # Write headers file_RC.write("#") @@ -773,13 +776,17 @@ def _post(self, fitted_x_list): file_RC.write("#fx(x) = {}\n".format(Rfactor)) file_RC.write("#cal_number = {}\n".format(cal_number)) file_RC.write("#spot_weight = {}\n".format(spot_weight)) - file_RC.write("#NOTICE : Intensities are NOT multiplied by spot_weight.\n") - file_RC.write("#The intensity I_(spot) for each spot is normalized as in the following equation.\n") + file_RC.write( + "#NOTICE : Intensities are NOT multiplied by spot_weight.\n" + ) + file_RC.write( + "#The intensity I_(spot) for each spot is normalized as in the following equation.\n" + ) file_RC.write("#sum( I_(spot) ) = 1\n") file_RC.write("#\n") - + label_column = ["glancing_angle"] - fmt_rc = '%.5f' + fmt_rc = "%.5f" for i in range(len(cal_number)): label_column.append(f"cal_number={self.cal_number[i]}") fmt_rc += " %.15e" @@ -789,22 +796,20 @@ def _post(self, fitted_x_list): file_RC.write("\n") g_angle_for_rc = np.array([glancing_angle]) np.savetxt( - file_RC, - np.block( - [g_angle_for_rc.T, - conv_I_calculated_normalized_l.T] - ), - fmt=fmt_rc - ) - - if self.isLogmode : + file_RC, + np.block([g_angle_for_rc.T, conv_I_calculated_normalized_l.T]), + fmt=fmt_rc, + ) + + if self.isLogmode: time_end = time.perf_counter() - self.detail_timer["make_RockingCurve.txt"] += time_end - time_sta + self.detail_timer["make_RockingCurve.txt"] += time_end - time_sta return Rfactor def _calc_I_from_file(self): - if self.isLogmode : time_sta = time.perf_counter() - + if self.isLogmode: + time_sta = time.perf_counter() + surface_output_file = self.surface_output_file calculated_first_line = self.calculated_first_line calculated_last_line = self.calculated_last_line @@ -818,146 +823,155 @@ def _calc_I_from_file(self): if self.run_scheme == "connect_so": Clines = self.surf_output - + elif self.run_scheme == "subprocess": file_input = open(self.surface_output_file, "r") Clines = file_input.readlines() file_input.close() # Reads the number of glancing angles and beams - line = Clines[calculated_info_line-1] + line = Clines[calculated_info_line - 1] line = line.replace(",", "") data = line.split() - + calc_number_of_g_angles_org = int(data[1]) - calc_number_of_beams_org = int(data[2]) - + calc_number_of_beams_org = int(data[2]) + if calc_number_of_g_angles_org != calculated_nlines: - if self.mpirank == 0 and not self.isWarning_calcnline : - msg ="WARNING:\n" - msg+="The number of lines in the calculated file " - msg+="that you have set up does not match " - msg+="the number of glancing angles in the calculated file.\n" - msg+="The number of lines (nline) in the calculated file " - msg+="used for the R-factor calculation " - msg+="is determined by the following values\n" - msg+='nline = "calculated_last_line" - "calculated_first_line" + 1.' + if self.mpirank == 0 and not self.isWarning_calcnline: + msg = "WARNING:\n" + msg += "The number of lines in the calculated file " + msg += "that you have set up does not match " + msg += "the number of glancing angles in the calculated file.\n" + msg += "The number of lines (nline) in the calculated file " + msg += "used for the R-factor calculation " + msg += "is determined by the following values\n" + msg += ( + 'nline = "calculated_last_line" - "calculated_first_line" + 1.' + ) print(msg) self.isWarning_calcnline = True calc_number_of_g_angles = calculated_nlines - + # Define the array for the original calculated data. - RC_data_org = np.zeros((calc_number_of_g_angles, calc_number_of_beams_org+1)) + RC_data_org = np.zeros( + (calc_number_of_g_angles, calc_number_of_beams_org + 1) + ) for g_angle_index in range(calc_number_of_g_angles): line_index = (calculated_first_line - 1) + g_angle_index - line = Clines[ line_index ] + line = Clines[line_index] line = line.replace(",", "") data = line.split() - RC_data_org[g_angle_index,0]=float(data[0]) + RC_data_org[g_angle_index, 0] = float(data[0]) for beam_index in range(calc_number_of_beams_org): - RC_data_org[g_angle_index, beam_index+1] = data[beam_index+1] + RC_data_org[g_angle_index, beam_index + 1] = data[beam_index + 1] - if self.isLogmode : + if self.isLogmode: time_end = time.perf_counter() self.detail_timer["load_STR_result"] += time_end - time_sta - - if self.isLogmode : time_sta = time.perf_counter() + + if self.isLogmode: + time_sta = time.perf_counter() verbose_mode = False # convolution data_convolution = lib_make_convolution.calc( - RC_data_org, - calc_number_of_beams_org, - calc_number_of_g_angles, - omega, - verbose_mode, - ) - - if self.isLogmode : + RC_data_org, + calc_number_of_beams_org, + calc_number_of_g_angles, + omega, + verbose_mode, + ) + + if self.isLogmode: time_end = time.perf_counter() self.detail_timer["convolution"] += time_end - time_sta - if self.isLogmode : time_sta = time.perf_counter() - + if self.isLogmode: + time_sta = time.perf_counter() + angle_number_convolution = data_convolution.shape[0] - if self.angle_number_experiment !=angle_number_convolution: + if self.angle_number_experiment != angle_number_convolution: raise exception.InputError( "ERROR: The number of glancing angles in the calculation data does not match the number of glancing angles in the experimental data." ) - glancing_angle = data_convolution[:,0] - + glancing_angle = data_convolution[:, 0] + beam_number_reference = len(cal_number) # Normalization of calculated data. for loop_index in range(beam_number_reference): cal_index = cal_number[loop_index] - conv_I_calculated = data_convolution[:,cal_index] + conv_I_calculated = data_convolution[:, cal_index] if self.normalization == "TOTAL": conv_I_calculated_norm = np.sum(conv_I_calculated) - conv_I_calculated_normalized = conv_I_calculated/conv_I_calculated_norm - conv_I_calculated_norm_l = np.array( - [conv_I_calculated_norm] - ) + conv_I_calculated_normalized = ( + conv_I_calculated / conv_I_calculated_norm + ) + conv_I_calculated_norm_l = np.array([conv_I_calculated_norm]) conv_I_calculated_normalized_l = np.array( - [conv_I_calculated_normalized] - ) - elif self.normalization=="MANY_BEAM" and self.weight_type=="calc": + [conv_I_calculated_normalized] + ) + elif self.normalization == "MANY_BEAM" and self.weight_type == "calc": conv_I_calculated_norm = np.sum(conv_I_calculated) - conv_I_calculated_normalized = conv_I_calculated/conv_I_calculated_norm - if loop_index == 0: #first loop - conv_I_calculated_norm_l = np.array( - [conv_I_calculated_norm] - ) + conv_I_calculated_normalized = ( + conv_I_calculated / conv_I_calculated_norm + ) + if loop_index == 0: # first loop + conv_I_calculated_norm_l = np.array([conv_I_calculated_norm]) conv_I_calculated_normalized_l = np.array( - [conv_I_calculated_normalized] - ) - else: # N-th loop + [conv_I_calculated_normalized] + ) + else: # N-th loop conv_I_calculated_norm_l = np.block( - [conv_I_calculated_norm_l, - conv_I_calculated_norm] - ) + [conv_I_calculated_norm_l, conv_I_calculated_norm] + ) conv_I_calculated_normalized_l = np.block( - [[conv_I_calculated_normalized_l], - [conv_I_calculated_normalized]] - ) - if loop_index == beam_number_reference-1: #first loop - #calculate spot_weight - self.spot_weight = ( conv_I_calculated_norm_l - / sum(conv_I_calculated_norm_l) )**2 - elif self.normalization=="MANY_BEAM" and self.weight_type=="manual": + [ + [conv_I_calculated_normalized_l], + [conv_I_calculated_normalized], + ] + ) + if loop_index == beam_number_reference - 1: # first loop + # calculate spot_weight + self.spot_weight = ( + conv_I_calculated_norm_l / sum(conv_I_calculated_norm_l) + ) ** 2 + elif self.normalization == "MANY_BEAM" and self.weight_type == "manual": conv_I_calculated_norm = np.sum(conv_I_calculated) - conv_I_calculated_normalized = conv_I_calculated/conv_I_calculated_norm - if loop_index == 0: #first loop - conv_I_calculated_norm_l = np.array( - [conv_I_calculated_norm] - ) + conv_I_calculated_normalized = ( + conv_I_calculated / conv_I_calculated_norm + ) + if loop_index == 0: # first loop + conv_I_calculated_norm_l = np.array([conv_I_calculated_norm]) conv_I_calculated_normalized_l = np.array( - [conv_I_calculated_normalized] - ) - else: # N-th loop + [conv_I_calculated_normalized] + ) + else: # N-th loop conv_I_calculated_norm_l = np.block( - [conv_I_calculated_norm_l, - conv_I_calculated_norm] - ) + [conv_I_calculated_norm_l, conv_I_calculated_norm] + ) conv_I_calculated_normalized_l = np.block( - [[conv_I_calculated_normalized_l], - [conv_I_calculated_normalized]] - ) + [ + [conv_I_calculated_normalized_l], + [conv_I_calculated_normalized], + ] + ) else: - msg ="ERROR: normalization must be " - msg+="MANY_BEAM or TOTAL" + msg = "ERROR: normalization must be " + msg += "MANY_BEAM or TOTAL" raise exception.InputError(msg) - - if self.isLogmode : + + if self.isLogmode: time_end = time.perf_counter() self.detail_timer["normalize_calc_I"] += time_end - time_sta - + return ( glancing_angle, angle_number_convolution, conv_I_calculated_norm_l, - conv_I_calculated_normalized_l - ) - + conv_I_calculated_normalized_l, + ) + def _calc_Rfactor(self, n_g_angle, calc_result, exp_result): spot_weight = self.spot_weight n_spot = len(spot_weight) @@ -966,35 +980,37 @@ def _calc_Rfactor(self, n_g_angle, calc_result, exp_result): for spot_index in range(n_spot): R_spot = 0.0 for angle_index in range(n_g_angle): - R_spot += (exp_result[spot_index, angle_index] - - calc_result[spot_index, angle_index] )**2 + R_spot += ( + exp_result[spot_index, angle_index] + - calc_result[spot_index, angle_index] + ) ** 2 R_spot = spot_weight[spot_index] * R_spot - R += R_spot + R += R_spot R = np.sqrt(R) elif self.Rfactor_type == "A2": R = 0.0 for spot_index in range(n_spot): R_spot = 0.0 for angle_index in range(n_g_angle): - R_spot += (exp_result[spot_index, angle_index] - - calc_result[spot_index, angle_index] )**2 + R_spot += ( + exp_result[spot_index, angle_index] + - calc_result[spot_index, angle_index] + ) ** 2 R_spot = spot_weight[spot_index] * R_spot - R += R_spot + R += R_spot else: # self.Rfactor_type == "B" - all_exp_result = [] + all_exp_result = [] all_calc_result = [] for spot_index in range(n_spot): for angle_index in range(n_g_angle): - all_exp_result.append( - exp_result[spot_index, angle_index]) - all_calc_result.append( - calc_result[spot_index, angle_index]) + all_exp_result.append(exp_result[spot_index, angle_index]) + all_calc_result.append(calc_result[spot_index, angle_index]) y1 = 0.0 y2 = 0.0 y3 = 0.0 for I_exp, I_calc in zip(all_exp_result, all_calc_result): y1 += (I_exp - I_calc) ** 2 - y2 += I_exp ** 2 - y3 += I_calc ** 2 + y2 += I_exp**2 + y3 += I_calc**2 R = y1 / (y2 + y3) return R