source: git/gfanlib/gfanlib_symmetry.cpp @ 74a91c9

jengelh-datetimespielwiese
Last change on this file since 74a91c9 was 74a91c9, checked in by Frank Seelisch <seelisch@…>, 13 years ago
new gfan lib version by Anders Jensen git-svn-id: file:///usr/local/Singular/svn/trunk@13668 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 15.2 KB
Line 
1/*
2 * gfanlib_symmetry.cpp
3 *
4 *  Created on: Oct 22, 2010
5 *      Author: anders
6 */
7
8#include "gfanlib_symmetry.h"
9#include <map>
10
11using namespace std;
12
13namespace gfan{
14class Trie
15{
16        class TrieNode
17        {
18                typedef map<int,class TrieNode> Map;
19                Map m;
20        public:
21                TrieNode()
22                {
23
24                }
25                TrieNode(IntVector const &v, int i)
26                {
27                        if(i<v.size())
28                        m[v[i]]=TrieNode(v,i+1);
29                }
30                int stabilizerSize(ZVector const &v, int i)const
31                {
32                  int ret=0;
33                  if(i==v.size())return 1;
34                  for(Map::const_iterator j=m.begin();j!=m.end();j++)
35                    {
36                      if(v[i]==v[j->first])
37                        ret+=j->second.stabilizerSize(v,i+1);
38                    }
39                  return ret;
40                }
41                void search(ZVector const &v, ZVector  &building, Permutation &tempPerm, Permutation &ret, ZVector &optimal, int i, bool &isImproving)const
42                {
43                        if(i==v.size()){ret=tempPerm;optimal=building;isImproving=false;return;}
44                        if(isImproving)
45                                building[i]=-0x7fffffff;
46                        else
47                                building[i]=optimal[i];
48                        for(Map::const_iterator j=m.begin();j!=m.end();j++)
49                                if(building[i]<v[j->first])
50                                {
51                                        isImproving=true;
52                                        building[i]=v[j->first];
53                                }
54                                for(Map::const_iterator j=m.begin();j!=m.end();j++)
55                                if(v[j->first]==building[i])
56                                {
57                                        tempPerm[i]=j->first;
58                                        j->second.search(v,building,tempPerm,ret,optimal,i+1,isImproving);
59                                }
60                }
61                void searchStabalizer(ZVector const &v, ZVector  &building, Permutation &tempPerm, Permutation &ret, ZVector &optimal, int i, bool &isImproving, ZVector const &toBeFixed)const
62                {
63                        if(i==v.size())
64                                if(!(tempPerm.apply(v)<optimal))
65                                        {
66                                        ret=tempPerm;
67                                        optimal=tempPerm.apply(v);
68                                        return;
69                                        }
70                                for(Map::const_iterator j=m.begin();j!=m.end();j++)
71                                        if(toBeFixed[i]==toBeFixed[j->first])
72                                {
73                                        tempPerm[i]=j->first;
74                                        j->second.searchStabalizer(v,building,tempPerm,ret,optimal,i+1,isImproving,toBeFixed);
75                                }
76                }
77/* this code contains mistakes          void searchStabalizer(IntegerVector const &v, IntegerVector  &building, IntegerVector &tempPerm, IntegerVector &ret, IntegerVector &optimal, int i, bool &isImproving, IntegerVector const &toBeFixed)const
78                {
79                        if(i==v.size()){ret=tempPerm;optimal=building;isImproving=false;debug<<"DEEP";return;}
80                        if(isImproving)
81                                building[i]=-0x7fffffff;
82                        else
83                                building[i]=optimal[i];
84                        for(Map::const_iterator j=m.begin();j!=m.end();j++)
85                                if(toBeFixed[i]==toBeFixed[j->first])
86                                if(v[j->first]>building[i])
87                                {
88                                        isImproving=true;
89                                        building[i]=v[j->first];
90                                }
91                                for(Map::const_iterator j=m.begin();j!=m.end();j++)
92                                        if(toBeFixed[i]==toBeFixed[j->first])
93                                if(v[j->first]==building[i])
94                                {
95                                        debug.printInteger(i);debug<<":";
96                                        debug.printInteger(j->first);debug<<" ";
97                                        tempPerm[i]=j->first;
98                                        j->second.searchStabalizer(v,building,tempPerm,ret,optimal,i+1,isImproving,toBeFixed);
99                                }
100                }*/
101        //      void doubleSearch();
102                void insert(Permutation const &v, int i)
103                {
104                        if(i==v.size())return;
105                        if(m.count(v[i]))
106                                m[v[i]].insert(v,i+1);
107                        else
108                                m[v[i]]=                TrieNode(v,i+1);
109                }
110/*                void print(int i, int n)const
111                {
112                        if(i==n)return;
113                        for(Map::const_iterator j=m.begin();j!=m.end();j++)
114                        {
115                                {for(int j=0;j<2*i;j++)debug<<" ";}
116                                debug.printInteger(j->first);
117                                debug<<"\n";
118                                j->second.print(i+1,n);
119                        }
120                        }*/
121                int size(int i,int n)const
122                {
123                        if(i==n)return 1;
124                        int ret=0;
125                        for(Map::const_iterator j=m.begin();j!=m.end();j++)
126                                ret+=j->second.size(i+1,n);
127                        return ret;
128                }
129        };
130public:
131        TrieNode theTree;
132        int n;
133        Trie(int n_):
134                n(n_),
135                theTree(Permutation(n_),0)
136        {
137        }
138        int size()const
139        {
140                return theTree.size(0,n);
141        }
142        void insert(Permutation const &v)
143        {
144                theTree.insert(v,0);
145//              debug<<v;
146//              theTree.print(0,v.size());
147
148//              debug<<"---------------------------------------------\n";
149        }
150        /**
151         * returns the sigma from the set with sigma(v) maximal in the lexicographic ordering.
152         */
153        Permutation search(ZVector const &v)
154        {
155                Permutation tempPerm(v.size());
156                Permutation ret(v.size());
157                ZVector building(v.size());
158                ZVector optimal=v;//the identity is always in the trie
159                bool isImproving=true;
160                theTree.search(v,building,tempPerm,ret,optimal,0,isImproving);
161                return ret;
162        }
163        Permutation searchStabalizer(ZVector const &v, ZVector const &toBeFixed)
164        {
165                Permutation tempPerm(v.size());
166                Permutation ret(v.size());
167                ZVector building(v.size());
168                ZVector optimal=v;//the identity is always in the trie
169                bool isImproving=true;
170                theTree.searchStabalizer(v,building,tempPerm,ret,optimal,0,isImproving,toBeFixed);
171                return ret;
172        }
173        int stabilizerSize(ZVector const &v)const
174        {
175          return theTree.stabilizerSize(v,0);
176        }
177};
178
179
180
181
182
183
184
185/*IntegerVector SymmetryGroup::identity(int n)
186{
187  IntegerVector v(n);
188  for(int i=0;i<n;i++)v[i]=i;
189
190  return v;
191}
192*/
193
194Permutation Permutation::inverse()const
195{
196  return applyInverse(Permutation(size()));
197}
198
199
200SymmetryGroup::SymmetryGroup(int n):
201//  byteTable(0),
202  trie(0)
203{
204  elements.insert(Permutation(n));
205}
206
207
208int SymmetryGroup::sizeOfBaseSet()const
209{
210  assert(!elements.empty());
211  return elements.begin()->size();
212}
213
214void SymmetryGroup::computeClosure(Permutation const &v) //does this work??
215{
216  ElementContainer newOnes;
217
218  newOnes.insert(v);
219
220  while(!newOnes.empty())
221    {
222      static int i;
223      i++;
224
225      Permutation v=*newOnes.begin();
226      for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)
227        {
228          {
229            Permutation n=i->apply(v);
230            if(0==elements.count(n))
231              newOnes.insert(n);
232          }
233          {
234            Permutation n=v.apply(*i);
235            if(0==elements.count(n))
236              newOnes.insert(n);
237          }
238        }
239      newOnes.erase(v);
240      elements.insert(v);
241    }
242}
243
244
245/*
246
247void SymmetryGroup::computeClosure(IntegerVectorList const &l)
248{
249  //  for(IntegerVectorList::const_iterator i=l.begin();i!=l.end();i++)
250  //  computeClosure(*i);
251
252  bool growing=true;
253  while(growing)
254    {
255      growing=false;
256      for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)
257        {
258          for(IntegerVectorList::const_iterator j=l.begin();j!=l.end();j++)
259            {
260              {
261                IntegerVector n(compose(*i,*j));
262                growing|=(0==elements.count(n));
263                elements.insert(n);
264              }
265              {
266                IntegerVector n(compose(*i,*j));
267                growing|=(0==elements.count(n));
268                elements.insert(n);
269              }
270            }
271        }
272    }
273}
274*/
275
276/*
277void SymmetryGroup::print(FILE *f)
278{
279  AsciiPrinter P(f);
280  P.printString("Printing SymmetryGroup\n");
281  IntegerVectorList l;
282  for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)
283    {
284      //      P.printVector(*i);
285      //      P.printNewLine();
286      l.push_back(*i);
287    }
288  P.printVectorList(l);
289  fprintf(f,"Group order=%i\n",elements.size());
290  P.printString("Done printing SymmetryGroup.\n");
291}
292*/
293
294
295Permutation Permutation::apply(Permutation const &b)const
296{
297  IntVector ret(size());
298  assert(size()==b.size());
299  for(int i=0;i<size();i++)ret[i]=b[(*this)[i]];
300  return Permutation(ret);
301}
302
303Permutation Permutation::applyInverse(Permutation const &b)const
304{
305  IntVector ret(size());
306  assert(size()==b.size());
307  for(int i=0;i<size();i++)ret[(*this)[i]]=b[i];
308  return Permutation(ret);
309}
310
311IntVector Permutation::apply(IntVector const &v)const
312{
313  IntVector ret(size());
314  assert(size()==v.size());
315  for(int i=0;i<size();i++)ret[i]=v[(*this)[i]];
316  return ret;
317}
318
319ZVector Permutation::apply(ZVector const &v)const
320{
321  ZVector ret(size());
322  assert(size()==v.size());
323  for(int i=0;i<size();i++)ret[i]=v[(*this)[i]];
324  return ret;
325}
326
327
328ZVector Permutation::applyInverse(ZVector const &v)const
329{
330  ZVector ret(size());
331  assert(size()==v.size());
332  for(int i=0;i<size();i++)ret[(*this)[i]]=v[i];
333  return ret;
334}
335//ZVector apply(ZVector const &v)const;
336//ZMatrix apply(ZMatrix const &m)const;
337//IntVector applyInverse(IntVector const &v)const;
338//ZMatrix applyInverse(ZMatrix const &m)const;
339
340ZVector SymmetryGroup::orbitRepresentative(ZVector const &v, Permutation *usedPermutation)const
341{
342  if(trie){
343    if(usedPermutation)
344      {
345        *usedPermutation=trie->search(v);
346        return usedPermutation->apply(v);
347      }
348    return trie->search(v).apply(v);
349  }
350  ZVector ret=v;
351  ElementContainer::const_iterator usedPerm;
352  for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)
353    {
354      ZVector q=i->apply(v);
355      if(! (q<ret))//negation to make sure that usedPerm is set
356        {
357          usedPerm=i;
358          ret=q;
359        }
360    }
361
362  if(usedPermutation)*usedPermutation=*usedPerm;
363
364  if(trie)
365  {
366//        debug<<"Input"<<v<<"\n";
367//        debug<<"Bruteforce"<<ret<<"\n";
368          Permutation triePerm=trie->search(v);
369//        debug<<"Trie"<<compose(triePerm,v)<<"\n";
370          assert((triePerm.apply(v)-ret).isZero());
371  }
372
373  return ret;
374}
375#if 0
376IntegerVector SymmetryGroup::orbitRepresentativeFixing(IntegerVector const &v, IntegerVector const &fixed)const
377{
378        if(trie){
379                return compose(trie->searchStabalizer(v,fixed),v);
380        }
381  IntegerVector ret=v;
382
383  for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)
384    if(compose(*i,fixed)==fixed)
385      {
386        IntegerVector q=compose(*i,v);
387        if(ret<q)ret=q;
388      }
389        if(trie){
390                IntegerVector temp=compose(trie->searchStabalizer(v,fixed),v);
391//              debug<<"Input"<<v;
392//              debug<<"Brute"<<ret;
393//              debug<<"Quick"<<temp;
394                assert((temp-ret).isZero());
395//              return compose(trie->searchStabalizer(v,fixed),v);
396        }
397  return ret;
398}
399#endif
400
401bool Permutation::isPermutation(IntVector const &a)
402{
403  int n=a.size();
404  IntVector temp(n);
405  for(int i=0;i<n;i++)temp[i]=-1;
406  for(int i=0;i<n;i++)
407    {
408      if(a[i]<0 || a[i]>=n)return false;
409      temp[i]=i;
410    }
411  for(int i=0;i<n;i++)if(temp[i]<0)return false;
412  return true;
413}
414
415
416int SymmetryGroup::orbitSize(ZVector const &stable)const
417{
418  int groupSize=elements.size();
419
420  int n=stable.size();
421  int numFixed=0;
422
423  if(trie)
424    {
425      numFixed=trie->stabilizerSize(stable);
426    }
427  else
428    {
429      for(SymmetryGroup::ElementContainer::const_iterator j=elements.begin();j!=elements.end();j++)
430        {
431          bool doesFix=true;
432
433          for(int i=0;i<n;i++)
434            if(stable[i]!=stable[(*j)[i]])
435              {
436                doesFix=false;
437                break;
438              }
439          if(doesFix)numFixed++;
440        }
441    }
442  return groupSize/numFixed;
443}
444
445
446ZVector Permutation::fundamentalDomainInequality()const
447{
448  for(int i=0;i<size();i++)
449    if((*this)[i]!=i)
450      return ZVector::standardVector(size(),i)-ZVector::standardVector(size(),(*this)[i]);
451  return ZVector(size());
452}
453
454
455ZMatrix SymmetryGroup::fundamentalDomainInequalities()const
456{
457  set<ZVector> ret2;
458  for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)
459    ret2.insert(i->fundamentalDomainInequality());
460
461  ZMatrix ret(0,sizeOfBaseSet());
462  for(set<ZVector>::const_iterator i=ret2.begin();i!=ret2.end();i++)
463    if(!i->isZero())ret.appendRow(*i);
464
465  return ret;
466}
467
468
469void SymmetryGroup::createTrie()
470{
471  trie=new Trie(sizeOfBaseSet());
472  for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)
473    trie->insert(*i);
474}
475
476
477
478bool SymmetryGroup::isTrivial()const
479{
480  ElementContainer::const_iterator i=elements.begin();
481  assert(i!=elements.end());
482  i++;
483  return i==elements.end();
484}
485#if 0
486static int mergeSortRek(IntegerVector &v, int begin, int end, IntegerVector &temp)
487{
488  if(end-begin<2)return 0;
489  int med=(begin+end)>>1;
490  int nswaps=mergeSortRek(v,begin,med,temp);
491  nswaps+=mergeSortRek(v,med,end,temp);
492
493  {
494    int Astart=begin;
495    int Alength=med-begin;
496    int Bstart=med;
497    int Blength=end-med;
498    int nextFree=begin;
499    while(nextFree!=end)
500      {
501//        debug<<"Astart:"<<Astart<<"Alength:"<<Alength<<"Bstart:"<<Bstart<<"Blength:"<<Blength<<"nextFree:"<<nextFree<<"\n";
502        if(Blength==0 || (Alength!=0 && v[Astart]<v[Bstart]))
503          {
504            temp[nextFree++]=v[Astart++];
505            Alength--;
506          }
507        else
508          {
509            temp[nextFree++]=v[Bstart++];
510            nswaps+=Alength;
511            Blength--;
512          }
513      }
514    for(int i=begin;i!=end;i++)v[i]=temp[i];
515  }
516//debug<<"return\n";
517  return nswaps;
518}
519
520int mergeSort(IntegerVector &v)
521{
522  IntegerVector temp(v.size());
523  return mergeSortRek(v,0,v.size(),temp);
524}
525#endif
526}
Note: See TracBrowser for help on using the repository browser.