DynamicSymmetry.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2013 Christian Seiler <christian@iwakd.de>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H
00011 #define EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H
00012 
00013 namespace Eigen {
00014 
00015 class DynamicSGroup
00016 {
00017   public:
00018     inline explicit DynamicSGroup() : m_numIndices(1), m_elements(), m_generators(), m_globalFlags(0) { m_elements.push_back(ge(Generator(0, 0, 0))); }
00019     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) { }
00020     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); }
00021     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; }
00022     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; }
00023 
00024     void add(int one, int two, int flags = 0);
00025 
00026     template<typename Gen_>
00027     inline void add(Gen_) { add(Gen_::One, Gen_::Two, Gen_::Flags); }
00028     inline void addSymmetry(int one, int two) { add(one, two, 0); }
00029     inline void addAntiSymmetry(int one, int two) { add(one, two, NegationFlag); }
00030     inline void addHermiticity(int one, int two) { add(one, two, ConjugationFlag); }
00031     inline void addAntiHermiticity(int one, int two) { add(one, two, NegationFlag | ConjugationFlag); }
00032 
00033     template<typename Op, typename RV, typename Index, std::size_t N, typename... Args>
00034     inline RV apply(const std::array<Index, N>& idx, RV initial, Args&&... args) const
00035     {
00036       eigen_assert(N >= m_numIndices && "Can only apply symmetry group to objects that have at least the required amount of indices.");
00037       for (std::size_t i = 0; i < size(); i++)
00038         initial = Op::run(h_permute(i, idx, typename internal::gen_numeric_list<int, N>::type()), m_elements[i].flags, initial, std::forward<Args>(args)...);
00039       return initial;
00040     }
00041 
00042     template<typename Op, typename RV, typename Index, typename... Args>
00043     inline RV apply(const std::vector<Index>& idx, RV initial, Args&&... args) const
00044     {
00045       eigen_assert(idx.size() >= m_numIndices && "Can only apply symmetry group to objects that have at least the required amount of indices.");
00046       for (std::size_t i = 0; i < size(); i++)
00047         initial = Op::run(h_permute(i, idx), m_elements[i].flags, initial, std::forward<Args>(args)...);
00048       return initial;
00049     }
00050 
00051     inline int globalFlags() const { return m_globalFlags; }
00052     inline std::size_t size() const { return m_elements.size(); }
00053 
00054     template<typename Tensor_, typename... IndexTypes>
00055     inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(Tensor_& tensor, typename Tensor_::Index firstIndex, IndexTypes... otherIndices) const
00056     {
00057       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.");
00058       return operator()(tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices>{{firstIndex, otherIndices...}});
00059     }
00060 
00061     template<typename Tensor_>
00062     inline internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup> operator()(Tensor_& tensor, std::array<typename Tensor_::Index, Tensor_::NumIndices> const& indices) const
00063     {
00064       return internal::tensor_symmetry_value_setter<Tensor_, DynamicSGroup>(tensor, *this, indices);
00065     }
00066   private:
00067     struct GroupElement {
00068       std::vector<int> representation;
00069       int flags;
00070       bool isId() const
00071       {
00072         for (std::size_t i = 0; i < representation.size(); i++)
00073           if (i != (size_t)representation[i])
00074             return false;
00075         return true;
00076       }
00077     };
00078     struct Generator {
00079       int one;
00080       int two;
00081       int flags;
00082       constexpr inline Generator(int one_, int two_, int flags_) : one(one_), two(two_), flags(flags_) {}
00083     };
00084 
00085     std::size_t m_numIndices;
00086     std::vector<GroupElement> m_elements;
00087     std::vector<Generator> m_generators;
00088     int m_globalFlags;
00089 
00090     template<typename Index, std::size_t N, int... n>
00091     inline std::array<Index, N> h_permute(std::size_t which, const std::array<Index, N>& idx, internal::numeric_list<int, n...>) const
00092     {
00093       return std::array<Index, N>{{ idx[n >= m_numIndices ? n : m_elements[which].representation[n]]... }};
00094     }
00095 
00096     template<typename Index>
00097     inline std::vector<Index> h_permute(std::size_t which, std::vector<Index> idx) const
00098     {
00099       std::vector<Index> result;
00100       result.reserve(idx.size());
00101       for (auto k : m_elements[which].representation)
00102         result.push_back(idx[k]);
00103       for (std::size_t i = m_numIndices; i < idx.size(); i++)
00104         result.push_back(idx[i]);
00105       return result;
00106     }
00107 
00108     inline GroupElement ge(Generator const& g) const
00109     {
00110       GroupElement result;
00111       result.representation.reserve(m_numIndices);
00112       result.flags = g.flags;
00113       for (std::size_t k = 0; k < m_numIndices; k++) {
00114         if (k == (std::size_t)g.one)
00115           result.representation.push_back(g.two);
00116         else if (k == (std::size_t)g.two)
00117           result.representation.push_back(g.one);
00118         else
00119           result.representation.push_back(int(k));
00120       }
00121       return result;
00122     }
00123 
00124     GroupElement mul(GroupElement, GroupElement) const;
00125     inline GroupElement mul(Generator g1, GroupElement g2) const
00126     {
00127       return mul(ge(g1), g2);
00128     }
00129 
00130     inline GroupElement mul(GroupElement g1, Generator g2) const
00131     {
00132       return mul(g1, ge(g2));
00133     }
00134 
00135     inline GroupElement mul(Generator g1, Generator g2) const
00136     {
00137       return mul(ge(g1), ge(g2));
00138     }
00139 
00140     inline int findElement(GroupElement e) const
00141     {
00142       for (auto ee : m_elements) {
00143         if (ee.representation == e.representation)
00144           return ee.flags ^ e.flags;
00145       }
00146       return -1;
00147     }
00148 
00149     void updateGlobalFlags(int flagDiffOfSameGenerator);
00150 };
00151 
00152 // dynamic symmetry group that auto-adds the template parameters in the constructor
00153 template<typename... Gen>
00154 class DynamicSGroupFromTemplateArgs : public DynamicSGroup
00155 {
00156   public:
00157     inline DynamicSGroupFromTemplateArgs() : DynamicSGroup()
00158     {
00159       add_all(internal::type_list<Gen...>());
00160     }
00161     inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs const& other) : DynamicSGroup(other) { }
00162     inline DynamicSGroupFromTemplateArgs(DynamicSGroupFromTemplateArgs&& other) : DynamicSGroup(other) { }
00163     inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(const DynamicSGroupFromTemplateArgs<Gen...>& o) { DynamicSGroup::operator=(o); return *this; }
00164     inline DynamicSGroupFromTemplateArgs<Gen...>& operator=(DynamicSGroupFromTemplateArgs<Gen...>&& o) { DynamicSGroup::operator=(o); return *this; }
00165   
00166   private:
00167     template<typename Gen1, typename... GenNext>
00168     inline void add_all(internal::type_list<Gen1, GenNext...>)
00169     {
00170       add(Gen1());
00171       add_all(internal::type_list<GenNext...>());
00172     }
00173 
00174     inline void add_all(internal::type_list<>)
00175     {
00176     }
00177 };
00178 
00179 inline DynamicSGroup::GroupElement DynamicSGroup::mul(GroupElement g1, GroupElement g2) const
00180 {
00181   eigen_internal_assert(g1.representation.size() == m_numIndices);
00182   eigen_internal_assert(g2.representation.size() == m_numIndices);
00183 
00184   GroupElement result;
00185   result.representation.reserve(m_numIndices);
00186   for (std::size_t i = 0; i < m_numIndices; i++) {
00187     int v = g2.representation[g1.representation[i]];
00188     eigen_assert(v >= 0);
00189     result.representation.push_back(v);
00190   }
00191   result.flags = g1.flags ^ g2.flags;
00192   return result;
00193 }
00194 
00195 inline void DynamicSGroup::add(int one, int two, int flags)
00196 {
00197   eigen_assert(one >= 0);
00198   eigen_assert(two >= 0);
00199   eigen_assert(one != two);
00200 
00201   if ((std::size_t)one >= m_numIndices || (std::size_t)two >= m_numIndices) {
00202     std::size_t newNumIndices = (one > two) ? one : two + 1;
00203     for (auto& gelem : m_elements) {
00204       gelem.representation.reserve(newNumIndices);
00205       for (std::size_t i = m_numIndices; i < newNumIndices; i++)
00206         gelem.representation.push_back(i);
00207     }
00208     m_numIndices = newNumIndices;
00209   }
00210 
00211   Generator g{one, two, flags};
00212   GroupElement e = ge(g);
00213 
00214   /* special case for first generator */
00215   if (m_elements.size() == 1) {
00216     while (!e.isId()) {
00217       m_elements.push_back(e);
00218       e = mul(e, g);
00219     }
00220 
00221     if (e.flags > 0)
00222       updateGlobalFlags(e.flags);
00223 
00224     // only add in case we didn't have identity
00225     if (m_elements.size() > 1)
00226       m_generators.push_back(g);
00227     return;
00228   }
00229 
00230   int p = findElement(e);
00231   if (p >= 0) {
00232     updateGlobalFlags(p);
00233     return;
00234   }
00235 
00236   std::size_t coset_order = m_elements.size();
00237   m_elements.push_back(e);
00238   for (std::size_t i = 1; i < coset_order; i++)
00239     m_elements.push_back(mul(m_elements[i], e));
00240   m_generators.push_back(g);
00241 
00242   std::size_t coset_rep = coset_order;
00243   do {
00244     for (auto g : m_generators) {
00245       e = mul(m_elements[coset_rep], g);
00246       p = findElement(e);
00247       if (p < 0) {
00248         // element not yet in group
00249         m_elements.push_back(e);
00250         for (std::size_t i = 1; i < coset_order; i++)
00251           m_elements.push_back(mul(m_elements[i], e));
00252       } else if (p > 0) {
00253         updateGlobalFlags(p);
00254       }
00255     }
00256     coset_rep += coset_order;
00257   } while (coset_rep < m_elements.size());
00258 }
00259 
00260 inline void DynamicSGroup::updateGlobalFlags(int flagDiffOfSameGenerator)
00261 {
00262     switch (flagDiffOfSameGenerator) {
00263       case 0:
00264       default:
00265         // nothing happened
00266         break;
00267       case NegationFlag:
00268         // every element is it's own negative => whole tensor is zero
00269         m_globalFlags |= GlobalZeroFlag;
00270         break;
00271       case ConjugationFlag:
00272         // every element is it's own conjugate => whole tensor is real
00273         m_globalFlags |= GlobalRealFlag;
00274         break;
00275       case (NegationFlag | ConjugationFlag):
00276         // every element is it's own negative conjugate => whole tensor is imaginary
00277         m_globalFlags |= GlobalImagFlag;
00278         break;
00279       /* NOTE:
00280        *   since GlobalZeroFlag == GlobalRealFlag | GlobalImagFlag, if one generator
00281        *   causes the tensor to be real and the next one to be imaginary, this will
00282        *   trivially give the correct result
00283        */
00284     }
00285 }
00286 
00287 } // end namespace Eigen
00288 
00289 #endif // EIGEN_CXX11_TENSORSYMMETRY_DYNAMICSYMMETRY_H
00290 
00291 /*
00292  * kate: space-indent on; indent-width 2; mixedindent off; indent-mode cstyle;
00293  */
 All Classes Functions Variables Typedefs Enumerator