Orbit_impl.hh

Go to the documentation of this file.

#ifndef CASM_Orbit_impl
#define CASM_Orbit_impl
 
#include <boost/iterator/transform_iterator.hpp>
#include <set>
 
#include "casm/crystallography/Lattice.hh"
#include "casm/global/errors.hh"
#include "casm/misc/algorithm.hh"
#include "casm/symmetry/Orbit.hh"
#include "casm/symmetry/SymGroup.hh"
 
// #include "casm/casm_io/Log.hh"
// #include "casm/casm_io/json/jsonParser.hh"
// #include "casm/casm_io/container/json_io.hh"
 
namespace CASM {
 
namespace Orbit_impl {
 
// -- Sort using indices ---
 
struct EqMapRow {
 
  EqMapRow(Index _a, Index _b,
           const std::vector<std::vector<Index>> &tmp_eq_map, const SymGroup &g)
      : a(_a), b(_b) {
    try {
      Index a2proto = g.ind_inverse(tmp_eq_map[a][0]);
      for (Index j = 0; j < tmp_eq_map[b].size(); ++j) {
        values.emplace(g.ind_prod(tmp_eq_map[b][j], a2proto));
      }
    } catch (const std::exception &e) {
      throw libcasm_runtime_error(
          "Error in Orbit constructor: Failed constructing EqMapRow.");
    }
  };
 
  Index a;                 // 'a' in tmp_eq_map[a][j]
  Index b;                 // 'b' in tmp_eq_map[b][j]
  std::set<Index> values;  // equivalence_map for b, relative to a
 
  bool operator<(const EqMapRow &other) const {
    return this->values < other.values;
  }
};
 
struct RelEqMap {
  RelEqMap(Index _a, const std::vector<std::vector<Index>> &tmp_eq_map,
           const SymGroup &g)
      : a(_a) {
    try {
      for (Index b = 0; b < tmp_eq_map.size(); ++b) {
        map.emplace(a, b, tmp_eq_map, g);
      }
    } catch (const std::exception &e) {
      throw libcasm_runtime_error(
          "Error in Orbit constructor: Failed constructing RelEqMap.");
    }
  }
 
  Index a;                 // 'a' in tmp_eq_map[a]
  std::set<EqMapRow> map;  // equivalence_map for b, relative to a
 
  bool operator<(const RelEqMap &other) const { return this->map < other.map; }
};
 
std::ostream &operator<<(std::ostream &sout, const RelEqMap &map);
 
}  // namespace Orbit_impl
 
/* -- Orbit Definitions ------------------------------------- */
 
template <typename _SymCompareType>
Orbit<_SymCompareType>::Orbit(Element generating_element,
                              const SymGroup &generating_group,
                              const SymCompareType &sym_compare)
    : m_generating_group(generating_group),
      m_sym_compare(sym_compare),
      m_invariants(sym_compare.make_invariants(generating_element)) {
  // element(i) compares equivalent to
  // prototype().copy_apply(equivalence_map[i][j]) for all j
 
  //  X   :   op(0)   op(1)   ...
  // proto:    0       ?       ?   (<-- proto.invariant_subgroup(g,
  // sym_compare)) el(1):  (the element generated by the first op that doesn't
  // result in proto) el(2):  (the element generated by the first op that
  // doesn't result in proto or el(1))
  //  ... :  ( etc. )
 
  // get tmp_elements w/ generating_element & g+sym_compare
  // get tmp_equivalence_map from tmp_elements...
  // use tmp_equivalence_map to generate equivalence_maps using the other
  //   elements as the prototype
  // use the element with the 'lowest' equivalence_map as the prototype
 
  const SymGroup &g = generating_group;
 
  auto compare = [&](const Element &A, const Element &B) {
    return m_sym_compare.compare(A, B);
  };
 
  // generate equivalents using std::map to remove duplicates
  // store ordered symop indices corresponding to transformations for each
  // equivalent
  std::map<Element, std::set<Index>, decltype(compare)> t_equiv(compare);
  try {
    for (Index i = 0; i < g.size(); i++) {
      t_equiv[m_sym_compare.prepare(
                  m_sym_compare.copy_apply(g[i], generating_element))]
          .insert(i);  // not g[i].index()!!;
    }
  } catch (const std::exception &e) {
    throw libcasm_runtime_error(
        "Error in Orbit constructor: Failed generating unique equivalents.");
  }
 
  // generate equivalence_map w/ tmp ordering
  std::vector<Element> tmp_element;
  std::vector<std::vector<Index>> tmp_equivalence_map;
  try {
    Index first2proto = g.ind_inverse(*(t_equiv.begin()->second.begin()));
 
    for (auto const &_equiv : t_equiv) {
      tmp_element.push_back(_equiv.first);
      tmp_equivalence_map.emplace_back();
      for (Index op_i : _equiv.second) {
        tmp_equivalence_map.back().push_back(g.ind_prod(op_i, first2proto));
      }
    }
 
    // sanity check that equivalence_map is rectangular
    for (Index j = 1; j < tmp_equivalence_map.size(); ++j) {
      if (tmp_equivalence_map[0].size() != tmp_equivalence_map[j].size()) {
        throw libcasm_runtime_error(
            "Error in Orbit constructor: equivalence map is not rectangular");
      }
    }
  } catch (const std::exception &e) {
    throw libcasm_runtime_error(
        std::string("Error in Orbit constructor: ") +
        "Failed generating initial equivalence map: " + e.what());
  }
 
  Index best_a = 0;
  Orbit_impl::RelEqMap best(best_a, tmp_equivalence_map, g);
  try {
    for (Index a = 1; a < tmp_equivalence_map.size(); ++a) {
      Orbit_impl::RelEqMap test(a, tmp_equivalence_map, g);
 
      if (test < best) {
        best_a = a;
        best = std::move(test);
      }
    }
  } catch (const std::exception &e) {
    throw libcasm_runtime_error(
        "Error in Orbit constructor: Failed generating sorted equivalence "
        "map.");
  }
 
  try {
    Index newproto_i = best.map.begin()->b;
 
    // Loop over equivalents 'i'
    for (const auto &row : best.map) {
      // Index of group element that maps proto to equiv 'i'
      Index proto2i = *(row.values.begin());
 
      Element tequiv = m_sym_compare.prepare(
          m_sym_compare.copy_apply(g[proto2i], tmp_element[newproto_i]));
      SymOp ttrans = m_sym_compare.spatial_transform();
 
      m_element.push_back(m_sym_compare.copy_apply(ttrans * g[proto2i],
                                                   tmp_element[newproto_i]));
      m_equivalence_map.emplace_back();
      // Loop over elements 'j' of new, sorted clust_group (0 row of best.map)
      for (const auto &cg_j : best.map.begin()->values) {
        Index value = g.ind_prod(proto2i, cg_j);
 
        tequiv = m_sym_compare.prepare(
            m_sym_compare.copy_apply(g[value], tmp_element[newproto_i]));
        ttrans = m_sym_compare.spatial_transform();
 
        m_equivalence_map.back().push_back(ttrans * g[value]);
      }
    }
  } catch (const std::exception &e) {
    throw libcasm_runtime_error(
        "Error in Orbit constructor: Failed copying sorted elements and "
        "equivalence map.");
  }
 
  if (m_equivalence_map[0][0].index() != 0) {
    throw libcasm_runtime_error(
        "Error in Orbit constructor: First equivalence map element is not "
        "identity.");
  }
}
 
template <typename _SymCompareType>
Orbit<_SymCompareType> &Orbit<_SymCompareType>::apply_sym(const SymOp &op) {
  // transform elements
  for (auto it = m_element.begin(); it != m_element.end(); ++it) {
    *it = m_sym_compare.copy_apply(op, *it);
  }
 
  // transform equivalence map: std::vector<std::vector<SymOp> >
  for (auto it = m_equivalence_map.begin(); it != m_equivalence_map.end();
       ++it) {
    for (auto op_it = it->begin(); op_it != it->end(); ++op_it) {
      op_it->apply_sym(op);
    }
  }
 
  // transform sym_compare functor
  m_sym_compare.apply_sym(op);
 
  return *this;
}
 
template <typename _SymCompareType>
bool Orbit<_SymCompareType>::operator<(const Orbit &B) const {
  return m_sym_compare.inter_orbit_compare(prototype(), invariants(),
                                           B.prototype(), B.invariants());
}
 
template <typename _SymCompareType>
void Orbit<_SymCompareType>::_construct_canonization_rep() const {
  if (equivalence_map().size() == 0)
    throw libcasm_runtime_error(
        "In Orbit::_construct_canonization_rep(), equivalence_map is "
        "uninitialized or empty! Cannot continue.");
 
  if (size() == 0) {
    m_canonization_rep_ID = SymGroupRepID::identity(0);
    return;
  }
 
  m_canonization_rep_ID =
      equivalence_map()[0][0].master_group().allocate_representation();
 
  for (Index j = 0; j < equivalence_map()[0].size(); j++) {
    std::unique_ptr<SymOpRepresentation> new_rep =
        m_sym_compare
            .canonical_transform(
                m_sym_compare.copy_apply(equivalence_map()[0][j], prototype()))
            ->inverse();
 
    for (Index i = 0; i < equivalence_map().size(); i++) {
      equivalence_map()[i][j].set_rep(m_canonization_rep_ID, *new_rep);
    }
  }
  return;
}
 
template <typename OrbitIterator, typename Element>
OrbitIterator find_orbit(OrbitIterator begin, OrbitIterator end, Element e) {
  typedef typename std::iterator_traits<OrbitIterator>::value_type orbit_type;
  typedef typename orbit_type::InvariantsType InvariantsType;
  const auto &sym_compare = begin->sym_compare();
  auto e_invariants = sym_compare.make_invariants(e);
 
  // first find range of possible orbit by checking invariants
  auto compare = [&](const InvariantsType &A, const InvariantsType &B) {
    return sym_compare.invariants_compare(A, B);
  };
  auto _range =
      std::equal_range(invariants_iterator(begin), invariants_iterator(end),
                       e_invariants, compare);
 
  // find if any of the orbits in range [_range.first, _range.second) contain
  // equivalent
  auto contains = [&](const orbit_type &orbit) { return orbit.contains(e); };
  auto res = std::find_if(_range.first.base(), _range.second.base(), contains);
  if (res == _range.second.base()) {
    return end;
  }
  return res;
}
 
}  // namespace CASM
 
#endif