Eigen-Contrib  5.0.1
 
Loading...
Searching...
No Matches
DynamicSymmetry.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2013 Christian Seiler <christian@iwakd.de>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9// SPDX-License-Identifier: MPL-2.0
10
11#ifndef EIGEN_TENSORSYMMETRY_DYNAMICSYMMETRY_H
12#define EIGEN_TENSORSYMMETRY_DYNAMICSYMMETRY_H
13
14// IWYU pragma: private
15#include "./InternalHeaderCheck.h"
16
17namespace Eigen {
18
19class DynamicSGroup {
20 public:
21 inline explicit DynamicSGroup() : m_numIndices(1), m_elements(), m_generators(), m_globalFlags(0) {
22 m_elements.push_back(ge(Generator(0, 0, 0)));
23 }
24 inline DynamicSGroup(const DynamicSGroup& o)
25 : m_numIndices(o.m_numIndices),
26 m_elements(o.m_elements),
27 m_generators(o.m_generators),
28 m_globalFlags(o.m_globalFlags) {}
29 inline DynamicSGroup(DynamicSGroup&& o)
30 : m_numIndices(o.m_numIndices), m_elements(), m_generators(o.m_generators), m_globalFlags(o.m_globalFlags) {
31 std::swap(m_elements, o.m_elements);
32 }
33 inline DynamicSGroup& operator=(const DynamicSGroup& o) {
34 m_numIndices = o.m_numIndices;
35 m_elements = o.m_elements;
36 m_generators = o.m_generators;
37 m_globalFlags = o.m_globalFlags;
38 return *this;
39 }
40 inline DynamicSGroup& operator=(DynamicSGroup&& o) {
41 m_numIndices = o.m_numIndices;
42 std::swap(m_elements, o.m_elements);
43 m_generators = o.m_generators;
44 m_globalFlags = o.m_globalFlags;
45 return *this;
46 }
47
48 void add(int one, int two, int flags = 0);
49
50 template <typename Gen_>
51 inline void add(Gen_) {
52 add(Gen_::One, Gen_::Two, Gen_::Flags);
53 }
54 inline void addSymmetry(int one, int two) { add(one, two, 0); }
55 inline void addAntiSymmetry(int one, int two) { add(one, two, NegationFlag); }
56 inline void addHermiticity(int one, int two) { add(one, two, ConjugationFlag); }
57 inline void addAntiHermiticity(int one, int two) { add(one, two, NegationFlag | ConjugationFlag); }
58
59 template <typename Op, typename RV, typename Index, std::size_t N, typename... Args>
60 inline RV apply(const std::array<Index, N>& idx, RV initial, Args&&... args) const {
61 eigen_assert(N >= m_numIndices &&
62 "Can only apply symmetry group to objects that have at least the required amount of indices.");
63 for (std::size_t i = 0; i < size(); i++)
64 initial = Op::run(h_permute(i, idx, std::make_integer_sequence<int, N>{}), m_elements[i].flags, initial,
65 std::forward<Args>(args)...);
66 return initial;
67 }
68
69 template <typename Op, typename RV, typename Index, typename... Args>
70 inline RV apply(const std::vector<Index>& idx, RV initial, Args&&... args) const {
71 eigen_assert(idx.size() >= m_numIndices &&
72 "Can only apply symmetry group to objects that have at least the required amount of indices.");
73 for (std::size_t i = 0; i < size(); i++)
74 initial = Op::run(h_permute(i, idx), m_elements[i].flags, initial, std::forward<Args>(args)...);
75 return initial;
76 }
77
78 inline int globalFlags() const { return m_globalFlags; }
79 inline std::size_t size() const { return m_elements.size(); }
80
81 template <typename Tensor_, typename... IndexTypes>
82 inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(Tensor_& tensor,
83 typename Tensor_::Index firstIndex,
84 IndexTypes... otherIndices) const {
85 static_assert(sizeof...(otherIndices) + 1 == Tensor_::NumIndices,
86 "Number of indices used to access a tensor coefficient must be equal to the rank of the tensor.");
87 return operator()(tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices>{{firstIndex, otherIndices...}});
88 }
89
90 template <typename Tensor_>
91 inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(
92 Tensor_& tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices> const& indices) const {
93 return internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup>(tensor, *this, indices);
94 }
95
96 private:
97 struct GroupElement {
98 std::vector<int> representation;
99 int flags;
100 bool isId() const {
101 for (std::size_t i = 0; i < representation.size(); i++)
102 if (i != (size_t)representation[i]) return false;
103 return true;
104 }
105 };
106 struct Generator {
107 int one;
108 int two;
109 int flags;
110 constexpr Generator(int one_, int two_, int flags_) : one(one_), two(two_), flags(flags_) {}
111 };
112
113 std::size_t m_numIndices;
114 std::vector<GroupElement> m_elements;
115 std::vector<Generator> m_generators;
116 int m_globalFlags;
117
118 template <typename Index, std::size_t N, int... n>
119 inline std::array<Index, N> h_permute(std::size_t which, const std::array<Index, N>& idx,
120 std::integer_sequence<int, n...>) const {
121 return std::array<Index, N>{{idx[n >= m_numIndices ? n : m_elements[which].representation[n]]...}};
122 }
123
124 template <typename Index>
125 inline std::vector<Index> h_permute(std::size_t which, const std::vector<Index>& idx) const {
126 std::vector<Index> result;
127 result.reserve(idx.size());
128 for (auto k : m_elements[which].representation) result.push_back(idx[k]);
129 for (std::size_t i = m_numIndices; i < idx.size(); i++) result.push_back(idx[i]);
130 return result;
131 }
132
133 inline GroupElement ge(Generator const& g) const {
134 GroupElement result;
135 result.representation.reserve(m_numIndices);
136 result.flags = g.flags;
137 for (std::size_t k = 0; k < m_numIndices; k++) {
138 if (k == (std::size_t)g.one)
139 result.representation.push_back(g.two);
140 else if (k == (std::size_t)g.two)
141 result.representation.push_back(g.one);
142 else
143 result.representation.push_back(int(k));
144 }
145 return result;
146 }
147
148 GroupElement mul(GroupElement, GroupElement) const;
149 inline GroupElement mul(Generator g1, GroupElement g2) const { return mul(ge(g1), g2); }
150
151 inline GroupElement mul(GroupElement g1, Generator g2) const { return mul(g1, ge(g2)); }
152
153 inline GroupElement mul(Generator g1, Generator g2) const { return mul(ge(g1), ge(g2)); }
154
155 inline int findElement(GroupElement e) const {
156 for (auto ee : m_elements) {
157 if (ee.representation == e.representation) return ee.flags ^ e.flags;
158 }
159 return -1;
160 }
161
162 void updateGlobalFlags(int flagDiffOfSameGenerator);
163};
164
165// dynamic symmetry group that auto-adds the template parameters in the constructor
166template <typename... Gen>
167class DynamicSGroupFromTemplateArgs : public DynamicSGroup {
168 public:
169 inline DynamicSGroupFromTemplateArgs() : DynamicSGroup() { add_all(internal::type_list<Gen...>()); }
170 inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs const& other) = default;
171 inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs&& other) : DynamicSGroup(other) {}
172 inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(const DynamicSGroupFromTemplateArgs<Gen...>& o) {
173 DynamicSGroup::operator=(o);
174 return *this;
175 }
176 inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(DynamicSGroupFromTemplateArgs<Gen...>&& o) {
177 DynamicSGroup::operator=(o);
178 return *this;
179 }
180
181 private:
182 template <typename Gen1, typename... GenNext>
183 inline void add_all(internal::type_list<Gen1, GenNext...>) {
184 add(Gen1());
185 add_all(internal::type_list<GenNext...>());
186 }
187
188 inline void add_all(internal::type_list<>) {}
189};
190
191inline DynamicSGroup::GroupElement DynamicSGroup::mul(GroupElement g1, GroupElement g2) const {
192 eigen_internal_assert(g1.representation.size() == m_numIndices);
193 eigen_internal_assert(g2.representation.size() == m_numIndices);
194
195 GroupElement result;
196 result.representation.reserve(m_numIndices);
197 for (std::size_t i = 0; i < m_numIndices; i++) {
198 int v = g2.representation[g1.representation[i]];
199 eigen_assert(v >= 0);
200 result.representation.push_back(v);
201 }
202 result.flags = g1.flags ^ g2.flags;
203 return result;
204}
205
206inline void DynamicSGroup::add(int one, int two, int flags) {
207 eigen_assert(one >= 0);
208 eigen_assert(two >= 0);
209 eigen_assert(one != two);
210
211 if ((std::size_t)one >= m_numIndices || (std::size_t)two >= m_numIndices) {
212 std::size_t newNumIndices = (one > two) ? one : two + 1;
213 for (auto& gelem : m_elements) {
214 gelem.representation.reserve(newNumIndices);
215 for (std::size_t i = m_numIndices; i < newNumIndices; i++) gelem.representation.push_back(static_cast<int>(i));
216 }
217 m_numIndices = newNumIndices;
218 }
219
220 Generator g{one, two, flags};
221 GroupElement e = ge(g);
222
223 /* special case for first generator */
224 if (m_elements.size() == 1) {
225 while (!e.isId()) {
226 m_elements.push_back(e);
227 e = mul(e, g);
228 }
229
230 if (e.flags > 0) updateGlobalFlags(e.flags);
231
232 // only add in case we didn't have identity
233 if (m_elements.size() > 1) m_generators.push_back(g);
234 return;
235 }
236
237 int p = findElement(e);
238 if (p >= 0) {
239 updateGlobalFlags(p);
240 return;
241 }
242
243 std::size_t coset_order = m_elements.size();
244 m_elements.push_back(e);
245 for (std::size_t i = 1; i < coset_order; i++) m_elements.push_back(mul(m_elements[i], e));
246 m_generators.push_back(g);
247
248 std::size_t coset_rep = coset_order;
249 do {
250 for (auto gen : m_generators) {
251 e = mul(m_elements[coset_rep], gen);
252 p = findElement(e);
253 if (p < 0) {
254 // element not yet in group
255 m_elements.push_back(e);
256 for (std::size_t i = 1; i < coset_order; i++) m_elements.push_back(mul(m_elements[i], e));
257 } else if (p > 0) {
258 updateGlobalFlags(p);
259 }
260 }
261 coset_rep += coset_order;
262 } while (coset_rep < m_elements.size());
263}
264
265inline void DynamicSGroup::updateGlobalFlags(int flagDiffOfSameGenerator) {
266 switch (flagDiffOfSameGenerator) {
267 case 0:
268 default:
269 // nothing happened
270 break;
271 case NegationFlag:
272 // every element is its own negative => whole tensor is zero
273 m_globalFlags |= GlobalZeroFlag;
274 break;
275 case ConjugationFlag:
276 // every element is its own conjugate => whole tensor is real
277 m_globalFlags |= GlobalRealFlag;
278 break;
279 case (NegationFlag | ConjugationFlag):
280 // every element is its own negative conjugate => whole tensor is imaginary
281 m_globalFlags |= GlobalImagFlag;
282 break;
283 /* NOTE:
284 * since GlobalZeroFlag == GlobalRealFlag | GlobalImagFlag, if one generator
285 * causes the tensor to be real and the next one to be imaginary, this will
286 * trivially give the correct result
287 */
288 }
289}
290
291} // end namespace Eigen
292
293#endif // EIGEN_TENSORSYMMETRY_DYNAMICSYMMETRY_H
294
295/*
296 * kate: space-indent on; indent-width 2; mixedindent off; indent-mode cstyle;
297 */
Eigen::DynamicSGroup
Dynamic symmetry group.
Definition DynamicSymmetry.h:19
Eigen
Namespace containing all symbols from the Eigen library.