source: git/gfanlib/gfanlib_symmetry.cpp @ 5417ff

spielwiese
Last change on this file since 5417ff was 43317d, checked in by Frank Seelisch <seelisch@…>, 13 years ago
creation of fans and insertion of cones therein git-svn-id: file:///usr/local/Singular/svn/trunk@13679 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 14.8 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
214IntMatrix SymmetryGroup::getGenerators()const
215{
216  IntMatrix ret(0,this->sizeOfBaseSet());
217  for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)ret.appendRow(i->toIntVector());
218  return ret;
219}
220
221void SymmetryGroup::computeClosure(Permutation const &v) //does this work??
222{
223  ElementContainer newOnes;
224
225  newOnes.insert(v);
226
227  while(!newOnes.empty())
228    {
229      Permutation v=*newOnes.begin();
230      for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)
231        {
232          {
233            Permutation n=i->apply(v);
234            if(0==elements.count(n))
235              newOnes.insert(n);
236          }
237          {
238            Permutation n=v.apply(v);
239            if(0==elements.count(n))
240              newOnes.insert(n);
241          }
242        }
243      newOnes.erase(v);
244      elements.insert(v);
245    }
246}
247
248
249void SymmetryGroup::computeClosure(IntMatrix const &l)
250{
251  for(int i=0;i<l.getHeight();i++)computeClosure(Permutation(l[i]));
252}
253
254
255/*
256void SymmetryGroup::print(FILE *f)
257{
258  AsciiPrinter P(f);
259  P.printString("Printing SymmetryGroup\n");
260  IntegerVectorList l;
261  for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)
262    {
263      //      P.printVector(*i);
264      //      P.printNewLine();
265      l.push_back(*i);
266    }
267  P.printVectorList(l);
268  fprintf(f,"Group order=%i\n",elements.size());
269  P.printString("Done printing SymmetryGroup.\n");
270}
271*/
272
273
274Permutation Permutation::apply(Permutation const &b)const
275{
276  IntVector ret(size());
277  assert(size()==b.size());
278  for(int i=0;i<size();i++)ret[i]=b[(*this)[i]];
279  return Permutation(ret);
280}
281
282Permutation Permutation::applyInverse(Permutation const &b)const
283{
284  IntVector ret(size());
285  assert(size()==b.size());
286  for(int i=0;i<size();i++)ret[(*this)[i]]=b[i];
287  return Permutation(ret);
288}
289
290IntVector Permutation::apply(IntVector const &v)const
291{
292  IntVector ret(size());
293  assert(size()==v.size());
294  for(int i=0;i<size();i++)ret[i]=v[(*this)[i]];
295  return ret;
296}
297
298ZVector Permutation::apply(ZVector const &v)const
299{
300  ZVector ret(size());
301  assert(size()==v.size());
302  for(int i=0;i<size();i++)ret[i]=v[(*this)[i]];
303  return ret;
304}
305
306
307ZVector Permutation::applyInverse(ZVector const &v)const
308{
309  ZVector ret(size());
310  assert(size()==v.size());
311  for(int i=0;i<size();i++)ret[(*this)[i]]=v[i];
312  return ret;
313}
314//ZVector apply(ZVector const &v)const;
315//ZMatrix apply(ZMatrix const &m)const;
316//IntVector applyInverse(IntVector const &v)const;
317//ZMatrix applyInverse(ZMatrix const &m)const;
318
319ZVector SymmetryGroup::orbitRepresentative(ZVector const &v, Permutation *usedPermutation)const
320{
321  if(trie){
322    if(usedPermutation)
323      {
324        *usedPermutation=trie->search(v);
325        return usedPermutation->apply(v);
326      }
327    return trie->search(v).apply(v);
328  }
329  ZVector ret=v;
330  ElementContainer::const_iterator usedPerm;
331  for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)
332    {
333      ZVector q=i->apply(v);
334      if(! (q<ret))//negation to make sure that usedPerm is set
335        {
336          usedPerm=i;
337          ret=q;
338        }
339    }
340
341  if(usedPermutation)*usedPermutation=*usedPerm;
342
343  if(trie)
344  {
345//        debug<<"Input"<<v<<"\n";
346//        debug<<"Bruteforce"<<ret<<"\n";
347          Permutation triePerm=trie->search(v);
348//        debug<<"Trie"<<compose(triePerm,v)<<"\n";
349          assert((triePerm.apply(v)-ret).isZero());
350  }
351
352  return ret;
353}
354#if 0
355IntegerVector SymmetryGroup::orbitRepresentativeFixing(IntegerVector const &v, IntegerVector const &fixed)const
356{
357        if(trie){
358                return compose(trie->searchStabalizer(v,fixed),v);
359        }
360  IntegerVector ret=v;
361
362  for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)
363    if(compose(*i,fixed)==fixed)
364      {
365        IntegerVector q=compose(*i,v);
366        if(ret<q)ret=q;
367      }
368        if(trie){
369                IntegerVector temp=compose(trie->searchStabalizer(v,fixed),v);
370//              debug<<"Input"<<v;
371//              debug<<"Brute"<<ret;
372//              debug<<"Quick"<<temp;
373                assert((temp-ret).isZero());
374//              return compose(trie->searchStabalizer(v,fixed),v);
375        }
376  return ret;
377}
378#endif
379
380bool Permutation::isPermutation(IntVector const &a)
381{
382  int n=a.size();
383  IntVector temp(n);
384  for(int i=0;i<n;i++)temp[i]=-1;
385  for(int i=0;i<n;i++)
386    {
387      if(a[i]<0 || a[i]>=n)return false;
388      temp[i]=i;
389    }
390  for(int i=0;i<n;i++)if(temp[i]<0)return false;
391  return true;
392}
393
394
395bool Permutation::arePermutations(IntMatrix const &m)
396{
397  for(int i=0;i<m.getHeight();i++)if(!isPermutation(m[i]))return false;
398  return true;
399}
400
401int SymmetryGroup::orbitSize(ZVector const &stable)const
402{
403  int groupSize=elements.size();
404
405  int n=stable.size();
406  int numFixed=0;
407
408  if(trie)
409    {
410      numFixed=trie->stabilizerSize(stable);
411    }
412  else
413    {
414      for(SymmetryGroup::ElementContainer::const_iterator j=elements.begin();j!=elements.end();j++)
415        {
416          bool doesFix=true;
417
418          for(int i=0;i<n;i++)
419            if(stable[i]!=stable[(*j)[i]])
420              {
421                doesFix=false;
422                break;
423              }
424          if(doesFix)numFixed++;
425        }
426    }
427  return groupSize/numFixed;
428}
429
430
431ZVector Permutation::fundamentalDomainInequality()const
432{
433  for(int i=0;i<size();i++)
434    if((*this)[i]!=i)
435      return ZVector::standardVector(size(),i)-ZVector::standardVector(size(),(*this)[i]);
436  return ZVector(size());
437}
438
439
440ZMatrix SymmetryGroup::fundamentalDomainInequalities()const
441{
442  set<ZVector> ret2;
443  for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)
444    ret2.insert(i->fundamentalDomainInequality());
445
446  ZMatrix ret(0,sizeOfBaseSet());
447  for(set<ZVector>::const_iterator i=ret2.begin();i!=ret2.end();i++)
448    if(!i->isZero())ret.appendRow(*i);
449
450  return ret;
451}
452
453
454void SymmetryGroup::createTrie()
455{
456  trie=new Trie(sizeOfBaseSet());
457  for(ElementContainer::const_iterator i=elements.begin();i!=elements.end();i++)
458    trie->insert(*i);
459}
460
461
462
463bool SymmetryGroup::isTrivial()const
464{
465  ElementContainer::const_iterator i=elements.begin();
466  assert(i!=elements.end());
467  i++;
468  return i==elements.end();
469}
470#if 0
471static int mergeSortRek(IntegerVector &v, int begin, int end, IntegerVector &temp)
472{
473  if(end-begin<2)return 0;
474  int med=(begin+end)>>1;
475  int nswaps=mergeSortRek(v,begin,med,temp);
476  nswaps+=mergeSortRek(v,med,end,temp);
477
478  {
479    int Astart=begin;
480    int Alength=med-begin;
481    int Bstart=med;
482    int Blength=end-med;
483    int nextFree=begin;
484    while(nextFree!=end)
485      {
486//        debug<<"Astart:"<<Astart<<"Alength:"<<Alength<<"Bstart:"<<Bstart<<"Blength:"<<Blength<<"nextFree:"<<nextFree<<"\n";
487        if(Blength==0 || (Alength!=0 && v[Astart]<v[Bstart]))
488          {
489            temp[nextFree++]=v[Astart++];
490            Alength--;
491          }
492        else
493          {
494            temp[nextFree++]=v[Bstart++];
495            nswaps+=Alength;
496            Blength--;
497          }
498      }
499    for(int i=begin;i!=end;i++)v[i]=temp[i];
500  }
501//debug<<"return\n";
502  return nswaps;
503}
504
505int mergeSort(IntegerVector &v)
506{
507  IntegerVector temp(v.size());
508  return mergeSortRek(v,0,v.size(),temp);
509}
510#endif
511}
Note: See TracBrowser for help on using the repository browser.