Permutations.h

Go to the documentation of this file.

#ifndef INCLUDE_MOLASSEMBLER_TEMPLE_PERMUTATIONS_H
#define INCLUDE_MOLASSEMBLER_TEMPLE_PERMUTATIONS_H

#include <algorithm>
#include <stdexcept>

namespace Scine {
namespace Molassembler {
namespace Temple {

template<class Container>
constexpr std::size_t permutationIndex(const Container& container) {
  const std::size_t size = container.size();

  std::size_t index = 0;
  std::size_t position = 2;// position 1 is paired with factor 0 and so is skipped
  std::size_t factor = 1;

  for(std::size_t p = size - 2; p != std::numeric_limits<std::size_t>::max(); --p) {
    std::size_t largerSuccessors = 0;

    for(std::size_t q = p + 1; q < size; ++q) {
      if(container[q] < container[p]) {
        ++largerSuccessors;
      }
    }

    index += (largerSuccessors * factor);
    factor *= position;
    ++position;
  }

  return index;
}

template<typename Container>
constexpr void inPlaceSwap(
  Container& data,
  const std::size_t a,
  const std::size_t b
) {
  auto intermediate = std::move(data.at(b));
  data.at(b) = std::move(data.at(a));
  data.at(a) = std::move(intermediate);
}

template<typename Container>
constexpr void inPlaceReverse(
  Container& data,
  const std::size_t indexFrom,
  const std::size_t indexTo
) {
  std::size_t a = indexFrom;
  std::size_t b = indexTo;
  while(a != b && a != --b) {
    inPlaceSwap(data, a++, b);
  }
}

template<typename Container>
constexpr bool inPlaceNextPermutation(
  Container& data,
  const std::size_t first,
  const std::size_t last
) {
  std::size_t size = data.size();
  if(!(
    first < last
    && first < size
    && last <= size
  )) {
    throw "Call parameters to inPlaceNextPermutation make no sense!";
  }

  std::size_t i = last - 1;
  std::size_t j = 0;
  std::size_t k = 0;

  while(true) {
    j = i;

    if(
      i != 0
      && data.at(--i) < data.at(j)
    ) {
      k = last;

      while(
        k != 0
        && !(data.at(i) < data.at(--k))
      ) {}

      inPlaceSwap(data, i, k);
      inPlaceReverse(data, j, last);
      return true;
    }

    if(i == first) {
      inPlaceReverse(data, first, last);
      return false;
    }
  }
}

template<typename Container>
constexpr bool inPlaceNextPermutation(Container& data) {
  return inPlaceNextPermutation(data, 0, data.size());
}

template<typename Container>
constexpr bool inPlacePreviousPermutation(
  Container& data,
  const std::size_t first,
  const std::size_t last
) {
  unsigned size = data.size();
  if(!(
    first < last
    && first < size
    && last <= size
  )) {
    throw "Call parameters to inPlaceNextPermutation make no sense!";
  }

  std::size_t i = last - 1;
  std::size_t j = 0;
  std::size_t k = 0;

  while(true) {
    j = i;

    if(
      i != 0
      && data.at(j) < data.at(--i)
    ) {
      k = last;

      while(
        k != 0
        && !(data.at(--k) < data.at(i))
      ) {}

      inPlaceSwap(data, i, k);
      inPlaceReverse(data, j, last);
      return true;
    }

    if(i == first) {
      inPlaceReverse(data, first, last);
      return false;
    }
  }
}

template<typename Container>
constexpr bool inPlacePreviousPermutation(Container& data) {
  return inPlacePreviousPermutation(data, 0, data.size());
}


template<class Container>
bool next_permutation(Container& container) {
  return std::next_permutation(
    std::begin(container),
    std::end(container)
  );
}

template<class Container>
bool prev_permutation(Container& container) {
  return std::prev_permutation(
    std::begin(container),
    std::end(container)
  );
}

template<class Container>
bool nextCombinationPermutation(
  Container& toPermute,
  const Container& limits
) {
  assert(toPermute.size() == limits.size());
  const unsigned cols = toPermute.size();

  // Check if all columns are full
  bool allFull = true;
  for(unsigned i = 0; i < cols; ++i) {
    if(toPermute[i] != limits[i]) {
      allFull = false;
      break;
    }
  }

  if(allFull) {
    return false;
  }

  // Make next permutation
  for(int i = cols - 1; i >= 0; --i) {
    if(toPermute[i] == limits[i]) {
      toPermute[i] = 0;
    } else {
      ++toPermute[i];
      return true;
    }
  }

  return true;
}

template<typename Container>
struct Permutation {
  using T = std::decay_t<decltype(std::declval<Container>().at(0))>;
  static_assert(
    std::is_convertible<T, unsigned>::value,
    "Permutation container value type must be convertible to an integral type"
  );

  constexpr Permutation() = default;

  constexpr explicit Permutation(unsigned N) : sigma(N) {
    for(unsigned i = 0; i < N; ++i) {
      sigma.at(i) = i;
    }
  }

  constexpr explicit Permutation(Container p) : sigma(std::move(p)) {}

  constexpr Permutation(const unsigned N, unsigned i) : sigma(N) {
    Container factorials(N);
    factorials.at(0) = 1;
    for(unsigned k = 1; k < N; ++k) {
      factorials.at(k) = factorials.at(k - 1) * k;
    }

    for(unsigned k = 0; k < N; ++k) {
      const unsigned fac = factorials.at(N - 1 - k);
      sigma.at(k) = i / fac;
      i %= fac;
    }

    for(int k = static_cast<int>(N) - 1; k > 0; --k) {
      for(int j = k - 1; j >= 0; --j) {
        if(sigma.at(j) <= sigma.at(k)) {
          ++sigma.at(k);
        }
      }
    }
  }

  template<typename OtherContainer>
  static constexpr Permutation ordering(const OtherContainer& unordered) {
    Permutation p(unordered.size());
    std::sort(
      std::begin(p.sigma),
      std::end(p.sigma),
      [&](const auto i, const auto j) {
        return unordered.at(i) < unordered.at(j);
      }
    );
    return p;
  }

  constexpr bool next() {
    return inPlaceNextPermutation(sigma);
  }

  constexpr bool prev() {
    return inPlacePreviousPermutation(sigma);
  }

  constexpr std::size_t index() const {
    return permutationIndex(sigma);
  }

  constexpr Permutation inverse() const {
    auto inverted = sigma;
    const unsigned size = inverted.size();
    for(unsigned i = 0; i < size; ++i) {
      inverted.at(sigma.at(i)) = i;
    }
    return Permutation(inverted);
  }

  constexpr auto at(const std::size_t i) const -> decltype(auto) {
    return sigma.at(i);
  }

  constexpr auto operator() (const std::size_t i) const -> decltype(auto) {
    return sigma.at(i);
  }

  template<typename OtherContainer>
  constexpr OtherContainer apply(const OtherContainer& other) const {
    const unsigned size = other.size();
    OtherContainer result(size);
    for(std::size_t i = 0; i < size; ++i) {
      result.at(i) = other.at(sigma.at(i));
    }
    return result;
  }

  constexpr Permutation compose(const Permutation& other) const {
    return Permutation {apply(other.sigma)};
  }

  Container sigma;
};

template<typename Container>
Permutation<Container> make_permutation(Container p) {
  return Permutation<Container>(std::move(p));
}

} // namespace Temple
} // namespace Molassembler
} // namespace Scine

#endif