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