permutation.h

Go to the documentation of this file.

/*
 *  This file is a part of TiledArray.
 *  Copyright (C) 2013  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/>.
 *
 */
 
#ifndef TILEDARRAY_PERMUTATION_H__INCLUDED
#define TILEDARRAY_PERMUTATION_H__INCLUDED
 
#include <algorithm>
#include <array>
#include <numeric>
 
#include <TiledArray/error.h>
#include <TiledArray/tensor/type_traits.h>
#include <TiledArray/type_traits.h>
#include <TiledArray/util/vector.h>
#include <TiledArray/utility.h>
 
namespace TiledArray {
 
// Forward declarations
class Permutation;
bool operator==(const Permutation&, const Permutation&);
std::ostream& operator<<(std::ostream&, const Permutation&);
template <typename T, std::size_t N>
inline std::array<T, N> operator*(const Permutation&, const std::array<T, N>&);
template <typename T, std::size_t N>
inline std::array<T, N>& operator*=(std::array<T, N>&, const Permutation&);
template <typename T, typename A>
inline std::vector<T> operator*(const Permutation&, const std::vector<T, A>&);
template <typename T, typename A>
inline std::vector<T, A>& operator*=(std::vector<T, A>&, const Permutation&);
template <typename T>
inline std::vector<T> operator*(const Permutation&,
                                const T* MADNESS_RESTRICT const);
 
class BipartitePermutation;
bool operator==(const BipartitePermutation&, const BipartitePermutation&);
 
namespace detail {
 
 
template <typename Perm, typename Arg, typename Result,
          typename = std::enable_if_t<is_permutation_v<Perm>>>
inline void permute_array(const Perm& perm, const Arg& arg, Result& result) {
  using std::size;
  TA_ASSERT(size(result) == size(arg));
  const unsigned int n = size(arg);
  for (unsigned int i = 0u; i < n; ++i) {
    const typename Perm::index_type pi = perm[i];
    TA_ASSERT(i < size(arg));
    TA_ASSERT(pi < size(result));
    result[pi] = arg[i];
  }
}
}  // namespace detail
 
 
class Permutation {
 public:
  typedef Permutation Permutation_;
  typedef unsigned int index_type;
  template <typename T>
  using vector = container::svector<T>;
  typedef vector<index_type>::const_iterator const_iterator;
 
 private:
  template <typename InIter>
  bool valid_permutation(InIter first, InIter last) {
    bool result = true;
    using diff_type = typename std::iterator_traits<InIter>::difference_type;
    const diff_type n = std::distance(first, last);
    TA_ASSERT(n >= 0);
    for (; first != last; ++first) {
      const diff_type value = *first;
      result = result && value >= 0 && (value < n) &&
               (std::count(first, last, *first) == 1ul);
    }
    return result;
  }
 
  // Used to select the correct constructor based on input template types
  struct Enabler {};
 
 protected:
  vector<index_type> p_;
 
 public:
  Permutation() = default;  // constructs an invalid 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<detail::is_input_iterator<
                                 InIter>::value>::type* = nullptr>
  Permutation(InIter first, InIter last) : p_(first, last) {
    TA_ASSERT(valid_permutation(first, last));
  }
 
 
  template <typename Index,
            typename = std::enable_if_t<detail::is_integral_range_v<Index>>>
  explicit Permutation(const Index& a)
      : Permutation(std::cbegin(a), std::cend(a)) {}
 
 
  explicit Permutation(vector<index_type>&& a) : p_(std::move(a)) {
    TA_ASSERT(valid_permutation(p_.begin(), p_.end()));
  }
 
 
  template <typename Integer,
            std::enable_if_t<std::is_integral_v<Integer>>* = nullptr>
  explicit Permutation(std::initializer_list<Integer> list)
      : Permutation(list.begin(), list.end()) {}
 
 
  index_type size() const { return p_.size(); }
 
 
  [[deprecated("use Permutation::size()")]] index_type dim() 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[](unsigned int i) const { return p_[i]; }
 
 
  vector<vector<index_type>> cycles() const {
    vector<vector<index_type>> result;
 
    vector<bool> placed_in_cycle(p_.size(), false);
 
    // 1. for each i compute its orbit
    // 2. if the orbit is longer than 1, sort and add to the list of cycles
    for (index_type i = 0; i != p_.size(); ++i) {
      if (not placed_in_cycle[i]) {
        vector<index_type> cycle(1, i);
        placed_in_cycle[i] = true;
 
        index_type next_i = p_[i];
        while (next_i != i) {
          cycle.push_back(next_i);
          placed_in_cycle[next_i] = true;
          next_i = p_[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;
  }
 
 
  static Permutation identity(const unsigned int dim) {
    Permutation result;
    result.p_.reserve(dim);
    for (unsigned int i = 0u; i < dim; ++i) result.p_.emplace_back(i);
    return result;
  }
 
 
  Permutation identity() const { return identity(p_.size()); }
 
 
  Permutation mult(const Permutation& other) const {
    const unsigned int n = p_.size();
    TA_ASSERT(n == other.p_.size());
    Permutation result;
    result.p_.reserve(n);
 
    for (unsigned int i = 0u; i < n; ++i) {
      const index_type p_i = p_[i];
      const index_type result_i = other.p_[p_i];
      result.p_.emplace_back(result_i);
    }
 
    return result;
  }
 
 
  Permutation inv() const {
    const index_type n = p_.size();
    Permutation result;
    result.p_.resize(n, 0ul);
    for (index_type i = 0ul; i < n; ++i) {
      const index_type pi = p_[i];
      result.p_[pi] = i;
    }
    return 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 = identity(p_.size());
 
    // Compute the power of value with the exponentiation by squaring.
    while (power) {
      if (power & 1) result = result.mult(value);
      value = value.mult(value);
      power >>= 1;
    }
 
    return result;
  }
 
 
  explicit operator bool() const { return !p_.empty(); }
 
 
  bool operator!() const { return p_.empty(); }
 
 
  const auto& data() const { return p_; }
 
 
  template <typename Archive>
  void serialize(Archive& ar) {
    ar& p_;
  }
 
};  // class Permutation
 
 
inline bool operator==(const Permutation& p1, const Permutation& p2) {
  return (p1.size() == p2.size()) &&
         std::equal(p1.data().begin(), p1.data().end(), p2.data().begin());
}
 
 
inline bool operator!=(const Permutation& p1, const Permutation& p2) {
  return !operator==(p1, p2);
}
 
 
inline bool operator<(const Permutation& p1, const Permutation& p2) {
  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) {
  std::size_t n = p.size();
  output << "{";
  for (unsigned int dim = 0; dim < n - 1; ++dim)
    output << dim << "->" << p.data()[dim] << ", ";
  output << n - 1 << "->" << p.data()[n - 1] << "}";
  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) {
  TA_ASSERT(perm.size() == a.size());
  std::array<T, N> result;
  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) {
  TA_ASSERT(perm.size() == a.size());
  const std::array<T, N> temp = a;
  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) {
  TA_ASSERT(perm.size() == v.size());
  std::vector<T> result(perm.size());
  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;
  detail::permute_array(perm, temp, v);
  return v;
}
 
 
template <typename T, std::size_t N>
inline boost::container::small_vector<T, N> operator*(
    const Permutation& perm, const boost::container::small_vector<T, N>& v) {
  TA_ASSERT(perm.size() == v.size());
  boost::container::small_vector<T, N> result(perm.size());
  detail::permute_array(perm, v, result);
  return result;
}
 
 
template <typename T, std::size_t N>
inline boost::container::small_vector<T, N>& operator*=(
    boost::container::small_vector<T, N>& v, const Permutation& perm) {
  const boost::container::small_vector<T, N> temp = v;
  detail::permute_array(perm, temp, v);
  return v;
}
 
 
template <typename T>
inline std::vector<T> operator*(const Permutation& perm,
                                const T* MADNESS_RESTRICT const ptr) {
  const unsigned int n = perm.size();
  std::vector<T> result(n);
  for (unsigned int i = 0u; i < n; ++i) {
    const typename Permutation::index_type perm_i = perm[i];
    const T ptr_i = ptr[i];
    result[perm_i] = ptr_i;
  }
  return result;
}
 
 
class BipartitePermutation {
 public:
  using index_type = Permutation::index_type;
  template <typename T>
  using vector = Permutation::vector<T>;
 
  BipartitePermutation() = default;
  BipartitePermutation(const BipartitePermutation&) = default;
  BipartitePermutation(BipartitePermutation&&) = default;
  ~BipartitePermutation() = default;
  BipartitePermutation& operator=(const BipartitePermutation&) = default;
  BipartitePermutation& operator=(BipartitePermutation&& other) = default;
 
 
  explicit operator bool() const { return static_cast<bool>(base_); }
 
 
  auto size() const { return base_.size(); }
 
 
  [[deprecated("use BipartitePermutation::size()")]] auto dim() const {
    return base_.size();
  }
 
 
  auto begin() const { return base_.begin(); }
 
 
  auto cbegin() const { return base_.cbegin(); }
 
 
  auto end() const { return base_.end(); }
 
 
  auto cend() const { return base_.cend(); }
 
  operator Permutation&() & {
    TA_ASSERT(second_size_ == 0);
    return base_;
  }
  operator Permutation&&() && {
    TA_ASSERT(second_size_ == 0);
    return std::move(base_);
  }
  operator const Permutation&() const& {
    TA_ASSERT(second_size_ == 0);
    return base_;
  }
 
  BipartitePermutation(const Permutation& p,
                       index_type second_partition_size = 0)
      : base_(p), second_size_(second_partition_size) {
    init();
  }
 
  BipartitePermutation(const Permutation& first, const Permutation& second)
      : second_size_(second.size()) {
    vector<index_type> base;
    base.reserve(first.size() + second.size());
    for (auto&& v : first) base.emplace_back(v);
    for (auto&& v : second) base.emplace_back(v);
    base_ = Permutation(base);
    init();
  }
 
  // clang-format off
 
  // clang-format on
  template <typename InIter, typename std::enable_if<detail::is_input_iterator<
                                 InIter>::value>::type* = nullptr>
  BipartitePermutation(InIter first, InIter last,
                       index_type second_partition_size = 0)
      : base_(first, last), second_size_(second_partition_size) {
    init();
  }
 
  // clang-format off
 
  // clang-format on
  template <typename Index,
            typename = std::enable_if_t<detail::is_integral_range_v<Index>>>
  explicit BipartitePermutation(const Index& a,
                                index_type second_partition_size = 0)
      : BipartitePermutation(std::cbegin(a), std::cend(a),
                             second_partition_size) {}
 
  // clang-format off
 
  // clang-format on
  explicit BipartitePermutation(vector<index_type>&& a,
                                index_type second_partition_size = 0)
      : base_(std::move(a)), second_size_(second_partition_size) {
    init();
  }
 
  const Permutation& first() const { return first_; }
  const Permutation& second() const { return second_; }
 
  index_type first_size() const { return this->size() - second_size_; }
 
  index_type second_size() const { return second_size_; }
 
 
  template <typename Archive>
  void serialize(Archive& ar) {
    ar& base_& second_size_;
    if constexpr (madness::archive::is_input_archive<Archive>::value) {
      first_ = {};
      second_ = {};
    }
  }
 
 private:
  Permutation base_;
  index_type second_size_ = 0;
 
  Permutation first_;
  Permutation second_;
 
  // initializes first_ and second_
  void init() {
    first_ = Permutation{this->begin(), this->begin() + first_size()};
    const auto n_first = first_size();
    vector<index_type> temp(second_size());
    for (auto i = n_first; i < size(); ++i)
      temp[i - n_first] = base_[i] - n_first;
    second_ = Permutation(temp.begin(), temp.end());
  }
 
  friend bool operator==(const BipartitePermutation& p1,
                         const BipartitePermutation& p2);
};
 
 
inline bool operator==(const BipartitePermutation& p1,
                       const BipartitePermutation& p2) {
  return (p1.second_size() == p2.second_size()) && p1.base_ == p2.base_;
}
 
 
inline bool operator!=(const BipartitePermutation& p1,
                       const BipartitePermutation& p2) {
  return !(p1 == p2);
}
 
 
inline auto inner(const Permutation& p) {
  abort();
  return Permutation{};
}
 
// N.B. can't return ref here due to possible dangling ref when p is bound to
// temporary
inline auto outer(const Permutation& p) { return p; }
 
inline auto inner_size(const Permutation& p) {
  abort();
  return 0;
}
 
inline auto outer_size(const Permutation& p) { return p.size(); }
 
inline auto inner(const BipartitePermutation& p) { return p.second(); }
 
inline auto outer(const BipartitePermutation& p) { return p.first(); }
 
inline auto inner_size(const BipartitePermutation& p) {
  return p.second_size();
}
 
inline auto outer_size(const BipartitePermutation& p) { return p.first_size(); }
 
}  // namespace TiledArray
 
#endif  // TILEDARRAY_PERMUTATION_H__INCLUED