permutation.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 Valeev
 *  Department of Chemistry, Virginia Tech
 *
 *  permutation.h
 *  May 13, 2015
 *
 */

#ifndef TILEDARRAY_SYMM_PERMUTATION_H__INCLUDED
#define TILEDARRAY_SYMM_PERMUTATION_H__INCLUDED

#include <array>
#include <vector>
#include <map>
#include <set>

#include <algorithm>

#include <TiledArray/type_traits.h>
#include <TiledArray/error.h>

namespace TiledArray {

  namespace symmetry {

    namespace detail {


      template <typename Perm, typename Arg, typename Result>
      inline void permute_array(const Perm& perm, const Arg& arg, Result& result) {
        TA_ASSERT(result.size() == arg.size());
        if (perm.domain_size() < arg.size()) {
          std::copy(arg.begin(), arg.end(), result.begin()); // if perm does not map every element of arg, copy arg to result first
        }
        for(const auto& p: perm.data()) {
          TA_ASSERT(result.size() > p.second);
          TA_ASSERT(arg.size() > p.second);
          result[p.second] = arg[p.first];
        }
      }
    } // namespace detail


    class Permutation {
    public:
      typedef Permutation Permutation_;
      typedef unsigned int index_type;
      typedef std::map<index_type,index_type> Map;
      typedef Map::const_iterator const_iterator;

    private:

      Map p_;

      static std::ostream& print_map(std::ostream& output, const Map& p) {
        for (auto i=p.cbegin(); i!=p.cend();) {
          output << i->first << "->" << i->second;
          if (++i != p.cend())
            output << ", ";
        }
        return output;
      }
      friend inline std::ostream& operator<<(std::ostream& output, const Permutation& p);

      template <typename InIter>
      static bool valid_permutation(InIter first, InIter last) {
        bool result = true;
        for(; first != last; ++first) {
          const auto value = *first;
          result = result && value >= 0 && (std::count(first, last, *first) == 1ul);
        }
        return result;
      }

      template <typename Map>
      static bool valid_permutation(const Map& input) {
        std::set<index_type> keys;
        std::set<index_type> values;
        for(const auto& e: input) {
          const auto& key = e.first;
          const auto& value = e.second;
          if (keys.find(key) == keys.end())
            keys.insert(key);
          else {
            //print_map(std::cout, input);
            return false; // key is found more than once
          }
          if (values.find(value) == values.end())
            values.insert(value);
          else {
            //print_map(std::cout, input);
            return false; // value is found more than once
          }
        }
        return keys == values; // must map domain into itself
      }

      // Used to select the correct constructor based on input template types
      struct Enabler { };

    public:

      Permutation() = default; // makes an identity permutation
      Permutation(const Permutation&) = default;
      Permutation(Permutation&&) = default;
      ~Permutation() = default;
      Permutation& operator=(const Permutation&) = default;
      Permutation& operator=(Permutation&& other) = default;


      template <typename InIter,
          typename std::enable_if<TiledArray::detail::is_input_iterator<InIter>::value>::type* = nullptr>
      Permutation(InIter first, InIter last)
      {
        TA_ASSERT( valid_permutation(first, last) );
        size_t i=0;
        for(auto e=first; e!=last; ++e, ++i) {
          auto p_i = *e;
          if (i != p_i) p_[i] = p_i;
        }
      }


      explicit Permutation(const std::vector<index_type>& a) :
          Permutation(a.begin(), a.end())
      {
      }


      explicit Permutation(std::initializer_list<index_type> list) :
          Permutation(list.begin(), list.end())
      {
      }


      Permutation(Map p) : p_(std::move(p))
      {
        TA_ASSERT(valid_permutation(p_));
      }


      unsigned int domain_size() const { return p_.size(); }



      const_iterator begin() const { return p_.begin(); }


      const_iterator cbegin() const { return p_.cbegin(); }


      const_iterator end() const { return p_.end(); }


      const_iterator cend() const { return p_.cend(); }



      index_type operator[](index_type e) const {
        auto e_image_iter = p_.find(e);
        if (e_image_iter != p_.end())
          return e_image_iter->second;
        else
          return e;
      }


      template <typename Set>
      Set domain() const {
        Set result;
        for(const auto& e: this->p_) {
          result.insert(result.begin(), e.first);
        }
        //std::sort(result.begin(), result.end());
        return result;
      }


      template <typename Integer,
                typename std::enable_if<std::is_integral<Integer>::value>::type* = nullptr>
      bool is_in_domain(Integer i) const {
        return p_.find(i) != p_.end();
      }


      std::vector<std::vector<index_type> > cycles() const {

        std::vector<std::vector<index_type>> result;

        std::set<index_type> placed_in_cycle;

        // safe to use non-const operator[] due to properties of Permutation (domain==image)
        auto& p_nonconst_ = const_cast<Map&>(p_);

        // 1. for each i compute its orbit
        // 2. if the orbit is longer than 1, sort and add to the list of cycles
        for(const auto& e: p_) {
          auto i = e.first;
          if (placed_in_cycle.find(i) == placed_in_cycle.end()) { // not in a cycle yet?
            std::vector<index_type> cycle(1,i);
            placed_in_cycle.insert(i);

            index_type next_i = p_nonconst_[i];
            while (next_i != i) {
              cycle.push_back(next_i);
              placed_in_cycle.insert(next_i);
              next_i = p_nonconst_[next_i];
            }

            if (cycle.size() != 1) {
              std::sort(cycle.begin(), cycle.end());
              result.emplace_back(cycle);
            }

          } // this i already in a cycle
        } // loop over i

        return result;
      }


      Permutation mult(const Permutation& other) const {
        // 1. domain of product = domain of this U domain of other
        using iset = std::set<index_type>;
        auto product_domain = this->domain<iset>();
        auto other_domain = other.domain<iset>();
        product_domain.insert(other_domain.begin(), other_domain.end());

        Map result;
        for(const auto& d: product_domain) {
          const auto d_image = other[(*this)[d]];
          if (d_image != d)
            result[d] = d_image;
        }

        return Permutation(result);
      }


      Permutation inv() const {
        Map result;
        for(const auto& iter: p_) {
          const auto i = iter.first;
          const auto pi = iter.second;
          result.insert(std::make_pair(pi,i));
        }
        return Permutation(std::move(result));
      }


      Permutation pow(int n) const {
        // Initialize the algorithm inputs
        int power;
        Permutation value;
        if(n < 0) {
          value = inv();
          power = -n;
        } else {
          value = *this;
          power = n;
        }

        Permutation result;

        while(power) {
            if(power & 1)
              result = result.mult(value);
            value = value.mult(value);
            power >>= 1;
        }

        return result;
      }


      const Map& data() const { return p_; }


      Permutation identity() const {
        return Permutation();
      }


      template <typename Archive>
      void serialize(Archive& ar) {
        ar & p_;
      }

    }; // class Permutation


    inline bool operator==(const Permutation& p1, const Permutation& p2) {
      return (p1.domain_size() == p2.domain_size())
             && p1.data() == p2.data();
    }


    inline bool operator!=(const Permutation& p1, const Permutation& p2) {
      return ! operator==(p1, p2);
    }


    inline bool operator<(const Permutation& p1, const Permutation& p2) {
      if (&p1 == &p2) return false;
      return std::lexicographical_compare(p1.data().begin(), p1.data().end(),
          p2.data().begin(), p2.data().end());
    }


    inline std::ostream& operator<<(std::ostream& output, const Permutation& p) {
      output << "{";
      Permutation::print_map(output, p.data());
      output << "}";
      return output;
    }


    inline Permutation operator -(const Permutation& perm) {
      return perm.inv();
    }


    inline Permutation operator*(const Permutation& p1, const Permutation& p2) {
      return p1.mult(p2);
    }


    inline Permutation& operator*=(Permutation& p1, const Permutation& p2) {
      return (p1 = p1 * p2);
    }


    inline Permutation operator^(const Permutation& perm, int n) {
      return perm.pow(n);
    }


    template <typename T, std::size_t N>
    inline std::array<T,N> operator*(const Permutation& perm, const std::array<T, N>& a) {
      std::array<T,N> result;
      symmetry::detail::permute_array(perm, a, result);
      return result;
    }


    template <typename T, std::size_t N>
    inline std::array<T,N>& operator*=(std::array<T,N>& a, const Permutation& perm) {
      const std::array<T,N> temp = a;
      symmetry::detail::permute_array(perm, temp, a);
      return a;
    }


    template <typename T, typename A>
    inline std::vector<T> operator*(const Permutation& perm, const std::vector<T, A>& v) {
      std::vector<T> result(v.size());
      symmetry::detail::permute_array(perm, v, result);
      return result;
    }


    template <typename T, typename A>
    inline std::vector<T, A>& operator*=(std::vector<T, A>& v, const Permutation& perm) {
      const std::vector<T, A> temp = v;
      symmetry::detail::permute_array(perm, temp, v);
      return v;
    }

  } // namespace symmetry

} // namespace TiledArray

#endif // TILEDARRAY_SYMM_PERMUTATION_H__INCLUDED