permutation.h
Go to the documentation of this file.
1 /*
2  * This file is a part of TiledArray.
3  * Copyright (C) 2015 Virginia Tech
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Edward Valeev
19  * Department of Chemistry, Virginia Tech
20  *
21  * permutation.h
22  * May 13, 2015
23  *
24  */
25 
26 #ifndef TILEDARRAY_SYMM_PERMUTATION_H__INCLUDED
27 #define TILEDARRAY_SYMM_PERMUTATION_H__INCLUDED
28 
29 #include <array>
30 #include <map>
31 #include <set>
32 #include <vector>
33 
34 #include <algorithm>
35 
36 #include <TiledArray/error.h>
37 #include <TiledArray/type_traits.h>
38 #include <TiledArray/util/vector.h>
39 
40 namespace TiledArray {
41 
47 namespace symmetry {
48 
50 
117 class Permutation {
118  public:
120  typedef unsigned int index_type;
121  template <typename T>
123  typedef std::map<index_type, index_type> Map;
124  typedef Map::const_iterator const_iterator;
125 
126  private:
128  Map p_;
129 
130  static std::ostream& print_map(std::ostream& output, const Map& p) {
131  for (auto i = p.cbegin(); i != p.cend();) {
132  output << i->first << "->" << i->second;
133  if (++i != p.cend()) output << ", ";
134  }
135  return output;
136  }
137  friend inline std::ostream& operator<<(std::ostream& output,
138  const Permutation& p);
139 
143  template <typename InIter>
144  static bool valid_permutation(InIter first, InIter last) {
145  bool result = true;
146  for (; first != last; ++first) {
147  const auto value = *first;
148  result = result && value >= 0 && (std::count(first, last, *first) == 1ul);
149  }
150  return result;
151  }
152 
156  template <typename Map>
157  static bool valid_permutation(const Map& input) {
158  std::set<index_type> keys;
159  std::set<index_type> values;
160  for (const auto& e : input) {
161  const auto& key = e.first;
162  const auto& value = e.second;
163  if (keys.find(key) == keys.end())
164  keys.insert(key);
165  else {
166  // print_map(std::cout, input);
167  return false; // key is found more than once
168  }
169  if (values.find(value) == values.end())
170  values.insert(value);
171  else {
172  // print_map(std::cout, input);
173  return false; // value is found more than once
174  }
175  }
176  return keys == values; // must map domain into itself
177  }
178 
179  // Used to select the correct constructor based on input template types
180  struct Enabler {};
181 
182  public:
183  Permutation() = default; // makes an identity permutation
184  Permutation(const Permutation&) = default;
185  Permutation(Permutation&&) = default;
186  ~Permutation() = default;
187  Permutation& operator=(const Permutation&) = default;
188  Permutation& operator=(Permutation&& other) = default;
189 
191 
197  template <typename InIter,
198  typename std::enable_if<TiledArray::detail::is_input_iterator<
199  InIter>::value>::type* = nullptr>
200  Permutation(InIter first, InIter last) {
201  TA_ASSERT(valid_permutation(first, last));
202  index_type i = 0;
203  for (auto e = first; e != last; ++e, ++i) {
204  index_type p_i = *e;
205  if (i != p_i) p_[i] = p_i;
206  }
207  }
208 
210 
213  template <typename Index, typename std::enable_if<
214  TiledArray::detail::is_integral_range_v<Index>,
215  bool>::type* = nullptr>
216  explicit Permutation(Index&& a) : Permutation(begin(a), end(a)) {}
217 
219 
222  template <typename Integer,
223  std::enable_if_t<std::is_integral_v<Integer>>* = nullptr>
224  explicit Permutation(std::initializer_list<Integer> list)
225  : Permutation(list.begin(), list.end()) {}
226 
228 
230  Permutation(Map p) : p_(std::move(p)) { TA_ASSERT(valid_permutation(p_)); }
231 
233 
235  unsigned int domain_size() const { return p_.size(); }
236 
245 
247 
249  const_iterator begin() const { return p_.begin(); }
250 
252 
254  const_iterator cbegin() const { return p_.cbegin(); }
255 
257 
259  const_iterator end() const { return p_.end(); }
260 
262 
264  const_iterator cend() const { return p_.cend(); }
265 
267 
269 
274  auto e_image_iter = p_.find(e);
275  if (e_image_iter != p_.end())
276  return e_image_iter->second;
277  else
278  return e;
279  }
280 
282 
285  template <typename Set>
286  Set domain() const {
287  Set result;
288  for (const auto& e : this->p_) {
289  result.insert(result.begin(), e.first);
290  }
291  // std::sort(result.begin(), result.end());
292  return result;
293  }
294 
296 
300  template <typename Integer,
301  typename std::enable_if<std::is_integral<Integer>::value>::type* =
302  nullptr>
303  bool is_in_domain(Integer i) const {
304  return p_.find(i) != p_.end();
305  }
306 
308 
322 
323  std::set<index_type> placed_in_cycle;
324 
325  // safe to use non-const operator[] due to properties of Permutation
326  // (domain==image)
327  auto& p_nonconst_ = const_cast<Map&>(p_);
328 
329  // 1. for each i compute its orbit
330  // 2. if the orbit is longer than 1, sort and add to the list of cycles
331  for (const auto& e : p_) {
332  auto i = e.first;
333  if (placed_in_cycle.find(i) ==
334  placed_in_cycle.end()) { // not in a cycle yet?
335  vector<index_type> cycle(1, i);
336  placed_in_cycle.insert(i);
337 
338  index_type next_i = p_nonconst_[i];
339  while (next_i != i) {
340  cycle.push_back(next_i);
341  placed_in_cycle.insert(next_i);
342  next_i = p_nonconst_[next_i];
343  }
344 
345  if (cycle.size() != 1) {
346  std::sort(cycle.begin(), cycle.end());
347  result.emplace_back(cycle);
348  }
349 
350  } // this i already in a cycle
351  } // loop over i
352 
353  return result;
354  }
355 
357 
360  Permutation mult(const Permutation& other) const {
361  // 1. domain of product = domain of this U domain of other
362  using iset = std::set<index_type>;
363  auto product_domain = this->domain<iset>();
364  auto other_domain = other.domain<iset>();
365  product_domain.insert(other_domain.begin(), other_domain.end());
366 
367  Map result;
368  for (const auto& d : product_domain) {
369  const auto d_image = other[(*this)[d]];
370  if (d_image != d) result[d] = d_image;
371  }
372 
373  return Permutation(result);
374  }
375 
377 
381  Permutation inv() const {
382  Map result;
383  for (const auto& iter : p_) {
384  const auto i = iter.first;
385  const auto pi = iter.second;
386  result.insert(std::make_pair(pi, i));
387  }
388  return Permutation(std::move(result));
389  }
390 
392 
397  Permutation pow(int n) const {
398  // Initialize the algorithm inputs
399  int power;
400  Permutation value;
401  if (n < 0) {
402  value = inv();
403  power = -n;
404  } else {
405  value = *this;
406  power = n;
407  }
408 
409  Permutation result;
410 
411  while (power) {
412  if (power & 1) result = result.mult(value);
413  value = value.mult(value);
414  power >>= 1;
415  }
416 
417  return result;
418  }
419 
421 
425  const Map& data() const { return p_; }
426 
428 
430  Permutation identity() const { return Permutation(); }
431 
433 
437  template <typename Archive>
438  void serialize(Archive& ar) {
439  ar& p_;
440  }
441 
442 }; // class Permutation
443 
445 
450 inline bool operator==(const Permutation& p1, const Permutation& p2) {
451  return (p1.domain_size() == p2.domain_size()) && p1.data() == p2.data();
452 }
453 
455 
460 inline bool operator!=(const Permutation& p1, const Permutation& p2) {
461  return !operator==(p1, p2);
462 }
463 
465 
470 inline bool operator<(const Permutation& p1, const Permutation& p2) {
471  if (&p1 == &p2) return false;
472  return std::lexicographical_compare(p1.data().begin(), p1.data().end(),
473  p2.data().begin(), p2.data().end());
474 }
475 
477 
481 inline std::ostream& operator<<(std::ostream& output, const Permutation& p) {
482  output << "{";
483  Permutation::print_map(output, p.data());
484  output << "}";
485  return output;
486 }
487 
489 
492 inline Permutation operator-(const Permutation& perm) { return perm.inv(); }
493 
495 
500 inline Permutation operator*(const Permutation& p1, const Permutation& p2) {
501  return p1.mult(p2);
502 }
503 
506  return (p1 = p1 * p2);
507 }
508 
510 
516 inline Permutation operator^(const Permutation& perm, int n) {
517  return perm.pow(n);
518 }
519 
522 namespace detail {
523 
524 // clang-format off
526 
534 // clang-format on
535 template <typename Arg, typename Result>
537  const Arg& arg, Result& result) {
538  TA_ASSERT(result.size() == arg.size());
539  if (perm.domain_size() < arg.size()) {
540  std::copy(arg.begin(), arg.end(),
541  result.begin()); // if perm does not map every element of arg,
542  // copy arg to result first
543  }
544  for (const auto& p : perm.data()) {
545  TA_ASSERT(result.size() > p.second);
546  TA_ASSERT(arg.size() > p.second);
547  result[p.second] = arg[p.first];
548  }
549 }
550 } // namespace detail
551 
553 
561 template <typename T, std::size_t N>
562 inline std::array<T, N> operator*(const Permutation& perm,
563  const std::array<T, N>& a) {
564  std::array<T, N> result;
565  symmetry::detail::permute_array(perm, a, result);
566  return result;
567 }
568 
570 
578 template <typename T, std::size_t N>
579 inline std::array<T, N>& operator*=(std::array<T, N>& a,
580  const Permutation& perm) {
581  const std::array<T, N> temp = a;
582  symmetry::detail::permute_array(perm, temp, a);
583  return a;
584 }
585 
587 
595 template <typename T, typename A>
596 inline std::vector<T> operator*(const Permutation& perm,
597  const std::vector<T, A>& v) {
598  std::vector<T> result(v.size());
599  symmetry::detail::permute_array(perm, v, result);
600  return result;
601 }
602 
604 
612 template <typename T, typename A>
613 inline std::vector<T, A>& operator*=(std::vector<T, A>& v,
614  const Permutation& perm) {
615  const std::vector<T, A> temp = v;
616  symmetry::detail::permute_array(perm, temp, v);
617  return v;
618 }
619 
621 
629 template <typename T, std::size_t N>
630 inline boost::container::small_vector<T, N> operator*(
631  const Permutation& perm, const boost::container::small_vector<T, N>& v) {
632  boost::container::small_vector<T, N> result(v.size());
633  symmetry::detail::permute_array(perm, v, result);
634  return result;
635 }
636 
638 
646 template <typename T, std::size_t N>
647 inline boost::container::small_vector<T, N>& operator*=(
648  boost::container::small_vector<T, N>& v, const Permutation& perm) {
649  const boost::container::small_vector<T, N> temp = v;
650  symmetry::detail::permute_array(perm, temp, v);
651  return v;
652 }
653 
654 } // namespace symmetry
655 
658 } // namespace TiledArray
659 
660 #endif // TILEDARRAY_SYMM_PERMUTATION_H__INCLUDED
boost::container::small_vector< T, N > svector
Definition: vector.h:43
Permutation & operator*=(Permutation &p1, const Permutation &p2)
return *this ^ other
Definition: permutation.h:505
bool operator!=(const Permutation &p1, const Permutation &p2)
Permutation inequality operator.
Definition: permutation.h:460
Permutation(const Permutation &)=default
Permutation(std::initializer_list< Integer > list)
Construct permutation with an initializer list.
Definition: permutation.h:224
Permutation operator*(const Permutation &p1, const Permutation &p2)
Permutation multiplication operator.
Definition: permutation.h:500
void serialize(Archive &ar)
Serialize permutation.
Definition: permutation.h:438
std::ostream & operator<<(std::ostream &output, const Permutation &p)
Add permutation to an output stream.
Definition: permutation.h:481
void permute_array(const TiledArray::symmetry::Permutation &perm, const Arg &arg, Result &result)
Create a permuted copy of an array.
Definition: permutation.h:536
bool operator==(const Permutation &p1, const Permutation &p2)
Permutation equality operator.
Definition: permutation.h:450
std::map< index_type, index_type > Map
Definition: permutation.h:123
const_iterator cbegin() const
Begin element iterator factory function.
Definition: permutation.h:254
vector< vector< index_type > > cycles() const
Cycles decomposition.
Definition: permutation.h:320
bool is_in_domain(Integer i) const
Test if an index is in the domain of this permutation.
Definition: permutation.h:303
Permutation pow(int n) const
Raise this permutation to the n-th power.
Definition: permutation.h:397
Permutation(InIter first, InIter last)
Construct permutation using its 1-line form given by range [first,last)
Definition: permutation.h:200
Permutation(Permutation &&)=default
#define TA_ASSERT(EXPR,...)
Definition: error.h:39
bool operator<(const Permutation &p1, const Permutation &p2)
Permutation less-than operator.
Definition: permutation.h:470
index_type operator[](index_type e) const
Computes image of an element under this permutation.
Definition: permutation.h:273
Permutation of a sequence of objects indexed by base-0 indices.
Definition: permutation.h:117
Permutation operator-(const Permutation &perm)
Inverse permutation operator.
Definition: permutation.h:492
Map::const_iterator const_iterator
Definition: permutation.h:124
Permutation(Index &&a)
Construct permutation using 1-line form given as an integral range.
Definition: permutation.h:216
unsigned int domain_size() const
Permutation domain size.
Definition: permutation.h:235
Permutation(Map p)
Construct permutation using its compressed 2-line form given by std::map.
Definition: permutation.h:230
Permutation inv() const
Construct the inverse of this permutation.
Definition: permutation.h:381
Set domain() const
Computes the domain of this permutation.
Definition: permutation.h:286
const_iterator end() const
End element iterator factory function.
Definition: permutation.h:259
Permutation & operator=(const Permutation &)=default
Permutation mult(const Permutation &other) const
Product of this permutation by other.
Definition: permutation.h:360
Permutation & operator=(Permutation &&other)=default
const_iterator cend() const
End element iterator factory function.
Definition: permutation.h:264
Permutation identity() const
Idenity permutation factory method.
Definition: permutation.h:430
friend std::ostream & operator<<(std::ostream &output, const Permutation &p)
Add permutation to an output stream.
Definition: permutation.h:481
const Map & data() const
Data accessor.
Definition: permutation.h:425
Permutation operator^(const Permutation &perm, int n)
Raise perm to the n-th power.
Definition: permutation.h:516
container::svector< T > vector
Definition: permutation.h:122
const_iterator begin() const
Begin element iterator factory function.
Definition: permutation.h:249