Commit e939fed7 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Improve parsing scripts

parent 59f29e3d
Loading
Loading
Loading
Loading
+18 −4
Original line number Diff line number Diff line
@@ -101,10 +101,12 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
    rmi[ri] = vec_rmi + (_nsph * ri);
    rei[ri] = vec_rei + (_nsph * ri);
  }
  int nllt = (_nlemt == 0) ? 2 * _nsph * _li * (_li + 2) : _nlemt;
  vec_w = new dcomplex[nllt * 4]();
  w = new dcomplex*[nllt];
  for (int wi = 0; wi < nllt; wi++) w[wi] = vec_w + (4 * wi);
  // int nllt = (_nlemt == 0) ? 2 * _nsph * _li * (_li + 2) : _nlemt;
  // vec_w = new dcomplex[nllt * 4]();
  // w = new dcomplex*[nllt];
  // for (int wi = 0; wi < nllt; wi++) w[wi] = vec_w + (4 * wi);
  vec_w = NULL;
  w = NULL;
  vec_rc = new double[num_layers]();
  rc = new double*[num_configurations];
  int last_layer_index = 0, cur_layer_index = 0;
@@ -688,6 +690,10 @@ ParticleDescriptorCluster::ParticleDescriptorCluster(GeometryConfiguration *gcon
  _nv3j = (_lm * (_lm + 1) * (2 * _lm + 7)) / 6;
  _ndi = _nsph * _nlim;
  _ndit = 2 * _nsph * _nlim;
  int nllt = (_nlemt == 0) ? 2 * _nsph * _li * (_li + 2) : _nlemt;
  vec_w = new dcomplex[nllt * 4]();
  w = new dcomplex*[nllt];
  for (int wi = 0; wi < nllt; wi++) w[wi] = vec_w + (4 * wi);
  vec_am0m = new dcomplex[_nlemt * _nlemt]();
  vec_fsac = new dcomplex[4]();
  vec_sac = new dcomplex[4]();
@@ -1086,6 +1092,10 @@ ParticleDescriptorInclusion::ParticleDescriptorInclusion(GeometryConfiguration *
  _nv3j = (_lm * (_lm + 1) * (2 * _lm + 7)) / 6;
  _ndi = _nsph * _nlim;
  _ndit = 2 * _nsph * _nlim;
  int nllt = (_nlemt == 0) ? 2 * _nsph * _li * (_li + 2) : _nlemt;
  vec_w = new dcomplex[nllt * 4]();
  w = new dcomplex*[nllt];
  for (int wi = 0; wi < nllt; wi++) w[wi] = vec_w + (4 * wi);
  vec_am0m = new dcomplex[_nlemt * _nlemt]();
  vec_fsac = new dcomplex[4]();
  vec_sac = new dcomplex[4]();
@@ -1399,6 +1409,10 @@ ParticleDescriptorSphere::ParticleDescriptorSphere(GeometryConfiguration *gconf,
  vec_sas = new dcomplex[4 * _nsph]();
  vec_vints = new dcomplex[16 * _nsph]();

  int nllt = (_nlemt == 0) ? 2 * _nsph * _li * (_li + 2) : _nlemt;
  vec_w = new dcomplex[nllt * 4]();
  w = new dcomplex*[nllt];
  for (int wi = 0; wi < nllt; wi++) w[wi] = vec_w + (4 * wi);
  fsas = new dcomplex[_nsph];
  sscs = new double[_nsph]();
  sexs = new double[_nsph]();
+40 −40
Original line number Diff line number Diff line
@@ -1867,12 +1867,12 @@ int ClusterOutputInfo::write_legacy(const std::string &output) {
		real(vec_fsas[jxi * configurations + si]),
		imag(vec_fsas[jxi * configurations + si])
        );
	fprintf(
		p_outfile, "INSERTION: CS_SPHERE  %15.7lE%15.7lE%15.7lE%15.7lE\n",
		alamb, vec_sphere_scs[jxi * configurations + si],
		vec_sphere_abs[jxi * configurations + si],
		vec_sphere_exs[jxi * configurations + si]
        );
	// fprintf(
	// 	p_outfile, "INSERTION: CS_SPHERE  %15.7lE%15.7lE%15.7lE%15.7lE\n",
	// 	alamb, vec_sphere_scs[jxi * configurations + si],
	// 	vec_sphere_abs[jxi * configurations + si],
	// 	vec_sphere_exs[jxi * configurations + si]
        // );
	fprintf(
		p_outfile, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
		vec_qschus[jxi * configurations + si],
@@ -1943,10 +1943,10 @@ int ClusterOutputInfo::write_legacy(const std::string &output) {
	      vec_cosavc1[jxi], vec_raprc1[jxi]
      );
      fprintf(p_outfile, "  Fk=%15.7lE\n", vec_fkc1[jxi]);
      fprintf(
	      p_outfile, "INSERTION: CSM_CLUSTER  %15.7lE%15.7lE%15.7lE%15.7lE\n",
	      alamb, vec_scc1[jxi], vec_abc1[jxi], vec_exc1[jxi]
      );
      // fprintf(
      // 	      p_outfile, "INSERTION: CSM_CLUSTER  %15.7lE%15.7lE%15.7lE%15.7lE\n",
      // 	      alamb, vec_scc1[jxi], vec_abc1[jxi], vec_exc1[jxi]
      // );
      // Perpendicular polarization cluster average section
      if (inpol == 0)
	fprintf(p_outfile, "   LIN  1\n");
@@ -2148,13 +2148,13 @@ int ClusterOutputInfo::write_legacy(const std::string &output) {
		      vec_dir_pschuc1[sat_dir_index],
		      vec_dir_s0magc1[sat_dir_index]
	      );
	      fprintf(
		      p_outfile, "INSERTION: CS1_CLUSTER  %13.5le%10.3le%10.3le%15.7le%15.7le%15.7le\n",
		      alamb, th + jth * thstp, ths + jths * thsstp,
		      vec_dir_scc1[sat_dir_index],
		      vec_dir_abc1[sat_dir_index],
		      vec_dir_exc1[sat_dir_index]
	      );
	      // fprintf(
	      // 	      p_outfile, "INSERTION: CS1_CLUSTER  %13.5le%10.3le%10.3le%15.7le%15.7le%15.7le\n",
	      // 	      alamb, th + jth * thstp, ths + jths * thsstp,
	      // 	      vec_dir_scc1[sat_dir_index],
	      // 	      vec_dir_abc1[sat_dir_index],
	      // 	      vec_dir_exc1[sat_dir_index]
	      // );
	      if (!goto190) {
		fprintf(
			p_outfile, "  COSAV=%15.7lE, RAPRS=%15.7lE\n",
@@ -2249,13 +2249,13 @@ int ClusterOutputInfo::write_legacy(const std::string &output) {
		      vec_dir_pschuc2[sat_dir_index],
		      vec_dir_s0magc2[sat_dir_index]
	      );
	      fprintf(
		      p_outfile, "INSERTION: CS2_CLUSTER  %13.5le%10.3le%10.3le%15.7le%15.7le%15.7le\n",
		      alamb, th + jth * thstp, ths + jths * thsstp,
		      vec_dir_scc2[sat_dir_index],
		      vec_dir_abc2[sat_dir_index],
		      vec_dir_exc2[sat_dir_index]
	      );
	      // fprintf(
	      // 	      p_outfile, "INSERTION: CS2_CLUSTER  %13.5le%10.3le%10.3le%15.7le%15.7le%15.7le\n",
	      // 	      alamb, th + jth * thstp, ths + jths * thsstp,
	      // 	      vec_dir_scc2[sat_dir_index],
	      // 	      vec_dir_abc2[sat_dir_index],
	      // 	      vec_dir_exc2[sat_dir_index]
	      // );
	      if (!goto190) {
		fprintf(
			p_outfile, "  COSAV=%15.7lE, RAPRS=%15.7lE\n",
@@ -4196,9 +4196,9 @@ int InclusionOutputInfo::write_legacy(const std::string &output) {
	  p_outfile, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
	  vec_scs1[jxi], vec_abs1[jxi], vec_exs1[jxi], vec_albeds1[jxi]
	);
	fprintf(p_outfile, "INSERTION: SCASECM %5d%15.7E%15.7E%15.7E%15.7E\n",
		-1, alamb, vec_scs1[jxi], vec_abs1[jxi], vec_exs1[jxi]
	);
	// fprintf(p_outfile, "INSERTION: SCASECM %5d%15.7E%15.7E%15.7E%15.7E\n",
	// 	-1, alamb, vec_scs1[jxi], vec_abs1[jxi], vec_exs1[jxi]
	// );
	fprintf(p_outfile, " ---- SCS/GS -- ABC/GS -- EXS/GS ---\n");
	fprintf(
	  p_outfile, " %14.7lE%15.7lE%15.7lE\n",
@@ -4225,9 +4225,9 @@ int InclusionOutputInfo::write_legacy(const std::string &output) {
	  p_outfile, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
	  vec_scs2[jxi], vec_abs2[jxi], vec_exs2[jxi], vec_albeds2[jxi]
	);
	fprintf(p_outfile, "INSERTION: SCASECM %5d%15.7E%15.7E%15.7E%15.7E\n",
		1, alamb, vec_scs2[jxi], vec_abs2[jxi], vec_exs2[jxi]
	);
	// fprintf(p_outfile, "INSERTION: SCASECM %5d%15.7E%15.7E%15.7E%15.7E\n",
	// 	1, alamb, vec_scs2[jxi], vec_abs2[jxi], vec_exs2[jxi]
	// );
	fprintf(p_outfile, " ---- SCS/GS -- ABC/GS -- EXS/GS ---\n");
	fprintf(
	  p_outfile, " %14.7lE%15.7lE%15.7lE\n",
@@ -4314,11 +4314,11 @@ int InclusionOutputInfo::write_legacy(const std::string &output) {
		  vec_dir_exs1[dir_index],
		  vec_dir_albeds1[dir_index]
		);
		fprintf(
		  p_outfile, "INSERTION: SCASEC %5d%14.7lE%14.7lE%14.7lE%14.7lE\n",
		  -1, alamb, vec_dir_scs1[dir_index], vec_dir_abs1[dir_index],
		  vec_dir_exs1[dir_index]
		);
		// fprintf(
		//   p_outfile, "INSERTION: SCASEC %5d%14.7lE%14.7lE%14.7lE%14.7lE\n",
		//   -1, alamb, vec_dir_scs1[dir_index], vec_dir_abs1[dir_index],
		//   vec_dir_exs1[dir_index]
		// );
		fprintf(p_outfile, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
		fprintf(
		  p_outfile, " %14.7lE%15.7lE%15.7lE\n",
@@ -4392,11 +4392,11 @@ int InclusionOutputInfo::write_legacy(const std::string &output) {
		  vec_dir_exs2[dir_index],
		  vec_dir_albeds2[dir_index]
		);
		fprintf(
		  p_outfile, "INSERTION: SCASEC %5d%14.7lE%14.7lE%14.7lE%14.7lE\n",
		  1, alamb, vec_dir_scs2[dir_index], vec_dir_abs2[dir_index],
		  vec_dir_exs2[dir_index]
		);
		// fprintf(
		//   p_outfile, "INSERTION: SCASEC %5d%14.7lE%14.7lE%14.7lE%14.7lE\n",
		//   1, alamb, vec_dir_scs2[dir_index], vec_dir_abs2[dir_index],
		//   vec_dir_exs2[dir_index]
		// );
		fprintf(p_outfile, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
		fprintf(
		  p_outfile, " %14.7lE%15.7lE%15.7lE\n",
+42 −2
Original line number Diff line number Diff line
@@ -6,7 +6,17 @@ This directory contains scripts and tools to evaluate the code functionality.

The code migration stage can be considered successfully fulfilled with the solution of the single sphere and of the cluster of spheres cases in *C++*. To test the reliability of the code, the *C++* version needs to produce consistent output with respect to the original *FORTRAN* code. Since the output files are generally too large for manual inspection and they can be affected by different approximations or numeric noise, a series of scripts designed to compare the two versions and to pin-point possible differences is provided in this folder. The scripts are written in *python*, therefore they need an available *python* environment to work.

## Comparing the results of *FORTRAN* and *C++* codes
## Create models from YAML input

The execution of `NP_TMcode` applications requires machine readable input files to describe the particle model and the radiation field properties. This files are ASCII files that can be edited with common text editors, but they obey a strict format which complies with the requirements of the original code implementation. In order to make the model creation easier, `NP_TMcode` provides a model editor script named `model_maker.py`. This script reads a model description from a `YAML` file and creates the necessary input configuration for the `NP_TMcode` applications. To run the script, just issue:

   > $PATH_TO_SCRIPT/model_maker.py YAML_CONFIGURATION_FILE

where `YAML_CONFIGURATION_FILE` must be a valid `YAML` model description (examples are given in the `test_data` folder). More details are given issuing:

   > $PATH_TO_SCRIPT/model_maker.py --help

## Compare the results of *FORTRAN* and *C++* codes

1. Follow the instructions to build and run the *FORTRAN* and *C++* versions of the code.
2. Run the `pycompare.py` script providing the following minimal arguments:
@@ -27,7 +37,17 @@ or just:

will print a help screen, giving a brief explanation of all the possible options.

## Estimating the execution time from a terminal log
## Estimate the recommended orders for running a model

This script applies the Wiscombe criterion and its modifications to estimate the recommended internal and external orders to run a model. To use it, issue:

   > $PATH_TO_SCRIPT/pywiscombe.py --li|le --wave=WAVELENGTH --rad=RADIUS

Wavelength and radius must be expressed in meters, with the radius being equal to the minimum radius encircling the particle. More details are given with:

   > $PATH_TO_SCRIPT/pywiscombe.py --help

## Estimate the execution time from a terminal log

Performance estimates can be obtained from the code logging system, assuming the uses chose to save the terminal output to some log file. To obtain the time spent by the code in performing a specific operation, the syntax is:

@@ -35,6 +55,26 @@ Performance estimates can be obtained from the code logging system, assuming the

where `LOG_FILE` must be the name of a file containing the output that would normally go to terminal, `FILTER` must be the starting character sequence of the log line identifying the operation that should be taken into account, and `NUM_THREADS` is the number of processes that were used to perform the whole calculation loop. In case no filter is given, the script provides an estimate of the total amount of time spent in doing the calculation. This estimate, however, is known not to be reliable, because it ignores thread concurrency effects. A more accurate estimate of the total time spent in executing the code is always saved in a file named `c_timing.log`.

## Extract parameters from the code output

The legacy and `HDF5` outputs of `NP_TMcode` applications contain all the values computed by the code. However, in many cases the user is interested in extracting plots of cross-sections, forces and torques as functions of wavelength. The `parse_output.py` script performs this task (currently only on the legacy output). To use the script after running a model, issue:

   > $PATH_TO_SCRIPT/parse_output.py --in CODE_RESULT_FILE --out OUTPUT_BASE_NAME --app=[SPH|CLU|INCLU]

This will parse the output of `np_sphere` (`app=SPH`), `np_cluster` (`app=CLU`) or `np_inclusion` (`app=INCLU`) and put the results in `CSV` text files, named after the given `OUTPUT_BASE_NAME`, but with distinguishing suffixes such as `_ics` for integrated cross-sections, `_dcs` for differential ones, `_irp` for integrated radiation pressure forces, etc. More details are given issuing:

   > $PATH_TO_SCRIPT/parse_output.py --help

## Extract the dynamic range of a complex matrix

This script examines a complex matrix and extracts histograms of the orders of magnitude of the real parts, the imaginary parts and the absolute values. Optionally, if `matplotlib` is available, a plot of the dynamic range is drawn. To use this script, issue:

   > $PATH_TO_SCRIPT/pydynrange.py MATRIX_FILE

where `MATRIX_FILE` must be a `CSV` representation of the matrix. More details are given issuing:

   > $PATH_TO_SCRIPT/pydynrange.py --help

# License

   Copyright (C) 2025   INAF - Osservatorio Astronomico di Cagliari
+6 −2
Original line number Diff line number Diff line
@@ -879,9 +879,13 @@ def write_legacy_gconf(conf):
    output.write(str_line)
    if (conf['use_refinement']):
        str_line = "USE_REFINEMENT=1\n"
    else:
        str_line = "USE_REFINEMENT=0\n"
    output.write(str_line)
    if (not conf['dyn_orders']):
        str_line = "USE_DYN_ORDER=0\n"
    else:
        str_line = "USE_DYN_ORDER=1\n"
    output.write(str_line)
    output.close()
    return result
+57 −7
Original line number Diff line number Diff line
@@ -149,9 +149,9 @@ def parse_legacy_oclu(config):
            out40.write("Wavelength,CosAv,RaPr\n")
        if ('ALL' in config['selection'] or 'DRP' in config['selection']):
            out51 = open(root_name + "_drp1.csv", "w") # open a file for differential radiation pressure forces in state -1
            out51.write("Wavelength,THi,THs,PHi,PHs,CosAv,RaPr\n")
            out52 = open(root_name + "_drp1.csv", "w") # open a file for differential radiation pressure forces in state +1
            out52.write("Wavelength,THi,THs,PHi,PHs,CosAv,RaPr\n")
            out51.write("Wavelength,THi,THs,PHi,PHs,CosAv,RaPr,Fl,Fr,Fk,Fx,Fy,Fz,TQEl,TQEr,TQEk,TQEx,TQEy,TQEz,TQSl,TQSr,TQSk,TQSx,TQSy,TQSz\n")
            out52 = open(root_name + "_drp2.csv", "w") # open a file for differential radiation pressure forces in state +1
            out52.write("Wavelength,THi,THs,PHi,PHs,CosAv,RaPr,Fl,Fr,Fk,Fx,Fy,Fz,TQEl,TQEr,TQEk,TQEx,TQEy,TQEz,TQSl,TQSr,TQSk,TQSx,TQSy,TQSz\n")

        # Define the quantities that you need to extract
        alam = 0.0
@@ -269,12 +269,37 @@ def parse_legacy_oclu(config):
                                    file_line = oclu_file.readline()
                            cosav = float(file_line[8:23].replace("D", "E"))
                            rapr = float(file_line[31:46].replace("D", "E"))
                            # we read the forces and torques
                            file_line = oclu_file.readline()
                            fl = float(file_line[5:20].replace("D", "E"))
                            fr = float(file_line[25:40].replace("D", "E"))
                            fk = float(file_line[45:60].replace("D", "E"))
                            file_line = oclu_file.readline()
                            fx = float(file_line[5:20].replace("D", "E"))
                            fy = float(file_line[25:40].replace("D", "E"))
                            fz = float(file_line[45:60].replace("D", "E"))
                            file_line = oclu_file.readline()
                            TQEl = float(file_line[8:23].replace("D", "E"))
                            TQEr = float(file_line[31:46].replace("D", "E"))
                            TQEk = float(file_line[54:69].replace("D", "E"))
                            file_line = oclu_file.readline()
                            TQSl = float(file_line[8:23].replace("D", "E"))
                            TQSr = float(file_line[31:46].replace("D", "E"))
                            TQSk = float(file_line[54:69].replace("D", "E"))
                            file_line = oclu_file.readline()
                            TQEx = float(file_line[8:23].replace("D", "E"))
                            TQEy = float(file_line[31:46].replace("D", "E"))
                            TQEz = float(file_line[54:69].replace("D", "E"))
                            file_line = oclu_file.readline()
                            TQSx = float(file_line[8:23].replace("D", "E"))
                            TQSy = float(file_line[31:46].replace("D", "E"))
                            TQSz = float(file_line[54:69].replace("D", "E"))
                            # we can write the RAPRS values
                            output_line = "{0:.7E},{1:.3E},{2:.3E},{3:.3E},{4:.3E},{5:.7E},{6:.7E}\n".format(alam, tidg, tsdg, pidg, psdg, cosav, rapr)
                            output_line = "{0:.7E},{1:.3E},{2:.3E},{3:.3E},{4:.3E},{5:.7E},{6:.7E},{7:.7E},{8:.7E},{9:.7E},{10:.7E},{11:.7E},{12:.7E},{13:.7E},{14:.7E},{15:.7E},{16:.7E},{17:.7E},{18:.7E},{19:.7E},{20:.7E},{21:.7E},{22:.7E},{23:.7E},{24:.7E}\n".format(alam, tidg, tsdg, pidg, psdg, cosav, rapr, fl, fr, fk, fx, fy, fx, TQEl, TQEr, TQEk, TQEx, TQEy, TQEz, TQSl, TQSr, TQSk, TQSx, TQSy, TQSz)
                            if (out51 is not None): out51.write(output_line)
                            # we know the differential values for polarization
                            # state 1 are after 9 more lines
                            for i in range(9):
                            # state 1 are after 3 more lines
                            for i in range(3):
                                file_line = oclu_file.readline()
                                # the following check is needed to parse C++ output
                                if ("INSERTION" in file_line):
@@ -294,8 +319,33 @@ def parse_legacy_oclu(config):
                                    file_line = oclu_file.readline()
                            cosav = float(file_line[8:23].replace("D", "E"))
                            rapr = float(file_line[31:46].replace("D", "E"))
                            # we read the forces and torques
                            file_line = oclu_file.readline()
                            fl = float(file_line[5:20].replace("D", "E"))
                            fr = float(file_line[25:40].replace("D", "E"))
                            fk = float(file_line[45:60].replace("D", "E"))
                            file_line = oclu_file.readline()
                            fx = float(file_line[5:20].replace("D", "E"))
                            fy = float(file_line[25:40].replace("D", "E"))
                            fz = float(file_line[45:60].replace("D", "E"))
                            file_line = oclu_file.readline()
                            TQEl = float(file_line[8:23].replace("D", "E"))
                            TQEr = float(file_line[31:46].replace("D", "E"))
                            TQEk = float(file_line[54:69].replace("D", "E"))
                            file_line = oclu_file.readline()
                            TQSl = float(file_line[8:23].replace("D", "E"))
                            TQSr = float(file_line[31:46].replace("D", "E"))
                            TQSk = float(file_line[54:69].replace("D", "E"))
                            file_line = oclu_file.readline()
                            TQEx = float(file_line[8:23].replace("D", "E"))
                            TQEy = float(file_line[31:46].replace("D", "E"))
                            TQEz = float(file_line[54:69].replace("D", "E"))
                            file_line = oclu_file.readline()
                            TQSx = float(file_line[8:23].replace("D", "E"))
                            TQSy = float(file_line[31:46].replace("D", "E"))
                            TQSz = float(file_line[54:69].replace("D", "E"))
                            # we can write the RAPRS values
                            output_line = "{0:.7E},{1:.3E},{2:.3E},{3:.3E},{4:.3E},{5:.7E},{6:.7E}\n".format(alam, tidg, tsdg, pidg, psdg, cosav, rapr)
                            output_line = "{0:.7E},{1:.3E},{2:.3E},{3:.3E},{4:.3E},{5:.7E},{6:.7E},{7:.7E},{8:.7E},{9:.7E},{10:.7E},{11:.7E},{12:.7E},{13:.7E},{14:.7E},{15:.7E},{16:.7E},{17:.7E},{18:.7E},{19:.7E},{20:.7E},{21:.7E},{22:.7E},{23:.7E},{24:.7E}\n".format(alam, tidg, tsdg, pidg, psdg, cosav, rapr, fl, fr, fk, fx, fy, fx, TQEl, TQEr, TQEk, TQEx, TQEy, TQEz, TQSl, TQSr, TQSk, TQSx, TQSy, TQSz)
                            if (out52 is not None): out52.write(output_line)
                    found_differentials = True # terminate the inner loop
            # The parsing loop ends here
Loading