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

Optimize legacy configuration binary I/O format

parent d058d28e

src/libnptm/Configuration.cpp

+127 −155
Original line number Diff line number Diff line
@@ -234,21 +234,21 @@ ScattererConfiguration::ScattererConfiguration( 
}



ScattererConfiguration::~ScattererConfiguration() {

  int max_ici = 0;

  int configurations = 0;

  for (int i = 1; i <= number_of_spheres; i++) {

    if (iog_vec[i - 1] < i) continue;

    int nsh = nshl_vec[i - 1];

    if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;

    configurations++;

  }

  for (int i = 0; i < max_ici; i++) {

  for (int i = 0; i < configurations; i++) {

    for (int j = 0; j < number_of_spheres; j++) {

      delete[] dc0_matrix[i][j];

    }

  }

  delete[] dc0_matrix;

  for (int i = 0; i < number_of_spheres; i++) {

  for (int i = 0; i < configurations; i++) {

    delete[] rcf[i];

  }

  delete[] rcf;

  delete[] nshl_vec;

  delete[] radii_of_spheres;

  delete[] iog_vec;
@@ -430,88 +430,82 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam 
    }

  }

  int *iog_vector = new int[nsph]();

  double *ros_vector = new double[nsph]();

  double **rcf_vector = new double*[nsph];

  int *nshl_vector = new int[nsph]();

  int filled_iogs = 0;

  str_target = file_lines[++last_read_line];

  re = regex("[0-9]+");

  while (filled_iogs < nsph) {

    string str_target = file_lines[++last_read_line];

    while(regex_search(str_target, m, re)) {

      iog_vector[filled_iogs++] = stoi(m.str());

  int configurations = 0;

  for (int i = 1; i <= nsph; i++) {

    bool success = regex_search(str_target, m, re);

    if (success) {

      iog_vector[i - 1] = stoi(m.str());

      str_target = m.suffix().str();

      if (filled_iogs == nsph) break;

      if (iog_vector[i - 1] >= i) configurations++;

    } else {

      str_target = file_lines[++last_read_line];

      i--;

    }

  }

  int index = 0;

  double *ros_vector = new double[configurations]();

  int *nshl_vector = new int[configurations]();

  double **rcf_vector = new double*[configurations];

  for (int i113 = 1; i113 <= nsph; i113++) {

    int i_val, nsh;

    double ros_val;

    if (iog_vector[i113 - 1] < i113) {

      rcf_vector[i113 - 1] = new double[1]();

      continue;

    }

    re = regex("[0-9]+");

    if (iog_vector[i113 - 1] < i113) continue;

    str_target = file_lines[++last_read_line];

    regex_search(str_target, m, re);

    i_val = stoi(m.str());

    str_target = m.suffix();

    nshl_vector[i113 - 1] = stoi(m.str());

    str_target = m.suffix().str();

    re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");

    regex_search(str_target, m, re);

    string str_number = m.str();

    str_number = regex_replace(str_number, regex("D"), "e");

    str_number = regex_replace(str_number, regex("d"), "e");

    ros_val = stod(str_number);

    nshl_vector[i113 - 1] = i_val;

    if (max_ici < (i_val + 1) / 2) max_ici = (i_val + 1) / 2;

    ros_vector[i113 - 1] = ros_val;

    nsh = nshl_vector[i113 - 1];

    ros_vector[i113 - 1] = stod(m.str());

    int nsh = nshl_vector[i113 - 1];

    if (i113 == 1) nsh += ies;

    rcf_vector[i113 - 1] = new double[nsh]();

    for (int ns = 0; ns < nsh; ns++) {

      double ns_rcf;

    for (int ns112 = 0; ns112 < nsh; ns112++) {

      str_target = file_lines[++last_read_line];

      regex_search(str_target, m, re);

      str_number = m.str();

      str_number = regex_replace(str_number, regex("D"), "e");

      str_number = regex_replace(str_number, regex("d"), "e");

      ns_rcf = stod(str_number);

      rcf_vector[i113 - 1][ns] = ns_rcf;

    }

  }

  complex<double> ***dc0m = new complex<double>**[max_ici];

  int dim3 = (_idfc == 0) ? nxi : 1;

  for (int dim1 = 0; dim1 < max_ici; dim1++) {

    dc0m[dim1] = new complex<double>*[nsph];

    for (int dim2 = 0; dim2 < nsph; dim2++) {

      dc0m[dim1][dim2] = new complex<double>[dim3]();

    }

  }

  re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");

      rcf_vector[i113 - 1][ns112] = stod(str_number);

    } // ns112 loop

  } // i113 loop



  complex<double> ***dc0m = new complex<double>**[configurations];

  complex<double> *dc0 = new complex<double>[configurations]();

  for (int di = 0; di < configurations; di++) {

    dc0m[di] = new complex<double>*[nsph];

    for (int dj = 0; dj < nsph; dj++) dc0m[di][dj] = new complex<double>[nxi]();

  } // di loop

  for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {

    if (_idfc != 0 && jxi468 > 1) continue;

    if (_idfc != 0 && jxi468 > 1) continue; // jxi468 loop

    for (int i162 = 1; i162 <= nsph; i162++) {

      if (iog_vector[i162 - 1] < i162) continue;

      if (iog_vector[i162 - 1] >= i162) {

	int nsh = nshl_vector[i162 - 1];

      int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?

      if (i162 == 1) ici = ici + ies;

      for (int i157 = 0; i157 < ici; i157++) {

	double dc0_real, dc0_img;

	int ici = (nsh + 1) / 2;

	if (i162 == 1) ici += ies;

	for (int ic157 = 0; ic157 < ici; ic157++) {

	  str_target = file_lines[++last_read_line];

	  regex_search(str_target, m, re);

	  string str_number = m.str();

	  str_number = regex_replace(str_number, regex("D"), "e");

	  str_number = regex_replace(str_number, regex("d"), "e");

	dc0_real = stod(str_number);

	  double rval = stod(str_number);

	  str_target = m.suffix().str();

	  regex_search(str_target, m, re);

	  str_number = m.str();

	  str_number = regex_replace(str_number, regex("D"), "e");

	  str_number = regex_replace(str_number, regex("d"), "e");

	dc0_img = stod(str_number);

	dc0m[i157][i162 - 1][jxi468 - 1] = complex<double>(dc0_real, dc0_img);

      }

    }

  }

	  double ival = stod(str_number);

	  dc0[ic157] = complex<double>(rval, ival);

	  dc0m[ic157][i162 - 1][jxi468 - 1] = dc0[ic157];

	} // ic157 loop

      }

    } // i162 loop

  } // jxi468 loop

  delete[] dc0;



  ScattererConfiguration *config = new ScattererConfiguration(

							      nsph,
@@ -627,104 +621,83 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) { 
}



ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {

  int nsph;

  int *iog;

  int _nsph;

  int *_iog_vec;

  int _ies;

  double _exdc, _wp, _xip;

  int _idfc, nxi;

  int *nshl_vector;

  double *xi_vec;

  double *ros_vector;

  double **rcf_vector;

  complex<double> ***dc0m;

  int _idfc, _nxi;

  int *_nshl_vec;

  double *_xi_vec;

  double *_ros_vec;

  double **_rcf_vec;

  complex<double> ***_dc0m;



  int max_ici = 0;

  fstream input;

  input.open(file_name.c_str(), ios::in | ios::binary);

  if (input.is_open()) {

    input.read(reinterpret_cast<char *>(&nsph), sizeof(int));

    iog = new int[nsph]();

    for (int i = 0; i < nsph; i++) {

      input.read(reinterpret_cast<char *>(&(iog[i])), sizeof(int));

    }

    int configuration_count = 0;

    for (int i115 = 1; i115 <= nsph; i115++) {

      if (iog[i115 - 1] < i115) continue;

      configuration_count++;

    }

  input.read(reinterpret_cast<char *>(&_nsph), sizeof(int));

  input.read(reinterpret_cast<char *>(&_ies), sizeof(int));

  _iog_vec = new int[_nsph]();

  for (int i = 0; i < _nsph; i++)

    input.read(reinterpret_cast<char *>(&(_iog_vec[i])), sizeof(int));

  input.read(reinterpret_cast<char *>(&_exdc), sizeof(double));

  input.read(reinterpret_cast<char *>(&_wp), sizeof(double));

  input.read(reinterpret_cast<char *>(&_xip), sizeof(double));

  input.read(reinterpret_cast<char *>(&_idfc), sizeof(int));

    input.read(reinterpret_cast<char *>(&nxi), sizeof(int));

    try {

      xi_vec = new double[nxi]();

    } catch (const bad_alloc &ex) {

      throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of scales " + to_string(nxi));

    }

    for (int i = 0; i < nxi; i++) {

      input.read(reinterpret_cast<char *>(&(xi_vec[i])), sizeof(double));

    }

    nshl_vector = new int[configuration_count]();

    ros_vector = new double[configuration_count]();

    rcf_vector = new double*[configuration_count];

    for (int i115 = 1; i115 <= nsph; i115++) {

      if (iog[i115 - 1] < i115) {

	rcf_vector[i115 - 1] = new double[1]();

	continue;

      }

      input.read(reinterpret_cast<char *>(&(nshl_vector[i115 - 1])), sizeof(int)); // was not from IOG

      input.read(reinterpret_cast<char *>(&(ros_vector[i115 - 1])), sizeof(double)); // was not from IOG

      int nsh = nshl_vector[i115 - 1]; // was not from IOG

      if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;

      try {

	rcf_vector[i115 - 1] = new double[nsh](); // was not from IOG

      } catch (const bad_alloc &ex) {

	throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of layers " + to_string(nsh));

      }

      for (int nsi = 0; nsi < nsh; nsi++) {

	input.read(reinterpret_cast<char *>(&(rcf_vector[i115 - 1][nsi])), sizeof(double)); // was not from IOG

      }

    }

    dc0m = new complex<double>**[max_ici];

    int dim3 = (_idfc == 0) ? nxi : 1;

    for (int dim1 = 0; dim1 < max_ici; dim1++) {

      dc0m[dim1] = new complex<double>*[nsph];

      for (int dim2 = 0; dim2 < nsph; dim2++) {

	dc0m[dim1][dim2] = new complex<double>[dim3]();

      }

  input.read(reinterpret_cast<char *>(&_nxi), sizeof(int));

  _xi_vec = new double[_nxi]();

  for (int i = 0; i < _nxi; i++)

    input.read(reinterpret_cast<char *>(&(_xi_vec[i])), sizeof(double));

  int configurations = 0;

  for (int i = 1; i <= _nsph; i++) {

    if (_iog_vec[i - 1] >= i) configurations++;

  }

  _nshl_vec = new int[configurations]();

  _ros_vec = new double[configurations]();

  _rcf_vec = new double*[configurations];

  for (int i115 = 1; i115 <= _nsph; i115++) {

    if (_iog_vec[i115 - 1] < i115) continue;

    input.read(reinterpret_cast<char *>(&(_nshl_vec[i115 - 1])), sizeof(int));

    input.read(reinterpret_cast<char *>(&(_ros_vec[i115 - 1])), sizeof(double));

    int nsh = _nshl_vec[i115 - 1];

    if (i115 == 1) nsh += _ies;

    _rcf_vec[i115 - 1] = new double[nsh]();

    for (int nsi = 0; nsi < nsh; nsi++)

      input.read(reinterpret_cast<char *>(&(_rcf_vec[i115 - 1][nsi])), sizeof(double));

  }

    for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {

  _dc0m = new complex<double>**[configurations];

  for (int di = 0; di < configurations; di++) {

    _dc0m[di] = new complex<double>*[_nsph];

    for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new complex<double>[_nxi]();

  } // di loop

  for (int jxi468 = 1; jxi468 <= _nxi; jxi468++) {

    if (_idfc != 0 && jxi468 > 1) continue;

      for (int i162 = 1; i162 <= nsph; i162++) {

	if (iog[i162 - 1] < i162) continue;

	int nsh = nshl_vector[iog[i162 - 1] - 1]; // was not from IOG

	int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?

    for (int i162 = 1; i162 <= _nsph; i162++) {

      if (_iog_vec[i162 - 1] < i162) continue;

      int nsh = _nshl_vec[i162 - 1];

      int ici = (nsh + 1) / 2;

      if (i162 == 1) ici = ici + _ies;

      for (int i157 = 0; i157 < ici; i157++) {

	  double dc0_real, dc0_img;

	  input.read(reinterpret_cast<char *>(&dc0_real), sizeof(double));

	  input.read(reinterpret_cast<char *>(&dc0_img), sizeof(double));

	  dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * dc0_img;

	double rval, ival;

	input.read(reinterpret_cast<char *>(&rval), sizeof(double));

	input.read(reinterpret_cast<char *>(&ival), sizeof(double));

	_dc0m[i157][i162 - 1][jxi468 - 1] = complex<double>(rval, ival);

      }

    }

  }

  input.close();

  } else { // Opening of the input file did not succeed.

    OpenConfigurationFileException ex(file_name);

    throw ex;

  }

    

  ScattererConfiguration *conf = new ScattererConfiguration(

							    nsph,

							    xi_vec,

							    nxi,

							    _nsph,

							    _xi_vec,

							    _nxi,

							    "XIV",

							    iog,

							    ros_vector,

							    nshl_vector,

							    rcf_vector,

							    _iog_vec,

							    _ros_vec,

							    _nshl_vec,

							    _rcf_vec,

							    _idfc,

							    dc0m,

							    false,

							    _dc0m,

							    (_ies == 1),

							    _exdc,

							    _wp,

							    _xip
@@ -889,12 +862,11 @@ void ScattererConfiguration::write_hdf5(string file_name) { 
}



void ScattererConfiguration::write_legacy(string file_name) {

  const double two_pi = acos(0.0) * 4.0;

  const double evc = 6.5821188e-16;

  fstream output;

  int ies = (use_external_sphere)? 1 : 0;

  output.open(file_name.c_str(), ios::out | ios::binary);

  output.write(reinterpret_cast<char *>(&number_of_spheres), sizeof(int));

  output.write(reinterpret_cast<char *>(&ies), sizeof(int));

  for (int i = 0; i < number_of_spheres; i++)

    output.write(reinterpret_cast<char *>(&(iog_vec[i])), sizeof(int));

  output.write(reinterpret_cast<char *>(&exdc), sizeof(double));