permutation_group.h

Go to the documentation of this file.

/*
 *  This file is a part of TiledArray.
 *  Copyright (C) 2015  Virginia Tech
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *  Edward F Valeev, Justus Calvin
 *  Department of Chemistry, Virginia Tech
 *
 *  permutation_group.h
 *  May 13, 2015
 *
 */
 
#ifndef TILEDARRAY_SYMM_PERMUTATION_GROUP_H__INCLUDED
#define TILEDARRAY_SYMM_PERMUTATION_GROUP_H__INCLUDED
 
#include <algorithm>
#include <cassert>
#include <numeric>
 
#include <TiledArray/symm/permutation.h>
 
namespace TiledArray {
 
namespace symmetry {
 
 
class PermutationGroup {
 public:
  using Permutation = TiledArray::symmetry::Permutation;
  using element_type = Permutation;
 
 protected:
  std::vector<Permutation> generators_;
  std::vector<Permutation> elements_;
 
 public:
  // Compiler generated functions
  PermutationGroup(const PermutationGroup&) = default;
  PermutationGroup(PermutationGroup&&) = default;
  PermutationGroup& operator=(const PermutationGroup&) = default;
  PermutationGroup& operator=(PermutationGroup&&) = default;
 
 
  PermutationGroup(std::vector<Permutation> generators)
      : generators_(std::move(generators)) {
    init(generators_, elements_);
  }
 
 
  unsigned int order() const { return elements_.size(); }
 
 
  static Permutation identity() { return Permutation(); }
 
 
  const Permutation& operator[](unsigned int i) const {
    TA_ASSERT(i < elements_.size());
    return elements_[i];
  }
 
 
  const std::vector<Permutation>& elements() const { return elements_; }
 
 
  const std::vector<Permutation>& generators() const { return generators_; }
 
 
 
 
  std::vector<Permutation>::const_iterator begin() const {
    return elements_.begin();
  }
 
 
  std::vector<Permutation>::const_iterator cbegin() const {
    return elements_.cbegin();
  }
 
 
  std::vector<Permutation>::const_iterator end() const {
    return elements_.end();
  }
 
 
  std::vector<Permutation>::const_iterator cend() const {
    return elements_.cend();
  }
 
 
 
  template <typename Set>
  Set domain() const {
    Set result;
    // sufficient to loop over generators
    for (const auto& e : generators_) {
      const auto e_domain = e.domain<Set>();
      result.insert(e_domain.begin(), e_domain.end());
    }
    return result;
  }
 
 protected:
  PermutationGroup() {}  // makes uninitialized group, all initialization is
                         // left to the derived class
 
  static void init(std::vector<Permutation>& generators,
                   std::vector<Permutation>& elements) {
    sort(generators.begin(), generators.end());
    auto unique_last = std::unique(generators.begin(), generators.end());
    generators.erase(unique_last, generators.end());
    {  // eliminate identity from the generator list
      auto I_iter =
          std::find(generators.begin(), generators.end(), Permutation());
      if (I_iter != generators.end()) generators.erase(I_iter);
    }
 
    // add the identity element first
    elements.emplace_back();
 
    for (const auto& g : generators) {
      elements.push_back(g);
    }
 
    // Generate the remaining elements in the group by multiplying by generators
    for (unsigned int g = 1u; g < elements.size(); ++g) {
      for (const auto& G : generators) {
        Permutation e = elements[g] * G;
        if (std::find(elements.cbegin(), elements.cend(), e) ==
            elements.cend()) {
          elements.emplace_back(std::move(e));
        }
      }
    }
 
    sort(elements.begin(), elements.end());
  }
 
};  // class PermutationGroup
 
 
inline bool operator==(const PermutationGroup& p1, const PermutationGroup& p2) {
  return (p1.order() == p2.order()) && p1.elements() == p2.elements();
}
 
 
inline bool operator!=(const PermutationGroup& p1, const PermutationGroup& p2) {
  return !operator==(p1, p2);
}
 
 
inline bool operator<(const PermutationGroup& p1, const PermutationGroup& p2) {
  return std::lexicographical_compare(p1.cbegin(), p1.cend(), p2.cbegin(),
                                      p2.cend());
}
 
 
inline std::ostream& operator<<(std::ostream& output,
                                const PermutationGroup& p) {
  output << "{";
  for (auto i = p.cbegin(); i != p.cend();) {
    output << *i;
    if (++i != p.cend()) output << ", ";
  }
  output << "}";
  return output;
}
 
 
class SymmetricGroup final : public PermutationGroup {
 public:
  using index_type = Permutation::index_type;
 
  // Compiler generated functions
  SymmetricGroup() = delete;
  SymmetricGroup(const SymmetricGroup&) = default;
  SymmetricGroup(SymmetricGroup&&) = default;
  SymmetricGroup& operator=(const SymmetricGroup&) = default;
  SymmetricGroup& operator=(SymmetricGroup&&) = default;
 
  SymmetricGroup(unsigned int degree)
      : SymmetricGroup(SymmetricGroup::iota_vector(degree)) {}
 
  template <typename InputIterator,
            typename std::enable_if< ::TiledArray::detail::is_input_iterator<
                InputIterator>::value>::type* = nullptr>
  SymmetricGroup(InputIterator begin, InputIterator end)
      : PermutationGroup(), domain_(begin, end) {
    for (auto iter = begin; iter != end; ++iter) {
      TA_ASSERT(*iter >= 0);
    }
 
    const auto degree = domain_.size();
 
    // Add generators to the list of elements
    if (degree > 2u) {
      for (unsigned int i = 0u; i < degree; ++i) {
        // Construct the generator and add to the list
        unsigned int i1 = (i + 1u) % degree;
        generators_.emplace_back(Permutation::Map{{domain_[i], domain_[i1]},
                                                  {domain_[i1], domain_[i]}});
      }
    } else if (degree == 2u) {
      // Construct the generator
      generators_.emplace_back(
          Permutation::Map{{domain_[0], domain_[1]}, {domain_[1], domain_[0]}});
    }
 
    init(generators_, elements_);
  }
 
 
  template <typename Integer,
            typename std::enable_if<std::is_integral<Integer>::value>::type* =
                nullptr>
  explicit SymmetricGroup(std::initializer_list<Integer> list)
      : SymmetricGroup(list.begin(), list.end()) {}
 
 
  unsigned int degree() const { return domain_.size(); }
 
 private:
  std::vector<index_type> domain_;
 
  static std::vector<index_type> iota_vector(size_t n) {
    std::vector<index_type> result(n);
    std::iota(result.begin(), result.end(), 0);
    return result;
  }
 
  SymmetricGroup(const std::vector<index_type>& domain)
      : SymmetricGroup(domain.begin(), domain.end()) {}
};
 
template <typename MultiIndex>
bool is_lexicographically_smallest(const MultiIndex& idx,
                                   const PermutationGroup& pg) {
  const auto idx_size = idx.size();
  for (const auto& p : pg) {
    for (size_t i = 0; i != idx_size; ++i) {
      auto idx_i = idx[i];
      auto idx_p_i = idx[p[i]];
      if (idx_p_i < idx_i) return false;
      if (idx_p_i > idx_i) break;
    }
  }
  return true;
}
 
 
inline PermutationGroup conjugate(const PermutationGroup& G,
                                  const PermutationGroup::Permutation& h) {
  std::vector<PermutationGroup::Permutation> conjugate_generators;
  const auto h_inv = h.inv();
  for (const auto& generator : G.generators()) {
    conjugate_generators.emplace_back(h * generator * h_inv);
  }
  return PermutationGroup{std::move(conjugate_generators)};
}
 
inline PermutationGroup intersect(const PermutationGroup& G1,
                                  const PermutationGroup& G2) {
  std::vector<PermutationGroup::Permutation> intersect_elements;
  std::set_intersection(G1.begin(), G1.end(), G2.begin(), G2.end(),
                        std::back_inserter(intersect_elements));
  return PermutationGroup{std::move(intersect_elements)};
}
 
 
template <typename Set>
inline PermutationGroup stabilizer(const PermutationGroup& G, const Set& f) {
  std::vector<PermutationGroup::Permutation> stabilizer_generators;
  for (const auto& generator : G.generators()) {
    bool fixes_set = true;
    for (const auto& i : f) {
      if (generator.is_in_domain(i)) {
        fixes_set = false;
        break;
      }
    }
    if (fixes_set) stabilizer_generators.push_back(generator);
  }
  return PermutationGroup{std::move(stabilizer_generators)};
}
 
}  // namespace symmetry
}  // namespace TiledArray
 
#endif  // TILEDARRAY_SYMM_PERMUTATION_GROUP_H__INCLUDED