Implement filtering of multiple files to the same wavelength scale (db9a8d6b) · Commits · Giacomo Mulas / NP_TMcode

src/scripts/filter_constants.py

+165 −47

Original line number	Diff line number	Diff line
		@@ -27,7 +27,8 @@
		# The script requires python3.

		import math
		#import pdb
		import numpy as np
		import pdb
		from sys import argv

		## \cond
		@@ -86,10 +87,43 @@ def main():
		print(msg_distance1)
		print(msg_distance2)
		if (config['make_plots']):
		plot_data(scan_info, config['wl_units'])
		plot_data([scan_info], config['wl_units'])
		# end result == 0 check
		else:
		result = scan_multiple_files(config)
		split_output = config['output_files'].split(',')
		aligned_infos = match_multiple_files(config)
		for i in range(len(aligned_infos)):
		info = aligned_infos[i]
		ecode = info['exit_code']
		file_name = split_output[i]
		if (ecode == 0):
		write_output_file(file_name, info)
		# INFO section
		num_filtered_data = info['num_filtered_data']
		num_orig_data = info['num_orig_data']
		msg_filtered = "{0:d} filtered lines".format(num_filtered_data) if num_filtered_data != 1 else "1 filtered line"
		msg_orig = "{0:d} input lines".format(num_orig_data) if num_orig_data != 1 else "1 input line"
		max_distance = get_max_distance(info)
		msg_distance1 = "INFO: maximum absolute distance was {0:.5g} at {1:.5e} {2:s} in the {3:s} part".format(
		max_distance['max_difference'],
		max_distance['max_wavelength'],
		config['wl_units'],
		max_distance['differing_set']
		)
		msg_distance2 = " (fitted value is {0:.5g}, actual data value is {1:.5g}).".format(
		max_distance['max_fitted'],
		max_distance['max_value']
		)
		print("INFO: extracted %s out of %s."%(msg_filtered, msg_orig))
		print(msg_distance1)
		print(msg_distance2)
		# end INFO section
		else:
		print("WARNING: scanning {0:s} resulted in error code {1:d}.".format(file_name, ecode))
		result += ecode
		if (result == 0):
		if (config['make_plots']):
		plot_data(aligned_infos, config['wl_units'])
		return result

		## \brief Maximum distance between input data and linear interpolation of filtered data.
		@@ -180,6 +214,7 @@ def parse_arguments():
		'threshold': 0.1,
		'wl_start': 0.0,
		'wl_end': 0.0,
		'wl_tolerance': 0.0,
		'wl_units': 'micrometers'
		}
		skip_arg = False
		@@ -217,6 +252,9 @@ def parse_arguments():
		elif (arg.startswith("--wl-stop=")):
		split_arg = arg.split('=')
		config['wl_end'] = float(split_arg[1])
		elif (arg.startswith("--wl-tol=")):
		split_arg = arg.split('=')
		config['wl_tolerance'] = float(split_arg[1])
		elif (arg.startswith("--wl-units=")):
		known_units = [
		'nanometers', 'nm', 'micrometers', 'um', 'millimeters', 'mm',
		@@ -251,27 +289,40 @@ def parse_arguments():

		## \brief Make a quick-look plot with MATPLOTLIB.
		#
		# \param[in] scan_info: `dict` A dictionary with the scanned data file.
		# \param[in] scan_infos: `list` A list of dictionaries with the scanned data file.
		# \param[in] wl_units: `string` The name of the units for the wavelength scale.
		def plot_data(scan_info, wl_units):
		def plot_data(scan_infos, wl_units):
		cname_array = [
		'red',
		'blue',
		'green',
		'purple',
		'orange',
		'cyan',
		'grey',
		'black'
		]
		# Plot making section
		for i in range(len(scan_infos)):
		scan_info = scan_infos[i]
		wl_orig = scan_info['wl_orig']
		reps_orig = scan_info['reps_orig']
		ieps_orig = scan_info['ieps_orig']
		wl_filtered = scan_info['wl_filtered']
		reps_filtered = scan_info['reps_filtered']
		ieps_filtered = scan_info['ieps_filtered']
		plt.plot(wl_orig, reps_orig, color='red', marker='', ls='-', label=r"Original $\mathfrak{Re}(\varepsilon)$")
		plt.plot(wl_filtered, reps_filtered, color='red', marker='o', ls='', label=r"Filtered $\mathfrak{Re}(\varepsilon)$")
		plt.plot(wl_orig, ieps_orig, color='blue', marker='', ls='-', label=r"Original $\mathfrak{Im}(\varepsilon)$")
		plt.plot(wl_filtered, ieps_filtered, color='blue', marker='o', ls='', label=r"Filtered $\mathfrak{Im}(\varepsilon)$")
		rcname = cname_array[(2 * i) % 8]
		icname = cname_array[(2 * i + 1) % 8]
		plt.plot(wl_orig, reps_orig, color=rcname, marker='', ls='-', label=r"Original $\mathfrak{Re}(\varepsilon)$")
		plt.plot(wl_filtered, reps_filtered, color=rcname, marker='o', ls='', label=r"Filtered $\mathfrak{Re}(\varepsilon)$")
		plt.plot(wl_orig, ieps_orig, color=icname, marker='', ls='--', label=r"Original $\mathfrak{Im}(\varepsilon)$")
		plt.plot(wl_filtered, ieps_filtered, color=icname, marker='s', ls='', label=r"Filtered $\mathfrak{Im}(\varepsilon)$")
		plt.xlabel("Wavelength ({0:s})".format(wl_units))
		plt.ylabel(r"$\mathfrak{Re}(\varepsilon)$\|$\mathfrak{Im}(\varepsilon)$")
		plt.legend(loc="best")
		plt.show()
		# end plot making section


		## \brief Print a command-line help summary.
		def print_help():
		print(" ############################################### ")
		@@ -303,6 +354,9 @@ def print_help():
		print(" window. ")
		print("--wl_stop=VALUE Ending wavelength in meters for the filtering ")
		print(" window. ")
		print("--wl_tol=VALUE Minimum separation to considered two wavelength ")
		print(" values as distinct in multiple files (default is")
		print(" 0.001 times the shortest wavelength).")
		print("--wl_units=UNITS Name of the wavelength units ONLY FOR PLOTTING ")
		print(" PURPOSES (data must be always in meters, only ")
		print(" MATPLOTLIB uses this setting for formatting). ")
		@@ -319,13 +373,16 @@ def print_help():
		# \param[in] config: `dict` A dictionary containing the script configuration.
		# \param[in] file_name: `string` The name of the single input file.
		# \return result: `dict` A dictionary containing the results of the scan,
		# including "exit_code" (`int`, 0 if succesful), "wl_orig" (`array-like`, the
		# original wavelength scale), "reps_orig" (`array-like`, the original real
		# parts of the dielectric functions), "ieps_orig" (`array-like`, the original
		# imaginary parts of the dielectric functions), "wl_filtered" (`array-like`,
		# the filtered wavelength scale), "reps_filtered" (`array-like`, the filtered
		# real parts of the dielectric functions), and "ieps_filtered" (`array-like`,
		# the filtered imaginary parts of the dielectric functions).
		# including "exit_code" (`int`, 0 if succesful), "wl_orig" (`array-like`,
		# the original wavelength scale), "reps_orig" (`array-like`, the original
		# real parts of the dielectric functions), "ieps_orig" (`array-like`, the
		# original imaginary parts of the dielectric functions), "wl_filtered"
		# (`array-like`, the filtered wavelength scale), "reps_filtered"
		# (`array-like`, the filtered real parts of the dielectric functions),
		# "ieps_filtered" (`array-like`, the filtered imaginary parts of the
		# dielectric functions), and "reason" (`array-like`, containing a code
		# to track whether a point was collected for step reasons [1], for
		# threshold filter [2], or for being a peak point [3]).
		def scan_single_file(config, file_name):
		result = {
		'exit_code': -1,
		@@ -338,7 +395,8 @@ def scan_single_file(config, file_name):
		'ieps_orig': [],
		'wl_filtered': [],
		'reps_filtered': [],
		'ieps_filtered': []
		'ieps_filtered': [],
		'reason': []
		}
		try:
		input_file = open(file_name, 'r')
		@@ -364,6 +422,7 @@ def scan_single_file(config, file_name):
		wl_filtered = result['wl_filtered']
		reps_filtered = result['reps_filtered']
		ieps_filtered = result['ieps_filtered']
		reason = result['reason']
		step = config['step']
		threshold = config['threshold']
		wl0 = 0.0
		@@ -396,6 +455,7 @@ def scan_single_file(config, file_name):
		wl_filtered.append(wl0 * wl_factor)
		reps_filtered.append(reps0)
		ieps_filtered.append(ieps0)
		reason.append(3)
		num_filtered_data += 1
		else:
		can_write = True
		@@ -433,6 +493,7 @@ def scan_single_file(config, file_name):
		wl_filtered.append(wl * wl_factor)
		reps_filtered.append(reps)
		ieps_filtered.append(ieps)
		reason.append(1)
		num_filtered_data += 1
		can_write = False
		wl0 = wl
		@@ -448,6 +509,7 @@ def scan_single_file(config, file_name):
		wl_filtered.append(wl1 * wl_factor)
		reps_filtered.append(reps1)
		ieps_filtered.append(ieps1)
		reason.append(3)
		num_filtered_data += 1
		wl0 = wl1
		reps0 = reps1
		@@ -465,6 +527,7 @@ def scan_single_file(config, file_name):
		wl_filtered.append(wl1 * wl_factor)
		reps_filtered.append(reps1)
		ieps_filtered.append(ieps1)
		reason.append(2)
		num_filtered_data += 1
		wl0 = wl1
		reps0 = reps1
		@@ -480,6 +543,7 @@ def scan_single_file(config, file_name):
		wl_filtered.append(wl1 * wl_factor)
		reps_filtered.append(reps1)
		ieps_filtered.append(ieps1)
		reason.append(2)
		num_filtered_data += 1
		wl0 = wl1
		reps0 = reps1
		@@ -512,29 +576,83 @@ def scan_single_file(config, file_name):
		# simulation.
		#
		# \param[in] config: `dict` A dictionary containing the script configuration.
		# \return result: `dict` A dictionary containing the results of the scan,
		# including "exit_code" (`int`, 0 if succesful), "wl_orig" (`array-like`, the
		# original wavelength scale), "reps_orig" (`array-like`, the original real
		# parts of the dielectric functions), "ieps_orig" (`array-like`, the original
		# imaginary parts of the dielectric functions), "wl_filtered" (`array-like`,
		# the filtered wavelength scale), "reps_filtered" (`array-like`, the filtered
		# real parts of the dielectric functions), and "ieps_filtered" (`array-like`,
		# the filtered imaginary parts of the dielectric functions).
		def scan_multiple_files(config):
		result = {
		'exit_code': -1,
		'header': "",
		'num_read_lines': 0,
		'num_orig_data': 0,
		'num_filtered_data': 0,
		'wl_orig': [],
		'reps_orig': [],
		'ieps_orig': [],
		'wl_filtered': [],
		'reps_filtered': [],
		'ieps_filtered': []
		}
		return result
		# \return aligned_infos: `list` A list of dictionaries containing the results
		# of filtering aligned to a common scale.
		def match_multiple_files(config):
		scan_infos = []
		split_input = config['input_files'].split(',')
		wl_factor = 1.0e6
		wl_units = config['wl_units']
		if (wl_units in ["nanometers", "nm"]):
		wl_factor = 1.0e9
		elif (wl_units in ["millimeters", "mm"]):
		wl_factor = 1.0e3
		elif (wl_units in ["centimeters", "cm"]):
		wl_factor = 1.0e2
		elif (wl_units in ["decimeters", "dm"]):
		wl_factor = 1.0e1
		elif (wl_units in ["meters", "m"]):
		wl_factor = 1.0
		for file_name in split_input:
		scan_infos.append(scan_single_file(config, file_name))
		# end scan_infos loop
		# Find the global X range
		x_min = min(np.min(s['wl_filtered']) for s in scan_infos)
		x_max = max(np.max(s['wl_filtered']) for s in scan_infos)
		tolerance = config['wl_tolerance'] if config['wl_tolerance'] != 0.0 else (1.0e-3 * x_min / wl_factor)
		# Extract special points
		special_ieps = []
		for s in scan_infos:
		x_arr = np.array(s['wl_filtered'])
		info_arr = np.array(s['reason'])
		special_ieps.extend(x_arr[info_arr > 1])
		# Make a regular grid
		#breakpoint()
		num_regular_points = int((x_max - x_min) / (config['step'] * wl_factor))
		regular_x = np.linspace(x_min, x_max, num_regular_points)
		# Get a coarse global X vector
		coarse_x = np.sort(np.unique(np.concatenate([regular_x, special_ieps])))
		# Numerical tolerance filtering
		mask = np.insert(np.diff(coarse_x) > tolerance, 0, True)
		common_x = coarse_x[mask]
		# Re-align each series on the common scale
		aligned_infos = []
		for s in scan_infos:
		wl_old = np.array(s['wl_filtered'])
		reps_old = np.array(s['reps_filtered'])
		ieps_old = np.array(s['ieps_filtered'])
		info_old = np.array(s['reason'])

		# Value interpolation on the new grid
		reps_interp = np.interp(common_x, wl_old, reps_old)
		ieps_interp = np.interp(common_x, wl_old, ieps_old)

		# Mapping of original INFO on the new scale
		info_new = np.ones(len(common_x), dtype=int)
		for x_val, info_val in zip(wl_old, info_old):
		if info_val > 1:
		# Find corresponding index in common_x
		idx = np.argmin(np.abs(common_x - x_val))
		if np.abs(common_x[idx] - x_val) <= tolerance:
		info_new[idx] = info_val
		# end of x_val, info_val loop

		aligned_infos.append({
		'exit_code': 0,
		'header': s['header'],
		'num_read_lines': s['num_read_lines'],
		'num_orig_data': s['num_orig_data'],
		'num_filtered_data': s['num_filtered_data'],
		'wl_orig': np.array(s['wl_orig']),
		'reps_orig': np.array(s['reps_orig']),
		'ieps_orig': np.array(s['ieps_orig']),
		'wl_filtered': common_x,
		'reps_filtered': reps_interp,
		'ieps_filtered': ieps_interp,
		'reason': info_new
		})
		# end of scan_infos loop
		return aligned_infos

		## \brief Write the filtered data to an output file.
		#