{0:0%dd}"%num_len
@@ -428,9 +429,15 @@ def mismatch_severities(str_f_values, str_c_values, config):
fractional = scale * (f_values[i] - c_values[i]) / f_values[i]
if (fractional < 0.0): fractional *= -1.0
if (fractional <= config['warning_threshold']):
- result[i] = 2
+ if (log_f_value > config['data_order']):
+ result[i] = 2
+ else:
+ result[i] = 1
else:
- result[i] = 3
+ if (log_f_value > config['data_order']):
+ result[i] = 3
+ else:
+ result[i] = 1
else:
result[i] = 1
else: # f_values[i] == 0 and c_values[i] != 0
@@ -441,9 +448,15 @@ def mismatch_severities(str_f_values, str_c_values, config):
fractional = scale * (c_values[i] - f_values[i]) / c_values[i]
if (fractional < 0.0): fractional *= -1.0
if (fractional <= config['warning_threshold']):
- result[i] = 2
+ if (log_c_value > config['data_order']):
+ result[i] = 2
+ else:
+ result[i] = 1
else:
- result[i] = 3
+ if (log_c_value > config['data_order']):
+ result[i] = 3
+ else:
+ result[i] = 1
else:
result[i] = 1
# End number comparison
@@ -460,31 +473,42 @@ def mismatch_severities(str_f_values, str_c_values, config):
# \returns config: `dict` A dictionary containing the script configuration.
def parse_arguments():
config = {
- 'fortran_file_name': '',
'c_file_name': '',
+ 'check_all': True,
+ 'data_order': -99.0,
+ 'fortran_file_name': '',
'full_log': False,
+ 'help_mode': False,
+ 'hide_noise': True,
+ 'html_output': 'pycompare.html',
'linewise': True,
'log_html': False,
- 'html_output': 'pycompare.html',
'say_progress': True,
'warning_threshold': 0.005,
- 'help_mode': False,
- 'check_all': True,
}
+ arg_index = 1
+ skip_arg = False
for arg in argv[1:]:
+ if skip_arg:
+ skip_arg = False
+ continue
split_arg = arg.split("=")
if (arg.startswith("--ffile")):
- config['fortran_file_name'] = split_arg[1]
+ config['fortran_file_name'] = argv[arg_index + 1]
+ arg_index += 1
+ skip_arg = True
elif (arg.startswith("--cfile")):
- config['c_file_name'] = split_arg[1]
+ config['c_file_name'] = argv[arg_index + 1]
+ arg_index += 1
+ skip_arg = True
+ elif (arg.startswith("--data-order")):
+ config['data_order'] = float(split_arg[1])
elif (arg.startswith("--full")):
config['full_log'] = True
elif (arg.startswith("--html")):
config['log_html'] = True
if (len(split_arg) == 2):
config['html_output'] = split_arg[1]
- elif (arg.startswith("--warn")):
- config['warning_threshold'] = float(split_arg[1])
elif (arg.startswith("--help")):
config['help_mode'] = True
elif (arg.startswith("--linewise")):
@@ -493,8 +517,13 @@ def parse_arguments():
config['say_progress'] = False
elif (arg.startswith("--quick")):
config['check_all'] = False
+ elif (arg.startswith("--show-noise")):
+ config['hide_noise'] = False
+ elif (arg.startswith("--warn")):
+ config['warning_threshold'] = float(split_arg[1])
else:
raise Exception("Unrecognized argument \'{0:s}\'".format(arg))
+ arg_index += 1
return config
## \brief Print a command-line help summary.
@@ -507,14 +536,16 @@ def print_help():
print("Usage: \"./pycompare.py OPTIONS\" ")
print(" ")
print("Valid options are: ")
- print("--ffile=FORTRAN_OUTPUT File containing the output of the FORTRAN code (mandatory).")
- print("--cfile=C++_OUTPUT File containing the output of the C++ code (mandatory).")
+ print("--ffile FORTRAN_OUTPUT File containing the output of the FORTRAN code (mandatory).")
+ print("--cfile C++_OUTPUT File containing the output of the C++ code (mandatory).")
+ print("--data-order=ORDER Consider data only down to specified order (default order is -99).")
print("--full Print all lines to log file (default prints only mismatches).")
print("--help Print this help and exit.")
print("--html[=OPT_OUTPUT_NAME] Enable logging to HTML file (default logs to \"pycompare.html\").")
print("--linewise Load only one line at a time. Useful to compare big files (True by default).")
print("--no-progress Disable progress logging.")
print("--quick Stop on first mismatch (default is to perform a full check).")
+ print("--show-noise Show noise in reports (default is to hide it).")
print("--warn Set a fractional threshold for numeric warning (default = 0.005).")
print(" ")
diff --git a/src/scripts/pydynrange.py b/src/scripts/pydynrange.py
new file mode 100755
index 0000000000000000000000000000000000000000..b9c74ca8bb5faf7f8c5f7ce3a29495ee7fc63ca4
--- /dev/null
+++ b/src/scripts/pydynrange.py
@@ -0,0 +1,243 @@
+#!/bin/python3
+
+# Copyright (C) 2024 INAF - Osservatorio Astronomico di Cagliari
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# A copy of the GNU General Public License is distributed along with
+# this program in the COPYING file. If not, see: .
+
+## @package pydynrange
+# \brief Script to calculate the dynamic range of a complex matrix.
+#
+# The script execution requires python3.
+
+import cmath
+import math
+import os
+
+from sys import argv, stdout
+
+## \brief Class to represent the dynamic range memory structure.
+#
+# In order to implement a light-weight test that is able to choose
+# whether to produce a text-only output or a text + diagnostic plot
+# output, the `PlotData` class stores the memory structures to
+# produce the text log and the plots. The data for the plots are
+# filled with dummy place-holders if the required output is text
+# only, while they are filled with actual data, otherwise.
+#
+# \returns exit_code: `int` 0 on successful completion.
+class PlotData:
+ # \brief PlotData instance constructor.
+ def __init__(self):
+ self.max_real = -1.0
+ self.max_imag = -1.0
+ self.max_absl = -1.0
+ self.min_real = 1.0
+ self.min_imag = 1.0
+ self.sml_absl = 1.0
+ self.sml_real = 1.0
+ self.sml_imag = 1.0
+ self.absl_hist = [0 for i in range(-99, 100)]
+ self.real_hist = [0 for i in range(-99, 100)]
+ self.imag_hist = [0 for i in range(-99, 100)]
+ self.max_real_mag = -100
+ self.min_real_mag = 100
+ self.max_imag_mag = -100
+ self.min_imag_mag = 100
+ self.max_absl_mag = -100
+ self.min_absl_mag = 100
+
+ # \brief Print a text log of the dynamic range.
+ def log_dynamic_range(self):
+ print(" MAX( ABS[AM] ) = %14.7e"%self.max_absl)
+ print(" MIN( ABS[AM] ) = %14.7e"%self.sml_absl)
+ print(" MAX( REAL[AM] ) = %14.7e"%self.max_real)
+ print(" MIN( REAL[AM] ) = %14.7e"%self.min_real)
+ print(" MAX( IMAG[AM] ) = %14.7e"%self.max_imag)
+ print(" MIN( IMAG[AM] ) = %14.7e"%self.min_imag)
+ print("MIN( ABS( REAL[AM] ) ) = %14.7e"%self.sml_real)
+ print("MIN( ABS( IMAG[AM] ) ) = %14.7e"%self.sml_imag)
+ return
+
+ # \brief Make histograms of dynamic range with matplotlib.
+ def plot_dynamic_range(self, config):
+ import matplotlib.pyplot as plt
+ import numpy as np
+ self.absl_hist = np.array(self.absl_hist)
+ self.real_hist = np.array(self.real_hist)
+ self.imag_hist = np.array(self.imag_hist)
+ vec_x_pos = np.arange(-98.5, 100.0, 1.0)
+ fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, ncols=1)
+ ax1.bar(vec_x_pos, self.absl_hist, color = "grey", label = "Abs. values")
+ ax1.set_xlim(self.min_absl_mag - 1, self.max_absl_mag + 1)
+ ax1.set_xlabel("Order of magnitude")
+ ax1.set_ylabel("Element counts")
+ ax1.legend(loc = "best")
+ ax2.bar(vec_x_pos, self.real_hist, color = "red", label = "Real values")
+ ax2.set_xlim(self.min_real_mag - 1, self.max_real_mag + 1)
+ ax2.set_xlabel("Order of magnitude")
+ ax2.set_ylabel("Element counts")
+ ax2.legend(loc = "best")
+ ax3.bar(vec_x_pos, self.imag_hist, color = "blue", label = "Imag. values")
+ ax3.set_xlim(self.min_imag_mag - 1, self.max_imag_mag + 1)
+ ax3.set_xlabel("Order of magnitude")
+ ax3.set_ylabel("Element counts")
+ ax3.legend(loc = "best")
+ if (config['log_plot_y']):
+ ax1.set_yscale("log")
+ ax2.set_yscale("log")
+ ax3.set_yscale("log")
+ plt.tight_layout()
+ plt.show()
+ return
+
+## \brief Main execution code
+#
+# `main()` is the function that handles the creation of the script configuration
+# and the execution of the calculation. It returns 0 on successful completion.
+#
+# \returns exit_code: `int` 0 on successful completion.
+def main():
+ config = parse_arguments()
+ exit_code = 0
+ if config['help_mode'] or len(argv) == 1:
+ config['help_mode'] = True
+ print_help()
+ else:
+ if config['matrix_name'] == "":
+ exit_code = 1
+ else:
+ exit_code = get_dynamic_range(config)
+ return exit_code
+
+## \brief Compute the dynamic range of a matrix.
+#
+# \param config: `dict` A dictionary containing the script configuration.
+#
+# \returns result: `int` The exit code of the operation (0 for success).
+def get_dynamic_range(config):
+ result = 0
+ num_read_lines = 0
+ pd = PlotData()
+ stdout.write("INFO: scanning data file. Please, wait...")
+ stdout.flush()
+ matrix_file = open(config['matrix_name'], 'r')
+ str_line = matrix_file.readline()
+ str_line = matrix_file.readline()
+ str_line = matrix_file.readline()
+ while (str_line != ""):
+ str_value = str_line.split("(")[1][:-2]
+ split_value = str_value.split(",")
+ real = float(split_value[0])
+ imag = float(split_value[1])
+ if (real > pd.max_real): pd.max_real = real
+ # if (config['make_plots'] and real != 0.0):
+ # mag_value = int(math.log10(abs(real)).floor())
+ # if (mag_value > pd.max_real_mag): pd.max_real_mag = mag_value
+ # if (mag_value < pd.min_real_mag): pd.min_real_mag = mag_value
+ # pd.real_hist[mag_value + 99] += 1
+ if (imag > pd.max_imag): pd.max_imag = imag
+ if (real < pd.min_real): pd.min_real = real
+ if (imag < pd.min_imag): pd.min_imag = imag
+ cvalue = complex(real, imag)
+ absl_val = abs(cvalue)
+ if (absl_val > pd.max_absl): pd.max_absl = absl_val
+ if (absl_val < pd.sml_absl): pd.sml_absl = absl_val if absl_val > 0.0 else pd.sml_absl
+ if (real < 0): real *= -1.0
+ if (imag < 0): imag *= -1.0
+ if (real < pd.sml_real): pd.sml_real = real if real > 0.0 else pd.sml_real
+ if (imag < pd.sml_imag): pd.sml_imag = imag if imag > 0.0 else pd.sml_real
+ if (config['make_plots']):
+ if (real > 0.0):
+ mag_value = int(math.floor(math.log10(real)))
+ if (mag_value > pd.max_real_mag): pd.max_real_mag = mag_value
+ if (mag_value < pd.min_real_mag): pd.min_real_mag = mag_value
+ pd.real_hist[mag_value + 99] += 1
+ if (imag > 0.0):
+ mag_value = int(math.floor(math.log10(imag)))
+ if (mag_value > pd.max_imag_mag): pd.max_imag_mag = mag_value
+ if (mag_value < pd.min_imag_mag): pd.min_imag_mag = mag_value
+ pd.imag_hist[mag_value + 99] += 1
+ if (absl_val > 0.0):
+ mag_value = int(math.floor(math.log10(absl_val)))
+ if (mag_value > pd.max_absl_mag): pd.max_absl_mag = mag_value
+ if (mag_value < pd.min_absl_mag): pd.min_absl_mag = mag_value
+ pd.absl_hist[mag_value + 99] += 1
+ str_line = matrix_file.readline()
+ if (config['limit'] > 0):
+ num_read_lines += 1
+ if (num_read_lines >= config['limit']):
+ str_line = "" # Close the while loop
+ matrix_file.close()
+ print(" done!")
+ pd.log_dynamic_range()
+ if (config['make_plots']):
+ pd.plot_dynamic_range(config)
+ return result
+
+## \brief Parse the command line arguments.
+#
+# The script behaviour can be modified through a set of mandatory and optional
+# arguments. The only mandatory argument is the name of the log file to be
+# parsed. Additional optional arguments are an operation filter, which should
+# be the starting sequence of the log strings to pe included in the timing
+# calculation and the number of threads used during code execution.
+#
+# \returns config: `dict` A dictionary containing the script configuration.
+def parse_arguments():
+ config = {
+ 'help_mode': False,
+ 'limit': -1,
+ 'log_plot_y': True,
+ 'make_plots': False,
+ 'matrix_name': "",
+ }
+ for arg in argv[1:]:
+ if (arg.startswith("--help")):
+ config['help_mode'] = True
+ elif (arg.startswith("--limit=")):
+ split_arg = arg.split('=')
+ config['limit'] = int(split_arg[1])
+ elif (arg.startswith("--lin-y")):
+ config['log_plot_y'] = False
+ elif (arg.startswith("--make-plots")):
+ config['make_plots'] = True
+ else:
+ if (os.path.isfile(arg)):
+ config['matrix_name'] = arg
+ else:
+ raise Exception("Unrecognized argument \'{0:s}\'".format(arg))
+ return config
+
+## \brief Print a command-line help summary.
+def print_help():
+ print(" ")
+ print("**************************** PYDYNRANGE ****************************")
+ print(" ")
+ print("Get the dynamic range of a complex matrix. ")
+ print(" ")
+ print("Usage: \"./pydynrange.py FILE_NAME [OPTIONS]\" ")
+ print(" ")
+ print("Valid options are: ")
+ print("--help Print this help and exit.")
+ print("--limit=NUMBER Check only NUMBER file lines.")
+ print("--lin-y Use linear scale on plot y-axes.")
+ print("--make-plots Plot histograms of magnitudes (requires matplotlib).")
+ print(" ")
+
+
+# ### PROGRAM EXECUTION ###
+## \cond
+res = main()
+## \endcond
+exit(res)
diff --git a/src/trapping/cfrfme.cpp b/src/trapping/cfrfme.cpp
index c00772293c5d9f9832168237329e8138b615a6fb..d2c19f3cec6abf5358e09d9f91975d0739019bd4 100644
--- a/src/trapping/cfrfme.cpp
+++ b/src/trapping/cfrfme.cpp
@@ -242,7 +242,7 @@ void frfme(string data_file, string output_path) {
// Array initialization
long swap1_size, swap2_size, tfrfme_size;
double size_mb;
- printf("INFO: calculation memory requirements\n");
+ printf("INFO: calculating memory requirements\n");
swap1_size = Swap1::get_memory_requirement(lm, nkv);
size_mb = 1.0 * swap1_size / 1024.0 / 1024.0;
printf("Swap 1: %.2lg MB\n", size_mb);
diff --git a/src/trapping/clffft.cpp b/src/trapping/clffft.cpp
index b36383756520522b7ff8f7e10f7a65626b677150..46828b2ddfb1d5de61249f6a570fa704dfc1f9e5 100644
--- a/src/trapping/clffft.cpp
+++ b/src/trapping/clffft.cpp
@@ -326,15 +326,18 @@ void lffft(string data_file, string output_path) {
if (is != 2222) {
if (is != 1111) {
if (is > 0) { // Goes to 305
+ if (ac) delete[] ac;
ac = new dcomplex[cil->nlemt]();
camp(ac, am0m, ws, cil);
// Goes to 445
} else if (is < 0) { // Goes to 405
+ if (ac) delete[] ac;
ac = new dcomplex[cil->nlemt]();
czamp(ac, amd, indam, ws, cil);
// Goes to 445
}
} else {
+ if (ac) delete[] ac;
ac = new dcomplex[cil->nlemt]();
samp(ac, tmsm, tmse, ws, cil);
// Goes to 445