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