unsupported/Eigen/CXX11/src/TensorSymmetry/DynamicSymmetry.h - mirror - Git at Google

 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
 // Copyright (C) 2013 Christian Seiler <christian@iwakd.de>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

 #ifndef EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H
 #define EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H

 // IWYU pragma: private
 #include "./InternalHeaderCheck.h"

 namespace Eigen {

 class DynamicSGroup {
  public:
   inline explicit DynamicSGroup() : m_numIndices(1), m_elements(), m_generators(), m_globalFlags(0) {
     m_elements.push_back(ge(Generator(0, 0, 0)));
   }
   inline DynamicSGroup(const DynamicSGroup& o)
       : m_numIndices(o.m_numIndices),
         m_elements(o.m_elements),
         m_generators(o.m_generators),
         m_globalFlags(o.m_globalFlags) {}
   inline DynamicSGroup(DynamicSGroup&& o)
       : m_numIndices(o.m_numIndices), m_elements(), m_generators(o.m_generators), m_globalFlags(o.m_globalFlags) {
     std::swap(m_elements, o.m_elements);
   }
   inline DynamicSGroup& operator=(const DynamicSGroup& o) {
     m_numIndices = o.m_numIndices;
     m_elements = o.m_elements;
     m_generators = o.m_generators;
     m_globalFlags = o.m_globalFlags;
     return *this;
   }
   inline DynamicSGroup& operator=(DynamicSGroup&& o) {
     m_numIndices = o.m_numIndices;
     std::swap(m_elements, o.m_elements);
     m_generators = o.m_generators;
     m_globalFlags = o.m_globalFlags;
     return *this;
   }

   void add(int one, int two, int flags = 0);

   template <typename Gen_>
   inline void add(Gen_) {
     add(Gen_::One, Gen_::Two, Gen_::Flags);
   }
   inline void addSymmetry(int one, int two) { add(one, two, 0); }
   inline void addAntiSymmetry(int one, int two) { add(one, two, NegationFlag); }
   inline void addHermiticity(int one, int two) { add(one, two, ConjugationFlag); }
   inline void addAntiHermiticity(int one, int two) { add(one, two, NegationFlag | ConjugationFlag); }

   template <typename Op, typename RV, typename Index, std::size_t N, typename... Args>
   inline RV apply(const std::array<Index, N>& idx, RV initial, Args&&... args) const {
     eigen_assert(N >= m_numIndices &&
                  "Can only apply symmetry group to objects that have at least the required amount of indices.");
     for (std::size_t i = 0; i < size(); i++)
       initial = Op::run(h_permute(i, idx, typename internal::gen_numeric_list<int, N>::type()), m_elements[i].flags,
                         initial, std::forward<Args>(args)...);
     return initial;
   }

   template <typename Op, typename RV, typename Index, typename... Args>
   inline RV apply(const std::vector<Index>& idx, RV initial, Args&&... args) const {
     eigen_assert(idx.size() >= m_numIndices &&
                  "Can only apply symmetry group to objects that have at least the required amount of indices.");
     for (std::size_t i = 0; i < size(); i++)
       initial = Op::run(h_permute(i, idx), m_elements[i].flags, initial, std::forward<Args>(args)...);
     return initial;
   }

   inline int globalFlags() const { return m_globalFlags; }
   inline std::size_t size() const { return m_elements.size(); }

   template <typename Tensor_, typename... IndexTypes>
   inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(Tensor_& tensor,
                                                                                    typename Tensor_::Index firstIndex,
                                                                                    IndexTypes... otherIndices) const {
     static_assert(sizeof...(otherIndices) + 1 == Tensor_::NumIndices,
                   "Number of indices used to access a tensor coefficient must be equal to the rank of the tensor.");
     return operator()(tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices>{{firstIndex, otherIndices...}});
   }

   template <typename Tensor_>
   inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(
       Tensor_& tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices> const& indices) const {
     return internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup>(tensor, *this, indices);
   }

  private:
   struct GroupElement {
     std::vector<int> representation;
     int flags;
     bool isId() const {
       for (std::size_t i = 0; i < representation.size(); i++)
         if (i != (size_t)representation[i]) return false;
       return true;
     }
   };
   struct Generator {
     int one;
     int two;
     int flags;
     constexpr inline Generator(int one_, int two_, int flags_) : one(one_), two(two_), flags(flags_) {}
   };

   std::size_t m_numIndices;
   std::vector<GroupElement> m_elements;
   std::vector<Generator> m_generators;
   int m_globalFlags;

   template <typename Index, std::size_t N, int... n>
   inline std::array<Index, N> h_permute(std::size_t which, const std::array<Index, N>& idx,
                                         internal::numeric_list<int, n...>) const {
     return std::array<Index, N>{{idx[n >= m_numIndices ? n : m_elements[which].representation[n]]...}};
   }

   template <typename Index>
   inline std::vector<Index> h_permute(std::size_t which, std::vector<Index> idx) const {
     std::vector<Index> result;
     result.reserve(idx.size());
     for (auto k : m_elements[which].representation) result.push_back(idx[k]);
     for (std::size_t i = m_numIndices; i < idx.size(); i++) result.push_back(idx[i]);
     return result;
   }

   inline GroupElement ge(Generator const& g) const {
     GroupElement result;
     result.representation.reserve(m_numIndices);
     result.flags = g.flags;
     for (std::size_t k = 0; k < m_numIndices; k++) {
       if (k == (std::size_t)g.one)
         result.representation.push_back(g.two);
       else if (k == (std::size_t)g.two)
         result.representation.push_back(g.one);
       else
         result.representation.push_back(int(k));
     }
     return result;
   }

   GroupElement mul(GroupElement, GroupElement) const;
   inline GroupElement mul(Generator g1, GroupElement g2) const { return mul(ge(g1), g2); }

   inline GroupElement mul(GroupElement g1, Generator g2) const { return mul(g1, ge(g2)); }

   inline GroupElement mul(Generator g1, Generator g2) const { return mul(ge(g1), ge(g2)); }

   inline int findElement(GroupElement e) const {
     for (auto ee : m_elements) {
       if (ee.representation == e.representation) return ee.flags ^ e.flags;
     }
     return -1;
   }

   void updateGlobalFlags(int flagDiffOfSameGenerator);
 };

 // dynamic symmetry group that auto-adds the template parameters in the constructor
 template <typename... Gen>
 class DynamicSGroupFromTemplateArgs : public DynamicSGroup {
  public:
   inline DynamicSGroupFromTemplateArgs() : DynamicSGroup() { add_all(internal::type_list<Gen...>()); }
   inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs const& other) : DynamicSGroup(other) {}
   inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs&& other) : DynamicSGroup(other) {}
   inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(const DynamicSGroupFromTemplateArgs<Gen...>& o) {
     DynamicSGroup::operator=(o);
     return *this;
   }
   inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(DynamicSGroupFromTemplateArgs<Gen...>&& o) {
     DynamicSGroup::operator=(o);
     return *this;
   }

  private:
   template <typename Gen1, typename... GenNext>
   inline void add_all(internal::type_list<Gen1, GenNext...>) {
     add(Gen1());
     add_all(internal::type_list<GenNext...>());
   }

   inline void add_all(internal::type_list<>) {}
 };

 inline DynamicSGroup::GroupElement DynamicSGroup::mul(GroupElement g1, GroupElement g2) const {
   eigen_internal_assert(g1.representation.size() == m_numIndices);
   eigen_internal_assert(g2.representation.size() == m_numIndices);

   GroupElement result;
   result.representation.reserve(m_numIndices);
   for (std::size_t i = 0; i < m_numIndices; i++) {
     int v = g2.representation[g1.representation[i]];
     eigen_assert(v >= 0);
     result.representation.push_back(v);
   }
   result.flags = g1.flags ^ g2.flags;
   return result;
 }

 inline void DynamicSGroup::add(int one, int two, int flags) {
   eigen_assert(one >= 0);
   eigen_assert(two >= 0);
   eigen_assert(one != two);

   if ((std::size_t)one >= m_numIndices || (std::size_t)two >= m_numIndices) {
     std::size_t newNumIndices = (one > two) ? one : two + 1;
     for (auto& gelem : m_elements) {
       gelem.representation.reserve(newNumIndices);
       for (std::size_t i = m_numIndices; i < newNumIndices; i++) gelem.representation.push_back(i);
     }
     m_numIndices = newNumIndices;
   }

   Generator g{one, two, flags};
   GroupElement e = ge(g);

   /* special case for first generator */
   if (m_elements.size() == 1) {
     while (!e.isId()) {
       m_elements.push_back(e);
       e = mul(e, g);
     }

     if (e.flags > 0) updateGlobalFlags(e.flags);

     // only add in case we didn't have identity
     if (m_elements.size() > 1) m_generators.push_back(g);
     return;
   }

   int p = findElement(e);
   if (p >= 0) {
     updateGlobalFlags(p);
     return;
   }

   std::size_t coset_order = m_elements.size();
   m_elements.push_back(e);
   for (std::size_t i = 1; i < coset_order; i++) m_elements.push_back(mul(m_elements[i], e));
   m_generators.push_back(g);

   std::size_t coset_rep = coset_order;
   do {
     for (auto g : m_generators) {
       e = mul(m_elements[coset_rep], g);
       p = findElement(e);
       if (p < 0) {
         // element not yet in group
         m_elements.push_back(e);
         for (std::size_t i = 1; i < coset_order; i++) m_elements.push_back(mul(m_elements[i], e));
       } else if (p > 0) {
         updateGlobalFlags(p);
       }
     }
     coset_rep += coset_order;
   } while (coset_rep < m_elements.size());
 }

 inline void DynamicSGroup::updateGlobalFlags(int flagDiffOfSameGenerator) {
   switch (flagDiffOfSameGenerator) {
     case 0:
     default:
       // nothing happened
       break;
     case NegationFlag:
       // every element is it's own negative => whole tensor is zero
       m_globalFlags |= GlobalZeroFlag;
       break;
     case ConjugationFlag:
       // every element is it's own conjugate => whole tensor is real
       m_globalFlags |= GlobalRealFlag;
       break;
     case (NegationFlag | ConjugationFlag):
       // every element is it's own negative conjugate => whole tensor is imaginary
       m_globalFlags |= GlobalImagFlag;
       break;
       /* NOTE:
        *   since GlobalZeroFlag == GlobalRealFlag | GlobalImagFlag, if one generator
        *   causes the tensor to be real and the next one to be imaginary, this will
        *   trivially give the correct result
        */
   }
 }

 }  // end namespace Eigen

 #endif  // EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H

 /*
  * kate: space-indent on; indent-width 2; mixedindent off; indent-mode cstyle;
  */
	// This file is part of Eigen, a lightweight C++ template library
	// for linear algebra.
	//
	// Copyright (C) 2013 Christian Seiler <christian@iwakd.de>
	//
	// This Source Code Form is subject to the terms of the Mozilla
	// Public License v. 2.0. If a copy of the MPL was not distributed
	// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

	#ifndef EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H
	#define EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H

	// IWYU pragma: private
	#include "./InternalHeaderCheck.h"

	namespace Eigen {

	class DynamicSGroup {
	public:
	inline explicit DynamicSGroup() : m_numIndices(1), m_elements(), m_generators(), m_globalFlags(0) {
	m_elements.push_back(ge(Generator(0, 0, 0)));
	}
	inline DynamicSGroup(const DynamicSGroup& o)
	: m_numIndices(o.m_numIndices),
	m_elements(o.m_elements),
	m_generators(o.m_generators),
	m_globalFlags(o.m_globalFlags) {}
	inline DynamicSGroup(DynamicSGroup&& o)
	: m_numIndices(o.m_numIndices), m_elements(), m_generators(o.m_generators), m_globalFlags(o.m_globalFlags) {
	std::swap(m_elements, o.m_elements);
	}
	inline DynamicSGroup& operator=(const DynamicSGroup& o) {
	m_numIndices = o.m_numIndices;
	m_elements = o.m_elements;
	m_generators = o.m_generators;
	m_globalFlags = o.m_globalFlags;
	return *this;
	}
	inline DynamicSGroup& operator=(DynamicSGroup&& o) {
	m_numIndices = o.m_numIndices;
	std::swap(m_elements, o.m_elements);
	m_generators = o.m_generators;
	m_globalFlags = o.m_globalFlags;
	return *this;
	}

	void add(int one, int two, int flags = 0);

	template <typename Gen_>
	inline void add(Gen_) {
	add(Gen_::One, Gen_::Two, Gen_::Flags);
	}
	inline void addSymmetry(int one, int two) { add(one, two, 0); }
	inline void addAntiSymmetry(int one, int two) { add(one, two, NegationFlag); }
	inline void addHermiticity(int one, int two) { add(one, two, ConjugationFlag); }
	inline void addAntiHermiticity(int one, int two) { add(one, two, NegationFlag \| ConjugationFlag); }

	template <typename Op, typename RV, typename Index, std::size_t N, typename... Args>
	inline RV apply(const std::array<Index, N>& idx, RV initial, Args&&... args) const {
	eigen_assert(N >= m_numIndices &&
	"Can only apply symmetry group to objects that have at least the required amount of indices.");
	for (std::size_t i = 0; i < size(); i++)
	initial = Op::run(h_permute(i, idx, typename internal::gen_numeric_list<int, N>::type()), m_elements[i].flags,
	initial, std::forward<Args>(args)...);
	return initial;
	}

	template <typename Op, typename RV, typename Index, typename... Args>
	inline RV apply(const std::vector<Index>& idx, RV initial, Args&&... args) const {
	eigen_assert(idx.size() >= m_numIndices &&
	"Can only apply symmetry group to objects that have at least the required amount of indices.");
	for (std::size_t i = 0; i < size(); i++)
	initial = Op::run(h_permute(i, idx), m_elements[i].flags, initial, std::forward<Args>(args)...);
	return initial;
	}

	inline int globalFlags() const { return m_globalFlags; }
	inline std::size_t size() const { return m_elements.size(); }

	template <typename Tensor_, typename... IndexTypes>
	inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(Tensor_& tensor,
	typename Tensor_::Index firstIndex,
	IndexTypes... otherIndices) const {
	static_assert(sizeof...(otherIndices) + 1 == Tensor_::NumIndices,
	"Number of indices used to access a tensor coefficient must be equal to the rank of the tensor.");
	return operator()(tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices>{{firstIndex, otherIndices...}});
	}

	template <typename Tensor_>
	inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(
	Tensor_& tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices> const& indices) const {
	return internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup>(tensor, *this, indices);
	}

	private:
	struct GroupElement {
	std::vector<int> representation;
	int flags;
	bool isId() const {
	for (std::size_t i = 0; i < representation.size(); i++)
	if (i != (size_t)representation[i]) return false;
	return true;
	}
	};
	struct Generator {
	int one;
	int two;
	int flags;
	constexpr inline Generator(int one_, int two_, int flags_) : one(one_), two(two_), flags(flags_) {}
	};

	std::size_t m_numIndices;
	std::vector<GroupElement> m_elements;
	std::vector<Generator> m_generators;
	int m_globalFlags;

	template <typename Index, std::size_t N, int... n>
	inline std::array<Index, N> h_permute(std::size_t which, const std::array<Index, N>& idx,
	internal::numeric_list<int, n...>) const {
	return std::array<Index, N>{{idx[n >= m_numIndices ? n : m_elements[which].representation[n]]...}};
	}

	template <typename Index>
	inline std::vector<Index> h_permute(std::size_t which, std::vector<Index> idx) const {
	std::vector<Index> result;
	result.reserve(idx.size());
	for (auto k : m_elements[which].representation) result.push_back(idx[k]);
	for (std::size_t i = m_numIndices; i < idx.size(); i++) result.push_back(idx[i]);
	return result;
	}

	inline GroupElement ge(Generator const& g) const {
	GroupElement result;
	result.representation.reserve(m_numIndices);
	result.flags = g.flags;
	for (std::size_t k = 0; k < m_numIndices; k++) {
	if (k == (std::size_t)g.one)
	result.representation.push_back(g.two);
	else if (k == (std::size_t)g.two)
	result.representation.push_back(g.one);
	else
	result.representation.push_back(int(k));
	}
	return result;
	}

	GroupElement mul(GroupElement, GroupElement) const;
	inline GroupElement mul(Generator g1, GroupElement g2) const { return mul(ge(g1), g2); }

	inline GroupElement mul(GroupElement g1, Generator g2) const { return mul(g1, ge(g2)); }

	inline GroupElement mul(Generator g1, Generator g2) const { return mul(ge(g1), ge(g2)); }

	inline int findElement(GroupElement e) const {
	for (auto ee : m_elements) {
	if (ee.representation == e.representation) return ee.flags ^ e.flags;
	}
	return -1;
	}

	void updateGlobalFlags(int flagDiffOfSameGenerator);
	};

	// dynamic symmetry group that auto-adds the template parameters in the constructor
	template <typename... Gen>
	class DynamicSGroupFromTemplateArgs : public DynamicSGroup {
	public:
	inline DynamicSGroupFromTemplateArgs() : DynamicSGroup() { add_all(internal::type_list<Gen...>()); }
	inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs const& other) : DynamicSGroup(other) {}
	inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs&& other) : DynamicSGroup(other) {}
	inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(const DynamicSGroupFromTemplateArgs<Gen...>& o) {
	DynamicSGroup::operator=(o);
	return *this;
	}
	inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(DynamicSGroupFromTemplateArgs<Gen...>&& o) {
	DynamicSGroup::operator=(o);
	return *this;
	}

	private:
	template <typename Gen1, typename... GenNext>
	inline void add_all(internal::type_list<Gen1, GenNext...>) {
	add(Gen1());
	add_all(internal::type_list<GenNext...>());
	}

	inline void add_all(internal::type_list<>) {}
	};

	inline DynamicSGroup::GroupElement DynamicSGroup::mul(GroupElement g1, GroupElement g2) const {
	eigen_internal_assert(g1.representation.size() == m_numIndices);
	eigen_internal_assert(g2.representation.size() == m_numIndices);

	GroupElement result;
	result.representation.reserve(m_numIndices);
	for (std::size_t i = 0; i < m_numIndices; i++) {
	int v = g2.representation[g1.representation[i]];
	eigen_assert(v >= 0);
	result.representation.push_back(v);
	}
	result.flags = g1.flags ^ g2.flags;
	return result;
	}

	inline void DynamicSGroup::add(int one, int two, int flags) {
	eigen_assert(one >= 0);
	eigen_assert(two >= 0);
	eigen_assert(one != two);

	if ((std::size_t)one >= m_numIndices \|\| (std::size_t)two >= m_numIndices) {
	std::size_t newNumIndices = (one > two) ? one : two + 1;
	for (auto& gelem : m_elements) {
	gelem.representation.reserve(newNumIndices);
	for (std::size_t i = m_numIndices; i < newNumIndices; i++) gelem.representation.push_back(i);
	}
	m_numIndices = newNumIndices;
	}

	Generator g{one, two, flags};
	GroupElement e = ge(g);

	/* special case for first generator */
	if (m_elements.size() == 1) {
	while (!e.isId()) {
	m_elements.push_back(e);
	e = mul(e, g);
	}

	if (e.flags > 0) updateGlobalFlags(e.flags);

	// only add in case we didn't have identity
	if (m_elements.size() > 1) m_generators.push_back(g);
	return;
	}

	int p = findElement(e);
	if (p >= 0) {
	updateGlobalFlags(p);
	return;
	}

	std::size_t coset_order = m_elements.size();
	m_elements.push_back(e);
	for (std::size_t i = 1; i < coset_order; i++) m_elements.push_back(mul(m_elements[i], e));
	m_generators.push_back(g);

	std::size_t coset_rep = coset_order;
	do {
	for (auto g : m_generators) {
	e = mul(m_elements[coset_rep], g);
	p = findElement(e);
	if (p < 0) {
	// element not yet in group
	m_elements.push_back(e);
	for (std::size_t i = 1; i < coset_order; i++) m_elements.push_back(mul(m_elements[i], e));
	} else if (p > 0) {
	updateGlobalFlags(p);
	}
	}
	coset_rep += coset_order;
	} while (coset_rep < m_elements.size());
	}

	inline void DynamicSGroup::updateGlobalFlags(int flagDiffOfSameGenerator) {
	switch (flagDiffOfSameGenerator) {
	case 0:
	default:
	// nothing happened
	break;
	case NegationFlag:
	// every element is it's own negative => whole tensor is zero
	m_globalFlags \|= GlobalZeroFlag;
	break;
	case ConjugationFlag:
	// every element is it's own conjugate => whole tensor is real
	m_globalFlags \|= GlobalRealFlag;
	break;
	case (NegationFlag \| ConjugationFlag):
	// every element is it's own negative conjugate => whole tensor is imaginary
	m_globalFlags \|= GlobalImagFlag;
	break;
	/* NOTE:
	* since GlobalZeroFlag == GlobalRealFlag \| GlobalImagFlag, if one generator
	* causes the tensor to be real and the next one to be imaginary, this will
	* trivially give the correct result
	*/
	}
	}

	} // end namespace Eigen

	#endif // EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H

	/*
	* kate: space-indent on; indent-width 2; mixedindent off; indent-mode cstyle;
	*/