<commit_before>/*
 * Copyright 2009-2015 The VOTCA Development Team (http://www.votca.org)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

#include "h5mdtrajectoryreader.h"
#include "hdf5.h"

#include <vector>

namespace votca {
namespace csg {

H5MDTrajectoryReader::H5MDTrajectoryReader(){
  has_velocity_ = H5MDTrajectoryReader::NONE;
  has_force_ = H5MDTrajectoryReader::NONE;
  has_id_group_ = H5MDTrajectoryReader::NONE;
  has_box_ = H5MDTrajectoryReader::NONE;
}

H5MDTrajectoryReader::~H5MDTrajectoryReader() {
  if (file_opened_) {
    H5Fclose(file_id_);
    file_opened_ = false;
  }
}

bool H5MDTrajectoryReader::Open(const string &file) {
  // Checks if we deal with hdf5 file.
  if (!H5Fis_hdf5(file.c_str())) {
    cout << file << " is not recognise as HDF5 file format" << endl;
    return false;
  }
  file_id_ = H5Fopen(file.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
  file_opened_ = true;

  // Check the version of the file.
  hid_t g_h5md = H5Gopen(file_id_, "h5md", H5P_DEFAULT);
  CheckError(g_h5md, "Unable to open /h5md group.");

  hid_t at_version = H5Aopen(g_h5md, "version", H5P_DEFAULT);
  CheckError(at_version, "Unable to read version attribute.");
  int version[2];
  H5Aread(at_version, H5Aget_type(at_version), &version);
  if (version[0] != 1 || (version[0] == 1 && version[1] > 0)) {
    cout << "Major version " << version[0] << endl;
    cout << "Minor version " << version[1] << endl;
    throw ios_base::failure("Wrong version of H5MD file.");
  }

  // Clean up.
  H5Aclose(at_version);
  H5Gclose(g_h5md);

  // Checks if the particles group exists and what is the number of members.
  particle_group_ = H5Gopen(file_id_, "particles", H5P_DEFAULT);
  CheckError(particle_group_, "Unable to open /particles group.");
  hsize_t num_obj = 0;
  H5Gget_num_objs(particle_group_, &num_obj);
  if (num_obj == 0) {
    throw ios_base::failure("The particles group is empty.");
  }

  first_frame_ = true;

  // Handle errors by internal check up.
  H5Eset_auto2(H5E_DEFAULT, NULL, NULL);

  return true;
}

void H5MDTrajectoryReader::Close() {
  if (file_opened_) {
    H5Fclose(file_id_);
    file_opened_ = false;
  }
}

void H5MDTrajectoryReader::Initialize(Topology &top) {
  string *particle_group_name_ = new string(top.getParticleGroup());
  if(*particle_group_name_ == "")
    throw ios_base::failure(
        "Missing particle group in topology. Please set `h5md_particle_group` tag with `name` attribute set to the particle group.");
  string *position_group_name = new string(*particle_group_name_ + "/position");
  atom_position_group_ = H5Gopen(particle_group_, position_group_name->c_str(), H5P_DEFAULT);
  CheckError(atom_position_group_, "Unable to open " + *position_group_name + " group");

  idx_frame_ = -1;
  ds_atom_position_ = H5Dopen(atom_position_group_, "value", H5P_DEFAULT);
  CheckError(ds_atom_position_, "Unable to open " + *position_group_name + "/value dataset");

  // Reads the box information.
  string *box_gr_name = new string(*particle_group_name_ + "/box");
  hid_t g_box = H5Gopen(particle_group_, box_gr_name->c_str(), H5P_DEFAULT);
  CheckError(g_box, "Unable to open " + *box_gr_name + " group");
  hid_t at_box_dimension = H5Aopen(g_box, "dimension", H5P_DEFAULT);
  CheckError(at_box_dimension, "Unable to open dimension attribute.");
  int dimension;
  H5Aread(at_box_dimension, H5Aget_type(at_box_dimension), &dimension);
  if (dimension != 3 ) {
    throw ios_base::failure(
        "Wrong dimension " + boost::lexical_cast<string>(dimension));
  }
  // TODO: check if boundary is periodic.
  string *box_edges_name = new string(*particle_group_name_ + "/box/edges");
  if (GroupExists(particle_group_, *box_edges_name)) {
    g_box = H5Gopen(particle_group_, box_gr_name->c_str(), H5P_DEFAULT);
    edges_group_ = H5Gopen(g_box, "edges", H5P_DEFAULT);
    ds_edges_group_ = H5Dopen(edges_group_, "value", H5P_DEFAULT);
    cout << "H5MD: has /box/edges" << endl;
    has_box_ = H5MDTrajectoryReader::TIMEDEPENDENT;
  } else {
    cout << "H5MD: static box" << endl;
    hid_t ds_edges = H5Dopen(g_box, "edges", H5P_DEFAULT);
    CheckError(ds_edges, "Unable to open /box/edges");
    double *box = new double[3];
    ReadStaticData(ds_edges, H5T_NATIVE_DOUBLE, box);
    cout << "H5MD: Found box " << box[0] << " x " << box[1] <<
        " x " << box[2] << endl;
    // Sets box size.
    m.ZeroMatrix();
    m[0][0] = box[0];
    m[1][1] = box[1];
    m[2][2] = box[2];
    top.setBox(m);
    has_box_ = H5MDTrajectoryReader::STATIC;
  }
  H5Gclose(g_box);

  // Gets the force group.
  string *force_group_name = new string(*particle_group_name_ + "/force");
  if (GroupExists(particle_group_, *force_group_name)) {
    atom_force_group_ = H5Gopen(particle_group_, force_group_name->c_str(), H5P_DEFAULT);
    ds_atom_force_ = H5Dopen(atom_force_group_, "value", H5P_DEFAULT);
    has_force_ = H5MDTrajectoryReader::TIMEDEPENDENT;
    cout << "H5MD: has /force" << endl;
  } else {
    has_force_ = H5MDTrajectoryReader::NONE;
  }

  // Gets the velocity group.
  string *velocity_group_name = new string(*particle_group_name_ + "/velocity");
  if (GroupExists(particle_group_, *velocity_group_name)) {
    atom_velocity_group_ = H5Gopen(particle_group_,
                                   velocity_group_name->c_str(),
                                   H5P_DEFAULT);
    ds_atom_velocity_ = H5Dopen(atom_velocity_group_, "value", H5P_DEFAULT);
    has_velocity_ = H5MDTrajectoryReader::TIMEDEPENDENT;
    cout << "H5MD: has /velocity" << endl;
  } else {
    has_velocity_ = H5MDTrajectoryReader::NONE;
  }

  // Gets the id group so that the atom id is taken from this group.
  string *id_group_name = new string(*particle_group_name_  + "/id");
  if (GroupExists(particle_group_, *id_group_name)) {
    atom_id_group_ = H5Gopen(particle_group_, id_group_name->c_str(), H5P_DEFAULT);
    ds_atom_id_ = H5Dopen(atom_id_group_, "value", H5P_DEFAULT);
    has_id_group_ = H5MDTrajectoryReader::TIMEDEPENDENT;
    cout << "H5MD: has /id group" << endl;
  } else {
    has_id_group_ = H5MDTrajectoryReader::NONE;
  }

  // Gets number of particles and dimensions.
  hid_t fs_atom_position_ = H5Dget_space(ds_atom_position_);
  CheckError(fs_atom_position_, "Unable to open atom position space.");
  hsize_t dims[3];
  rank_ = H5Sget_simple_extent_dims(fs_atom_position_, dims, NULL);
  N_particles_ = dims[1];
  vec_components_ = dims[2];
  max_idx_frame_ = dims[0] - 1;

  // TODO: reads mass, charge and particle type.

  if (!(has_id_group_ != NONE) && top.BeadCount() > 0 && N_particles_ != top.BeadCount()) {
    cout << "Warning: The number of beads (" << N_particles_ << ")";
    cout << " in the trajectory is different than defined in the topology ("
        << top.BeadCount() << ")" << endl;<commit_msg>small fix, check on the number of particles<commit_after>/*
 * Copyright 2009-2015 The VOTCA Development Team (http://www.votca.org)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

#include "h5mdtrajectoryreader.h"
#include "hdf5.h"

#include <vector>

namespace votca {
namespace csg {

H5MDTrajectoryReader::H5MDTrajectoryReader(){
  has_velocity_ = H5MDTrajectoryReader::NONE;
  has_force_ = H5MDTrajectoryReader::NONE;
  has_id_group_ = H5MDTrajectoryReader::NONE;
  has_box_ = H5MDTrajectoryReader::NONE;
}

H5MDTrajectoryReader::~H5MDTrajectoryReader() {
  if (file_opened_) {
    H5Fclose(file_id_);
    file_opened_ = false;
  }
}

bool H5MDTrajectoryReader::Open(const string &file) {
  // Checks if we deal with hdf5 file.
  if (!H5Fis_hdf5(file.c_str())) {
    cout << file << " is not recognise as HDF5 file format" << endl;
    return false;
  }
  file_id_ = H5Fopen(file.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
  file_opened_ = true;

  // Check the version of the file.
  hid_t g_h5md = H5Gopen(file_id_, "h5md", H5P_DEFAULT);
  CheckError(g_h5md, "Unable to open /h5md group.");

  hid_t at_version = H5Aopen(g_h5md, "version", H5P_DEFAULT);
  CheckError(at_version, "Unable to read version attribute.");
  int version[2];
  H5Aread(at_version, H5Aget_type(at_version), &version);
  if (version[0] != 1 || (version[0] == 1 && version[1] > 0)) {
    cout << "Major version " << version[0] << endl;
    cout << "Minor version " << version[1] << endl;
    throw ios_base::failure("Wrong version of H5MD file.");
  }

  // Clean up.
  H5Aclose(at_version);
  H5Gclose(g_h5md);

  // Checks if the particles group exists and what is the number of members.
  particle_group_ = H5Gopen(file_id_, "particles", H5P_DEFAULT);
  CheckError(particle_group_, "Unable to open /particles group.");
  hsize_t num_obj = 0;
  H5Gget_num_objs(particle_group_, &num_obj);
  if (num_obj == 0) {
    throw ios_base::failure("The particles group is empty.");
  }

  first_frame_ = true;

  // Handle errors by internal check up.
  H5Eset_auto2(H5E_DEFAULT, NULL, NULL);

  return true;
}

void H5MDTrajectoryReader::Close() {
  if (file_opened_) {
    H5Fclose(file_id_);
    file_opened_ = false;
  }
}

void H5MDTrajectoryReader::Initialize(Topology &top) {
  string *particle_group_name_ = new string(top.getParticleGroup());
  if(*particle_group_name_ == "")
    throw ios_base::failure(
        "Missing particle group in topology. Please set `h5md_particle_group` tag with `name` attribute set to the particle group.");
  string *position_group_name = new string(*particle_group_name_ + "/position");
  atom_position_group_ = H5Gopen(particle_group_, position_group_name->c_str(), H5P_DEFAULT);
  CheckError(atom_position_group_, "Unable to open " + *position_group_name + " group");

  idx_frame_ = -1;
  ds_atom_position_ = H5Dopen(atom_position_group_, "value", H5P_DEFAULT);
  CheckError(ds_atom_position_, "Unable to open " + *position_group_name + "/value dataset");

  // Reads the box information.
  string *box_gr_name = new string(*particle_group_name_ + "/box");
  hid_t g_box = H5Gopen(particle_group_, box_gr_name->c_str(), H5P_DEFAULT);
  CheckError(g_box, "Unable to open " + *box_gr_name + " group");
  hid_t at_box_dimension = H5Aopen(g_box, "dimension", H5P_DEFAULT);
  CheckError(at_box_dimension, "Unable to open dimension attribute.");
  int dimension;
  H5Aread(at_box_dimension, H5Aget_type(at_box_dimension), &dimension);
  if (dimension != 3 ) {
    throw ios_base::failure(
        "Wrong dimension " + boost::lexical_cast<string>(dimension));
  }
  // TODO: check if boundary is periodic.
  string *box_edges_name = new string(*particle_group_name_ + "/box/edges");
  if (GroupExists(particle_group_, *box_edges_name)) {
    g_box = H5Gopen(particle_group_, box_gr_name->c_str(), H5P_DEFAULT);
    edges_group_ = H5Gopen(g_box, "edges", H5P_DEFAULT);
    ds_edges_group_ = H5Dopen(edges_group_, "value", H5P_DEFAULT);
    cout << "H5MD: has /box/edges" << endl;
    has_box_ = H5MDTrajectoryReader::TIMEDEPENDENT;
  } else {
    cout << "H5MD: static box" << endl;
    hid_t ds_edges = H5Dopen(g_box, "edges", H5P_DEFAULT);
    CheckError(ds_edges, "Unable to open /box/edges");
    double *box = new double[3];
    ReadStaticData(ds_edges, H5T_NATIVE_DOUBLE, box);
    cout << "H5MD: Found box " << box[0] << " x " << box[1] <<
        " x " << box[2] << endl;
    // Sets box size.
    m.ZeroMatrix();
    m[0][0] = box[0];
    m[1][1] = box[1];
    m[2][2] = box[2];
    top.setBox(m);
    has_box_ = H5MDTrajectoryReader::STATIC;
  }
  H5Gclose(g_box);

  // Gets the force group.
  string *force_group_name = new string(*particle_group_name_ + "/force");
  if (GroupExists(particle_group_, *force_group_name)) {
    atom_force_group_ = H5Gopen(particle_group_, force_group_name->c_str(), H5P_DEFAULT);
    ds_atom_force_ = H5Dopen(atom_force_group_, "value", H5P_DEFAULT);
    has_force_ = H5MDTrajectoryReader::TIMEDEPENDENT;
    cout << "H5MD: has /force" << endl;
  } else {
    has_force_ = H5MDTrajectoryReader::NONE;
  }

  // Gets the velocity group.
  string *velocity_group_name = new string(*particle_group_name_ + "/velocity");
  if (GroupExists(particle_group_, *velocity_group_name)) {
    atom_velocity_group_ = H5Gopen(particle_group_,
                                   velocity_group_name->c_str(),
                                   H5P_DEFAULT);
    ds_atom_velocity_ = H5Dopen(atom_velocity_group_, "value", H5P_DEFAULT);
    has_velocity_ = H5MDTrajectoryReader::TIMEDEPENDENT;
    cout << "H5MD: has /velocity" << endl;
  } else {
    has_velocity_ = H5MDTrajectoryReader::NONE;
  }

  // Gets the id group so that the atom id is taken from this group.
  string *id_group_name = new string(*particle_group_name_  + "/id");
  if (GroupExists(particle_group_, *id_group_name)) {
    atom_id_group_ = H5Gopen(particle_group_, id_group_name->c_str(), H5P_DEFAULT);
    ds_atom_id_ = H5Dopen(atom_id_group_, "value", H5P_DEFAULT);
    has_id_group_ = H5MDTrajectoryReader::TIMEDEPENDENT;
    cout << "H5MD: has /id group" << endl;
  } else {
    has_id_group_ = H5MDTrajectoryReader::NONE;
  }

  // Gets number of particles and dimensions.
  hid_t fs_atom_position_ = H5Dget_space(ds_atom_position_);
  CheckError(fs_atom_position_, "Unable to open atom position space.");
  hsize_t dims[3];
  rank_ = H5Sget_simple_extent_dims(fs_atom_position_, dims, NULL);
  N_particles_ = dims[1];
  vec_components_ = dims[2];
  max_idx_frame_ = dims[0] - 1;

  // TODO: reads mass, charge and particle type.

  if (has_id_group_ == H5MDTrajectoryReader::NONE && top.BeadCount() > 0 && N_particles_ != top.BeadCount()) {
    cout << "Warning: The number of beads (" << N_particles_ << ")";
    cout << " in the trajectory is different than defined in the topology ("
        << top.BeadCount() << ")" << endl;
{'fraction_non_alphanumeric': 0.12522995162499148, 'fraction_numerical': 0.018600531443755536, 'mean_word_length': 4.204964539007092, 'pattern_counts': {'":': 0, '<': 139, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2579', 'n_tokens_mistral': 5367, 'n_tokens_neox': 5223, 'n_words': 1210}