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 <numeric>
#include <algorithm>
#include <cassert>

#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