Commit 35ce77dc authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Implement single file filtering logic and inline documentation

parent d7f9bd25
Loading
Loading
Loading
Loading
+93 −32
Original line number Diff line number Diff line
@@ -27,7 +27,7 @@
#  The script requires python3.

import math
import pdb
# import pdb
from sys import argv

## \cond
@@ -66,8 +66,27 @@ def main():
                result = scan_multiple_files(config)
    return result

## \brief Maximum distance between input data and linear interpolation of filtered data.
#  
#  This function returns a dictionary object to make a quick estimate of the quality
#  of the chosen filter. The diagnostic is the maximum absolute offset between the
#  distribution of the input data and a segmented linear interpolation touching all
#  the filtered data.
#
#  \param[in] wl_orig: `array-like` Array of the input data wavelengths.
#  \param[in] reps_orig: `array-like` Array of the real part of the input optical functions.
#  \param[in] ieps_orig: `array-like` Array of the imaginary part of the input optical functions.
#  \param[in] wl_filtered: `array-like` Array of the filtered data wavelengths.
#  \param[in] reps_filtered: `array-like` Array of the real part of the filtered optical functions.
#  \param[in] ieps_filtered: `array-like` Array of the imaginary part of the filtered optical functions.
#  \returns result: `dict` A dictionary containing the wavelength of the maximum difference
#           (`max_wavelength`), the value of the input data at that wavelength (`mav_value`),
#           the value of the interpolated filtered functions at the same wavelength (`max_fitted`),
#           the offset between interpolation and data (`max_difference`), and the set of values
#           where the difference was observed (`differing_set`, being either REAL or IMAGINARY).
def get_max_distance(wl_orig, reps_orig, ieps_orig, wl_filtered, reps_filtered, ieps_filtered):
    max_difference = 0.0
    max_fitted = 0.0
    max_value = 0.0
    max_wavelength = 0.0
    differing_set = ""
@@ -91,11 +110,13 @@ def get_max_distance(wl_orig, reps_orig, ieps_orig, wl_filtered, reps_filtered,
                idiff *= -1.0
            if (rdiff > max_difference):
                max_difference = rdiff
                max_fitted = rp
                max_value = reps_orig[i]
                max_wavelength = wl
                differing_set = "real"
            if (idiff > max_difference):
                max_difference = idiff
                max_fitted = ip
                max_value = ieps_orig[i]
                max_wavelength = wl
                differing_set = "imaginary"
@@ -107,6 +128,7 @@ def get_max_distance(wl_orig, reps_orig, ieps_orig, wl_filtered, reps_filtered,
    result = {
        'max_wavelength': max_wavelength,
        'max_difference': max_difference,
        'max_fitted': max_fitted,
        'max_value': max_value,
        'differing_set': differing_set
    }
@@ -179,6 +201,7 @@ def parse_arguments():
                raise Exception("Unrecognized wavelength units %s!"%config['wl_units'])
        else:
            raise Exception("Unrecognized argument \"{0:s}\"!".format(arg))
    # end for loop
    return config

## \brief Print a command-line help summary.
@@ -194,14 +217,37 @@ def print_help():
    print("Usage: \"./filter_constants.py --in INPUT [options]\"                 ")
    print("                                                                      ")
    print("Valid options are:                                                    ")
    print("--in                  Comma separated list of input files.")
    print("                      (mandatory).                        ")
    print("--in                  Comma separated list of input files (mandatory).")
    print("--out                 Comma separated list of output files (optional, ")
    print("                      but, if given, must be as many as for --in).    ")
    print("--help                Print this help and exit.                       ")
    print("--no-peaks            Disable mandatory extraction of peaks (default  ")
    print("                      is enabling peaks).                             ")
    print("--no-plots            Disable plotting the solution with MATPLOTLIB   ")
    print("                      (default is enable plots).                      ")
    print("--step=VALUE          Regular step in meters to sample flat regions   ")
    print("                      of the functions (use <= 0 to disable; default  ")
    print("                      is 5e-8).                                       ")
    print("--threshold=VALUE     Relative tolerance threshold to pick a point in ")
    print("                      the filtering process (optional, default is     ")
    print("                      0.1).                                           ")
    print("--wl_start=VALUE      Starting wavelength in meters for the filtering ")
    print("                      window.                                         ")
    print("--wl_stop=VALUE       Ending wavelength in meters for the filtering   ")
    print("                      window.                                         ")
    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).   ")
    print("--version             Print script version and exit.                  ")
    print("                                                                      ")

## \brief Filter a single file based on the configuration options.
#
#  Perform filtering of a single file based on custom thresholda and step
#  configurations. The filtered data are written to a CSV file, then a
#  diagnostic log is printed to terminal. Optionally the filter selection
#  is shown as a plot, if MATPLOTLIB is available on the system.
#
#  \param[in] config: `dict` A dictionary containing the script configuration.
#  \return result: `int` An integer exit code (0 if successful).
def scan_single_file(config):
@@ -216,9 +262,10 @@ def scan_single_file(config):
        output_name = config['output_files']
        if (output_name == ''):
            output_name = file_name.split('.')[0] + "_filtered.csv"
        split_output = output_name.split('/')
        if (len(split_output) > 1):
            output_name = split_output[-1]
        if (output_name == config['input_files']):
            print("ERROR: overwriting of input files is not supported!")
            result = 2
            return result
        output_file = open(output_name, 'w')
        wl_factor = 1.0e6
        wl_units = config['wl_units']
@@ -267,7 +314,6 @@ def scan_single_file(config):
                    ieps_orig.append(ieps0)
                    num_orig_data += 1
                    if (wl0 >= config['wl_start']):
                        breakpoint()
                        output_file.write(file_line)
                        wl_filtered.append(wl0 * wl_factor)
                        reps_filtered.append(reps0)
@@ -292,17 +338,29 @@ def scan_single_file(config):
                        break # while loop
                    if (step > 0.0):
                        if (wl1 - wl0 >= step):
                            #breakpoint()
                            # compute the values at step location with linear interpolation
                            wl = wl0 + step
                            x0 = wl_orig[-2] / wl_factor if len(wl_orig) > 1 else wl0
                            x1 = wl_orig[-1] / wl_factor if len(wl_orig) > 1 else wl1
                            dx = wl - x0
                            ry0 = reps_orig[-2] if len(reps_orig) > 1 else reps0
                            ry1 = reps_orig[-1] if len(reps_orig) > 1 else reps1
                            dry = ry1 - ry0
                            iy0 = ieps_orig[-2] if len(ieps_orig) > 1 else ieps0
                            iy1 = ieps_orig[-1] if len(ieps_orig) > 1 else ieps1
                            diy = iy1 - iy0
                            reps = ry0 + dry * dx / (x1 - x0)
                            ieps = iy0 + diy * dx / (x1 - x0)
                            # write a line if step is enabled and satisfied
                            output_file.write(file_line)
                            wl_filtered.append(wl1 * wl_factor)
                            reps_filtered.append(reps1)
                            ieps_filtered.append(ieps1)
                            wl_filtered.append(wl * wl_factor)
                            reps_filtered.append(reps)
                            ieps_filtered.append(ieps)
                            num_filtered_data += 1
                            can_write = False
                            wl0 = wl1
                            reps0 = reps1
                            ieps0 = ieps1
                            wl0 = wl
                            reps0 = reps
                            ieps0 = ieps
                    # end of step > 0.0 check
                    if (config['force_peaks']):
                        rpeak = (dreps * last_dreps < 0.0)
@@ -374,7 +432,10 @@ def scan_single_file(config):
            config['wl_units'],
            max_distance['differing_set']
        )
        msg_distance2 = "      (actual data value is {0:.5g}).".format(max_distance['max_value'])
        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)