<commit_before>/*
 *            Copyright 2016 The MUSCET Development Team
 *
 *      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 <regex>
#include <votca/xtp/segmentmapper.h>

namespace votca {
namespace xtp {

template <class AtomContainer>
SegmentMapper<AtomContainer>::SegmentMapper(Logger& log) : _log(log) {
  FillMap();
};

template <class AtomContainer>
AtomContainer SegmentMapper<AtomContainer>::map(const Segment& seg,
                                                const SegId& segid) const {
  if (segid.hasFile()) {
    std::string filename = segid.FileName();
    filename = std::regex_replace(filename, std::regex("\\$SEGID"),
                                  std::to_string(seg.getId()));
    filename =
        std::regex_replace(filename, std::regex("\\$SEGNAME"), seg.getName());
    return map(seg, filename);
  } else {
    QMState state = segid.getQMState();
    return map(seg, state);
  }
}

template <class AtomContainer>
template <typename T>
Eigen::Vector3d SegmentMapper<AtomContainer>::CalcWeightedPos(
    const std::vector<double>& weights, const T& atoms) const {
  Eigen::Vector3d map_pos = Eigen::Vector3d::Zero();
  double map_weight = 0.0;
  for (unsigned i = 0; i < atoms.size(); i++) {
    map_pos += atoms[i]->getPos() * weights[i];
    map_weight += weights[i];
  }
  return map_pos = map_pos / map_weight;
}

template <class AtomContainer>
void SegmentMapper<AtomContainer>::ParseFragment(Seginfo& seginfo,
                                                 const tools::Property& frag) {
  tools::Tokenizer tok_map_atoms(
      frag.get(_mapatom_xml["atoms"]).as<std::string>(), " \t\n");
  std::vector<std::string> map_atoms = tok_map_atoms.ToVector();
  tools::Tokenizer tok_md_atoms(frag.get("mdatoms").as<std::string>(), " \t\n");
  std::vector<std::string> md_atoms = tok_md_atoms.ToVector();

  if (md_atoms.size() != map_atoms.size()) {
    throw std::runtime_error(
        "Mapping for segment " + seginfo.segname + " fragment " +
        frag.get("name").as<std::string>() +
        " does not have same numbers of md and " + _mapatom_xml["atoms"] +
        "."
        "If you want to leave a qmatom out, place a ':' instead");
  }

  tools::Tokenizer tok_weights(getWeights(frag), " \t\n");
  std::vector<double> weights;
  tok_weights.ConvertToVector(weights);

  if (md_atoms.size() != weights.size()) {
    throw std::runtime_error("Mapping for segment " + seginfo.segname +
                             " fragment " + frag.get("name").as<std::string>() +
                             " does not have same numbers of md and weights. "
                             "If you want to leave a " +
                             _mapatom_xml["name"] +
                             " out, place a '0' instead");
  }

  FragInfo mapfragment;
  std::vector<int> mapatom_ids;

  for (unsigned i = 0; i < map_atoms.size(); i++) {
    const std::string& map_string = map_atoms[i];
    const std::string& md_string = md_atoms[i];
    const double& weight = weights[i];
    MD_atom_id md_result = StringToMDIndex(md_string);
    seginfo.mdatoms.push_back(md_result);

    if (map_string == ":") {
      continue;
    }
    mapsite_id map_result = StringToMapIndex(map_string);
    mapatom_ids.push_back(map_result.first);
    seginfo.mapatoms.push_back(map_result);
    if (Atom::GetElementFromMDName(md_result.second) != map_result.second) {
      XTP_LOG_SAVE(logINFO, _log)
          << "WARNING: mdatom'" << md_result.second << "' and "
          << _mapatom_xml["name"] << " '" << map_result.second
          << "' do not have same element" << std::flush;
    }
    mapfragment._mapatom_ids.push_back(map_result);
    mapfragment._weights.push_back(weight);
    mapfragment._mdatom_ids.push_back(md_result);
  }

  tools::Tokenizer tok_frame(getFrame(frag), " \t\n");
  std::vector<int> frame;
  tok_frame.ConvertToVector(frame);
  if (frame.size() > 3) {
    throw std::runtime_error(
        "Local frame for segment " + seginfo.segname + " fragment " +
        frag.get("name").as<std::string>() +
        " has more than 3 atoms, please specify only up to three atoms");
  } else if (frame.empty()) {
    throw std::runtime_error("No local frame for segment " + seginfo.segname +
                             " fragment " + frag.get("name").as<std::string>() +
                             " specified");
  }
  for (int atomid : frame) {
    if (std::find(mapatom_ids.begin(), mapatom_ids.end(), atomid) ==
        mapatom_ids.end()) {
      throw std::runtime_error("Atom " + std::to_string(atomid) +
                               " in local frame cannot be found in " +
                               _mapatom_xml["atoms"] + ".");
    }
  }

  mapfragment._map_local_frame = frame;
  seginfo.fragments.push_back(mapfragment);
}

template <class AtomContainer>
void SegmentMapper<AtomContainer>::LoadMappingFile(const std::string& mapfile) {

  tools::Property topology_map;
  tools::load_property_from_xml(topology_map, mapfile);

  std::string molkey = "topology.molecules.molecule";
  std::vector<tools::Property*> molecules = topology_map.Select(molkey);
  std::string segkey = "segments.segment";
  for (tools::Property* mol : molecules) {
    std::vector<tools::Property*> segments = mol->Select(segkey);
    for (tools::Property* seg : segments) {
      Seginfo seginfo;

      std::string coordfile_key = _mapatom_xml["coords"] + "_*";
      std::vector<tools::Property*> files = seg->Select(coordfile_key);
      for (tools::Property* file : files) {
        seginfo.coordfiles[file->name()] = file->as<std::string>();
      }
      std::string segname = seg->get("name").as<std::string>();
      seginfo.segname = segname;
      std::string fragkey = "fragments.fragment";

      seginfo.map2md = seg->ifExistsReturnElseThrowRuntimeError<bool>("map2md");

      std::vector<tools::Property*> fragments = seg->Select(fragkey);
      for (tools::Property* frag : fragments) {
        ParseFragment(seginfo, *frag);
      }

      int map_atom_min_id =
          std::min_element(seginfo.mapatoms.begin(), seginfo.mapatoms.end(),
                           [](const mapsite_id& a, const mapsite_id& b) {
                             return a.first < b.first;
                           })
              ->first;
      if (map_atom_min_id != 0) {
        throw std::runtime_error(_mapatom_xml["atoms"] + " for segment " +
                                 seginfo.segname +
                                 " do not start at zero index. Each segment "
                                 "should have its own coordinate file");<commit_msg>added comment in cpp code for python tool<commit_after>/*
 *            Copyright 2016 The MUSCET Development Team
 *
 *      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 <regex>
#include <votca/xtp/segmentmapper.h>

namespace votca {
namespace xtp {

template <class AtomContainer>
SegmentMapper<AtomContainer>::SegmentMapper(Logger& log) : _log(log) {
  FillMap();
};

template <class AtomContainer>
AtomContainer SegmentMapper<AtomContainer>::map(const Segment& seg,
                                                const SegId& segid) const {
  if (segid.hasFile()) {
    std::string filename = segid.FileName();
    filename = std::regex_replace(filename, std::regex("\\$SEGID"),
                                  std::to_string(seg.getId()));
    filename =
        std::regex_replace(filename, std::regex("\\$SEGNAME"), seg.getName());
    return map(seg, filename);
  } else {
    QMState state = segid.getQMState();
    return map(seg, state);
  }
}

template <class AtomContainer>
template <typename T>
Eigen::Vector3d SegmentMapper<AtomContainer>::CalcWeightedPos(
    const std::vector<double>& weights, const T& atoms) const {
  Eigen::Vector3d map_pos = Eigen::Vector3d::Zero();
  double map_weight = 0.0;
  for (unsigned i = 0; i < atoms.size(); i++) {
    map_pos += atoms[i]->getPos() * weights[i];
    map_weight += weights[i];
  }
  return map_pos = map_pos / map_weight;
}

template <class AtomContainer>
void SegmentMapper<AtomContainer>::ParseFragment(Seginfo& seginfo,
                                                 const tools::Property& frag) {
  tools::Tokenizer tok_map_atoms(
      frag.get(_mapatom_xml["atoms"]).as<std::string>(), " \t\n");
  std::vector<std::string> map_atoms = tok_map_atoms.ToVector();
  tools::Tokenizer tok_md_atoms(frag.get("mdatoms").as<std::string>(), " \t\n");
  std::vector<std::string> md_atoms = tok_md_atoms.ToVector();

  if (md_atoms.size() != map_atoms.size()) {
    throw std::runtime_error(
        "Mapping for segment " + seginfo.segname + " fragment " +
        frag.get("name").as<std::string>() +
        " does not have same numbers of md and " + _mapatom_xml["atoms"] +
        "."
        "If you want to leave a qmatom out, place a ':' instead");
  }

  tools::Tokenizer tok_weights(getWeights(frag), " \t\n");
  std::vector<double> weights;
  tok_weights.ConvertToVector(weights);

  if (md_atoms.size() != weights.size()) {
    throw std::runtime_error("Mapping for segment " + seginfo.segname +
                             " fragment " + frag.get("name").as<std::string>() +
                             " does not have same numbers of md and weights. "
                             "If you want to leave a " +
                             _mapatom_xml["name"] +
                             " out, place a '0' instead");
  }

  FragInfo mapfragment;
  std::vector<int> mapatom_ids;

  for (unsigned i = 0; i < map_atoms.size(); i++) {
    const std::string& map_string = map_atoms[i];
    const std::string& md_string = md_atoms[i];
    const double& weight = weights[i];
    MD_atom_id md_result = StringToMDIndex(md_string);
    seginfo.mdatoms.push_back(md_result);

    if (map_string == ":") {
      continue;
    }
    mapsite_id map_result = StringToMapIndex(map_string);
    mapatom_ids.push_back(map_result.first);
    seginfo.mapatoms.push_back(map_result);
    if (Atom::GetElementFromMDName(md_result.second) != map_result.second) {
      XTP_LOG_SAVE(logINFO, _log)
          << "WARNING: mdatom'" << md_result.second << "' and "
          << _mapatom_xml["name"] << " '" << map_result.second
          << "' do not have same element" << std::flush;
    }
    mapfragment._mapatom_ids.push_back(map_result);
    mapfragment._weights.push_back(weight);
    mapfragment._mdatom_ids.push_back(md_result);
  }

  tools::Tokenizer tok_frame(getFrame(frag), " \t\n");
  std::vector<int> frame;
  tok_frame.ConvertToVector(frame);
  if (frame.size() > 3) {
    throw std::runtime_error(
        "Local frame for segment " + seginfo.segname + " fragment " +
        frag.get("name").as<std::string>() +
        " has more than 3 atoms, please specify only up to three atoms");
  } else if (frame.empty()) {
    throw std::runtime_error("No local frame for segment " + seginfo.segname +
                             " fragment " + frag.get("name").as<std::string>() +
                             " specified");
  }
  for (int atomid : frame) {
    if (std::find(mapatom_ids.begin(), mapatom_ids.end(), atomid) ==
        mapatom_ids.end()) {
      throw std::runtime_error("Atom " + std::to_string(atomid) +
                               " in local frame cannot be found in " +
                               _mapatom_xml["atoms"] + ".");
    }
  }

  mapfragment._map_local_frame = frame;
  seginfo.fragments.push_back(mapfragment);
}

template <class AtomContainer>
void SegmentMapper<AtomContainer>::LoadMappingFile(const std::string& mapfile) {

  tools::Property topology_map;
  tools::load_property_from_xml(topology_map, mapfile);

  std::string molkey = "topology.molecules.molecule";
  std::vector<tools::Property*> molecules = topology_map.Select(molkey);
  std::string segkey = "segments.segment";
  for (tools::Property* mol : molecules) {
    std::vector<tools::Property*> segments = mol->Select(segkey);
    for (tools::Property* seg : segments) {
      Seginfo seginfo;

      std::string coordfile_key = _mapatom_xml["coords"] + "_*";
      std::vector<tools::Property*> files = seg->Select(coordfile_key);
      for (tools::Property* file : files) {
        seginfo.coordfiles[file->name()] = file->as<std::string>();
      }
      std::string segname = seg->get("name").as<std::string>();
      seginfo.segname = segname;
      std::string fragkey = "fragments.fragment";

      seginfo.map2md = seg->ifExistsReturnElseThrowRuntimeError<bool>("map2md");

      std::vector<tools::Property*> fragments = seg->Select(fragkey);
      for (tools::Property* frag : fragments) {
        ParseFragment(seginfo, *frag);
      }

      int map_atom_min_id =
          std::min_element(seginfo.mapatoms.begin(), seginfo.mapatoms.end(),
                           [](const mapsite_id& a, const mapsite_id& b) {
                             return a.first < b.first;
                           })
              ->first;
      if (map_atom_min_id != 0) {
        throw std::runtime_error(
            _mapatom_xml["atoms"] + " for segment " + seginfo.segname +
            " do not start at zero index. Each segment "
            "should have its own coordinate file. If you use an old ctp "
            "mapping file run 'xtp_update_mapfile' on it.");
{'fraction_non_alphanumeric': 0.13685608443419306, 'fraction_numerical': 0.002935625917383099, 'mean_word_length': 2.7901986754966885, 'pattern_counts': {'":': 2, '<': 105, '<?xml version=': 0, '>': 87, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 16}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3683681', 'n_tokens_mistral': 4423, 'n_tokens_neox': 4354, 'n_words': 985}