Cadabra Computer algebra system for field theory problems
Combinatorics.hh
Go to the documentation of this file.
1 /*
2
3  Cadabra: a field-theory motivated computer algebra system.
4  Copyright (C) 2001-2014 Kasper Peeters <kasper.peeters@phi-sci.com>
5
6 This program is free software: you can redistribute it and/or
7  modify it under the terms of the GNU General Public License as
8 published by the Free Software Foundation, either version 3 of the
10
11 This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17  along with this program. If not, see <http://www.gnu.org/licenses/>.
18
19 */
20
21 /*
22  length vector
23  normal combinations: one element, value=total length.
24  normal permutations: n elements, each equal to 1.
25
26 */
27
28 #pragma once
29
30 #include <vector>
31 #include <cassert>
32 #include <algorithm>
33 #include <iomanip>
34 #include <iostream>
35 #include <map>
36
37 namespace combin {
38
39  typedef std::vector<unsigned int> range_t;
40  typedef std::vector<range_t> range_vector_t;
41  typedef std::vector<int> weights_t;
42
43  unsigned long factorial(unsigned int x);
45  long vector_sum(const std::vector<int>&);
47  unsigned long vector_prod(const std::vector<unsigned int>&);
49  unsigned long vector_prod_fact(const std::vector<unsigned int>&);
50
51  bool operator==(const std::vector<unsigned int>&, const std::vector<unsigned int>&);
52
54  long hash(const std::vector<unsigned int>&);
55
56  template<class T>
58  public:
60  combinations_base(const std::vector<T>&);
61  virtual ~combinations_base();
62
63  void permute(long start=-1, long end=-1);
64  virtual void clear();
65  virtual void clear_results();
66  unsigned int sum_of_sublengths() const;
67  void set_unit_sublengths();
68  unsigned int multiplier(const std::vector<T>&) const;
69  unsigned int total_permutations() const; // including the ones not stored
70
72
73  unsigned int block_length;
74  std::vector<unsigned int> sublengths;
76  std::vector<T> original;
78  std::vector<weights_t> weights;
79  std::vector<int> max_weights;
80  std::vector<weight_cond> weight_conditions;
81  unsigned int sub_problem_blocksize; // when non-zero, do permutations within
82
83  protected:
84  virtual void vector_generated(const std::vector<unsigned int>&)=0;
85  virtual bool entry_accepted(unsigned int current) const;
86  std::vector<unsigned int> temparr;
87
89  std::vector<int> current_weight;
90  private:
91  bool is_allowed_by_weight_constraints(unsigned int i);
92  bool final_weight_constraints_check() const;
93  void update_weights(unsigned int i);
94  void restore_weights(unsigned int i);
95  void nextstep(unsigned int current, unsigned int fromalgehad, unsigned int groupindex,
97
98  };
99
100  template<class T>
101  class combinations : public combinations_base<T> {
102  public:
103  typedef typename std::vector<std::vector<T> > permuted_sets_t;
104  typedef typename permuted_sets_t::const_iterator const_iterator;
105
107  combinations(const std::vector<T>&);
108  virtual ~combinations();
109
110  virtual void clear();
111  virtual void clear_results();
112  const std::vector<T>& operator[](unsigned int) const;
113  int ordersign(unsigned int) const;
114  unsigned int size() const;
115  unsigned int multiplier(unsigned int) const;
116
117  protected:
118  virtual void vector_generated(const std::vector<unsigned int>&);
119
120  private:
122  };
123
124  template<class T>
125  class symmetriser;
126
127  template<class T>
128  class symm_helper : public combinations_base<T> {
129  public:
131  virtual void clear();
132
134  protected:
135  bool first_one;
137  virtual void vector_generated(const std::vector<unsigned int>&);
138  };
139
140  template<class T>
141  class symm_val_helper : public combinations_base<T> {
142  public:
144  virtual void clear();
145
147  protected:
148  bool first_one;
150  virtual void vector_generated(const std::vector<unsigned int>&);
151  };
152
153  template<class T>
154  class symmetriser {
155  public:
157  void apply_symmetry(long start=-1, long end=-1);
158
159  std::vector<T> original;
160  unsigned int block_length;
161  std::vector<unsigned int> permute_blocks; // offset in unit elements! (not in blocks)
162  std::vector<T> value_permute;
164
165  std::vector<unsigned int> sublengths; // refers to position within permute_blocks
166  range_vector_t input_asym; // as in combinations_base
167  range_vector_t sublengths_scattered; // sublengths, in original, but not connected.
168
171  range_t values_to_locations(const std::vector<T>& values) const;
172
173  const std::vector<T>& operator[](unsigned int) const;
174  int signature(unsigned int) const;
175  void set_multiplicity(unsigned int pos, int val);
176  unsigned int size() const;
177  void clear();
179  void collect();
181
182  friend class symm_helper<T>;
183  friend class symm_val_helper<T>;
184  private:
187  unsigned int current_;
188  std::vector<std::vector<T> > originals;
189  std::vector<int> multiplicity;
190  };
191
193  const range_vector_t& indv,
194  range_vector_t& target);
195
196  template<class iterator1, class iterator2>
197  int ordersign(iterator1 b1, iterator1 e1, iterator2 b2, iterator2 e2, int stepsize=1);
198
199  template<class iterator1>
200  int ordersign(iterator1 b1, iterator1 e1);
201
202  template<class T>
203  T fact(T x);
204
205  template<class T>
206  std::ostream& operator<<(std::ostream& str, const symmetriser<T>& sym);
207
208
209  /*
210  I assume PI consists of the integers 1 to N.
211  It can be done with O(N) comparisons and transpositions of integers
212  in the list.
213
214  sign:= 1;
215  for i from 1 to N do
216  while PI[i] <> i do
217  interchange PI[i] and PI[PI[i]];
218  sign:= -sign
219  od
220  od
221  */
222
223  template<class iterator1, class iterator2>
224  int ordersign(iterator1 b1, iterator1 e1, iterator2 b2, iterator2 e2, int stepsize)
225  {
226  int sign=1;
227  std::vector<bool> crossedoff(std::distance(b1,e1),false);
228  while(b1!=e1) {
229  int otherpos=0;
230  iterator2 it=b2;
231  while(it!=e2) {
232  if( (*it)==(*b1) && crossedoff[otherpos]==false) {
233  crossedoff[otherpos]=true;
234  break;
235  }
236  else {
237  if(!crossedoff[otherpos])
238  sign=-sign;
239  }
240  it+=stepsize;
241  ++otherpos;
242  }
243  b1+=stepsize;
244  }
245  return sign;
246  }
247
248  //template<class iterator1, class iterator2, class comparator>
249  //int ordersign(iterator1 b1, iterator1 e1, iterator2 b2, iterator2 e2, comparator cmp, int stepsize)
250  // {
251  // int sign=1;
252  // std::vector<bool> crossedoff(std::distance(b1,e1),false);
253  // while(b1!=e1) {
254  // int otherpos=0;
255  // iterator2 it=b2;
256  // while(it!=e2) {
257  // if(cmp((*it), (*b1)) && crossedoff[otherpos]==false) {
258  // crossedoff[otherpos]=true;
259  // break;
260  // }
261  // else {
262  // if(!crossedoff[otherpos])
263  // sign=-sign;
264  // }
265  // it+=stepsize;
266  // ++otherpos;
267  // }
268  // b1+=stepsize;
269  // }
270  // return sign;
271  // }
272
273  template<class iterator1>
274  int ordersign(iterator1 b1, iterator1 e1)
275  {
276  std::vector<unsigned int> fil;
277  for(int k=0; k<distance(b1,e1); ++k)
278  fil.push_back(k);
279  return ordersign(fil.begin(), fil.end(), b1, e1);
280  }
281
282  template<class T>
283  T fact(T x)
284  {
285  T ret=1;
286  assert(x>=0);
287  while(x!=0) {
288  ret*=x--;
289  }
290  return ret;
291  }
292
293
294  // Implementations
295
296  template<class T>
298  : block_length(1), multiple_pick(false), sub_problem_blocksize(0)
299  {
300  }
301
302  template<class T>
303  combinations_base<T>::combinations_base(const std::vector<T>& oa)
304  : block_length(1), original(oa), multiple_pick(false), sub_problem_blocksize(0)
305  {
306  }
307
308  template<class T>
310  : combinations_base<T>()
311  {
312  }
313
314  template<class T>
315  combinations<T>::combinations(const std::vector<T>& oa)
316  : combinations_base<T>(oa)
317  {
318  }
319
320  template<class T>
322  {
323  }
324
325  template<class T>
327  {
328  }
329
330  template<class T>
331  void combinations<T>::vector_generated(const std::vector<unsigned int>& toadd)
332  {
333  ++this->vector_generated_called_;
334  if((this->start_==-1 || this->vector_generated_called_ >= this->start_) &&
335  (this->end_==-1 || this->vector_generated_called_ < this->end_)) {
337  for(unsigned int i=0; i<toadd.size(); ++i)
338  for(unsigned int bl=0; bl<this->block_length; ++bl)
340  storage.push_back(newone);
341  }
342  }
343
344  template<class T>
345  bool combinations_base<T>::entry_accepted(unsigned int) const
346  {
347  return true;
348  }
349
350  template<class T>
352  {
353  start_=start;
354  end_=end;
355  vector_generated_called_=-1;
356
357  // Initialise weight handling.
358  current_weight.clear();
359  current_weight.resize(weights.size(), 0);
360  for(unsigned int i=0; i<weights.size(); ++i)
361  assert(weights[i].size() == original.size()/block_length);
362  if(weights.size()>0) {
363  if(weight_conditions.size()==0)
364  weight_conditions.resize(weights.size(), weight_equals);
365  else assert(weight_conditions.size()==weights.size());
366  }
367  else assert(weight_conditions.size()==0);
368
369  // Sublength handling.
370  assert(sublengths.size()!=0);
371  unsigned int len=sum_of_sublengths();
372
373  // Consistency checks.
374  assert(original.size()%block_length==0);
375  if(!multiple_pick)
376  assert(len*block_length<=original.size());
377
378  for(unsigned int i=0; i<this->input_asym.size(); ++i)
379  std::sort(this->input_asym[i].begin(), this->input_asym[i].end());
380
381  temparr=std::vector<unsigned int>(len/* *block_length*/);
384  }
385
386  template<class T>
388  {
389  block_length=1;
390  sublengths.clear();
391  this->input_asym.clear();
392  original.clear();
393  weights.clear();
394  max_weights.clear();
395  weight_conditions.clear();
396  sub_problem_blocksize=0;
397  temparr.clear();
398  current_weight.clear();
399  }
400
401  template<class T>
403  {
404  temparr.clear();
405  }
406
407  template<class T>
409  {
410  storage.clear();
412  }
413
414  template<class T>
416  {
417  storage.clear();
419  }
420
421  template<class T>
422  const std::vector<T>& combinations<T>::operator[](unsigned int i) const
423  {
424  assert(i<storage.size());
425  return storage[i];
426  }
427
428  template<class T>
429  unsigned int combinations<T>::size() const
430  {
431  return storage.size();
432  }
433
434  template<class T>
436  {
437  unsigned int ret=0;
438  for(unsigned int i=0; i<sublengths.size(); ++i)
439  ret+=sublengths[i];
440  return ret;
441  }
442
443  template<class T>
445  {
446  return vector_generated_called_+1;
447  }
448
449  template<class T>
451  {
452  sublengths.clear();
453  for(unsigned int i=0; i<original.size()/block_length; ++i)
454  sublengths.push_back(1);
455  }
456
457  template<class T>
458  int combinations<T>::ordersign(unsigned int num) const
459  {
460  assert(num<storage.size());
461  return combin::ordersign(storage[0].begin(), storage[0].end(),
462  storage[num].begin(), storage[num].end(), this->block_length);
463  }
464
465  template<class T>
466  unsigned int combinations<T>::multiplier(unsigned int num) const
467  {
468  return combinations_base<T>::multiplier(this->storage[num]);
469  }
470
471  template<class T>
472  unsigned int combinations_base<T>::multiplier(const std::vector<T>& stor) const
473  {
474  unsigned long numerator=1;
475  for(unsigned int i=0; i<this->input_asym.size(); ++i)
476  numerator*=fact(this->input_asym[i].size());
477
478  unsigned long denominator=1;
479  for(unsigned int i=0; i<this->input_asym.size(); ++i) {
480  // for each input asym, and for each output asym, count
481  // the number of overlap elements.
482  unsigned int current=0;
483  for(unsigned int k=0; k<this->sublengths.size(); ++k) {
484  if(this->sublengths[k]>1) {
485  unsigned int overlap=0;
486  for(unsigned int slc=0; slc<this->sublengths[k]; ++slc) {
487  for(unsigned int j=0; j<this->input_asym[i].size(); ++j) {
488  unsigned int index=0;
489  while(!(stor[current]==this->original[index]))
490  ++index;
491  if(index==this->input_asym[i][j])
492  ++overlap;
493  }
494  ++current;
495  }
496  if(overlap>0)
497  denominator*=fact(overlap);
498  // FIXME: for each overlap thus found, divide out by a factor
499  // due to the fact that output asym ranges can overlap.
500  // well, that's not right either.
501  }
502  else ++current;
503  }
504  }
505
506  return numerator/denominator;
507  }
508
509  template<class T>
511  {
512  if(weights.size()==0) return true;
513  for(unsigned int cn=0; cn<current_weight.size(); ++cn) {
514  if(weight_conditions[cn]==weight_less)
515  if(current_weight[cn]+weights[cn][i] >= max_weights[cn])
516  return false;
517  }
518  return true;
519  }
520
521  template<class T>
523  {
524  for(unsigned int cn=0; cn<current_weight.size(); ++cn) {
525  switch(weight_conditions[cn]) {
526  case weight_equals:
527  if(current_weight[cn]!=max_weights[cn])
528  return false;
529  break;
530  case weight_less:
531  break;
532  case weight_greater:
533  if(current_weight[cn]<=max_weights[cn])
534  return false;
535  break;
536  }
537  }
538  return true;
539  }
540
541  template<class T>
543  {
544  if(weights.size()==0) return;
545  for(unsigned int cn=0; cn<current_weight.size(); ++cn)
546  current_weight[cn]+=weights[cn][i];
547  }
548
549  template<class T>
551  {
552  if(weights.size()==0) return;
553  for(unsigned int cn=0; cn<current_weight.size(); ++cn)
554  current_weight[cn]-=weights[cn][i];
555  }
556
557  template<class T>
558  void combinations_base<T>::nextstep(unsigned int current, unsigned int lowest_in_group, unsigned int groupindex,
560  {
561  unsigned int grouplen=0;
562  for(unsigned int i=0; i<=groupindex; ++i)
563  grouplen+=sublengths[i];
564  if(current==grouplen) { // group is filled
565  ++groupindex;
566  if(groupindex==sublengths.size()) {
567  if(final_weight_constraints_check())
568  vector_generated(temparr);
569  return;
570  }
571  lowest_in_group=0;
572  }
573
574  unsigned int starti=0, endi=original.size()/block_length;
575  if(sub_problem_blocksize>0) {
576  starti=current-current%sub_problem_blocksize;
577  endi=starti+sub_problem_blocksize;
578  }
579  for(unsigned int i=starti; i<endi; i++) {
582  if(is_allowed_by_weight_constraints(i)) {
583  // handle input_asym
584  for(unsigned k=0; k<this->input_asym.size(); ++k) {
585  for(unsigned int kk=0; kk<this->input_asym[k].size(); ++kk) {
586  if(i==this->input_asym[k][kk]) {
587  unsigned int k2=kk;
588  while(k2!=0) {
589  --k2;
591  // std::cout << "discarding " << std::endl;
593  break;
594  }
595  }
596  }
597  }
599  }
600  }
603  if(i+1>lowest_in_group) {
605  update_weights(i);
606  temparr[current]=i;
607  // for(unsigned bl=0; bl<block_length; ++bl)
608  // temparr[current*block_length+bl]=original[i*block_length+bl];
609  if(entry_accepted(current)) {
611  }
613  restore_weights(i);
614  }
615  }
616  }
617  }
618
619  template<class T>
621  : block_length(1), permutation_sign(1), sh_(*this), svh_(*this)
622  {
623  }
624
625  template<class T>
627  {
628  original.clear();
629  block_length=1;
630  permute_blocks.clear();
631  value_permute.clear();
632  permutation_sign=1;
633  sublengths.clear();
634  input_asym.clear();
635  sublengths_scattered.clear();
636  originals.clear();
637  multiplicity.clear();
638  }
639
640  template<class T>
642  {
643  std::cout << "collecting" << std::endl;
644  // Fill the hash map: entries which are equal have to sit in the same
645  // bin, but there may be other entries in that bin which still have to
646  // be separated.
647  std::multimap<long, unsigned int> hashmap;
648  for(unsigned int i=0; i<originals.size(); ++i)
649  hashmap.insert(std::pair<long, unsigned int>(hash(originals[i]), i));
650
651  // Collect equal vectors.
652  std::multimap<long, unsigned int>::iterator it=hashmap.begin(), thisbin1, thisbin2, tmpit;
653  while(it!=hashmap.end()) {
654  long current_hash=it->first;
655  thisbin1=it;
656  while(thisbin1!=hashmap.end() && thisbin1->first==current_hash) {
657  thisbin2=thisbin1;
658  ++thisbin2;
659  while(thisbin2!=hashmap.end() && thisbin2->first==current_hash) {
660  if(originals[(*thisbin1).second]==originals[(*thisbin2).second]) {
661  multiplicity[(*thisbin1).second]+=multiplicity[(*thisbin2).second];
662  multiplicity[(*thisbin2).second]=0;
663  tmpit=thisbin2;
664  ++tmpit;
665  hashmap.erase(thisbin2);
666  thisbin2=tmpit;
667  }
668  else ++thisbin2;
669  }
670  ++thisbin1;
671  }
672  it=thisbin1;
673  }
674
675  remove_multiplicity_zero();
676  }
677
678  template<class T>
680  {
681  std::vector<std::vector<T> > new_originals;
682  std::vector<int> new_multiplicity;
683  for(unsigned int k=0; k<originals.size(); ++k) {
684  if(multiplicity[k]!=0) {
685  new_originals.push_back(originals[k]);
686  new_multiplicity.push_back(multiplicity[k]);
687  }
688  }
689  originals=new_originals;
690  multiplicity=new_multiplicity;
691  }
692
693
694  template<class T>
696  {
697  unsigned int current_length=originals.size();
698  if(current_length==0) {
699  originals.push_back(original);
700  multiplicity.push_back(1);
701  current_length=1;
702  }
703
704  // Some options are mutually exclusive.
705  assert(permute_blocks.size()>0 || value_permute.size()>0);
706  assert(sublengths.size()==0 || sublengths_scattered.size()==0);
707
708  if(permute_blocks.size()==0) { // permute by value
709  assert(value_permute.size()!=0);
710
711  if(input_asym.size()==0 && sublengths_scattered.size()==0) {
712  // When permuting by value, we can do the permutation once,
713  // and then figure out (see vector_generated of symm_val_helper),
714  // for each permutation which is already stored in the symmetriser,
715  // how the objects are moved.
716
717  current_=current_length;
718  svh_.clear();
719  svh_.original=value_permute;
720  svh_.input_asym.clear();
721  svh_.sublengths=sublengths;
722  svh_.current_multiplicity=combin::vector_prod_fact(sublengths);
723  if(svh_.sublengths.size()==0)
724  svh_.set_unit_sublengths();
725
726  svh_.permute(start, end);
727  // Since we do not divide by the number of permutations, we need
728  // to adjust the multiplicity of all the originals.
729  // for(unsigned int i=0; i<current_; ++i)
730  // multiplicity[i] *= svh_.current_multiplicity;
731  }
732  else {
733  // However, when there is input_asym or sublength_scattered
734  // are present, we cannot just do the permutation on the
735  // values and then put them into all existing sets, since the
736  // overlap of input_asym with the objects to be permuted will
737  // be different for every set. Therefore, we have to apply
738  // the permutation algorithm separately to each and every set
739  // which is already stored in the symmetriser. We convert
740  // the problem to a permute-by-location problem.
741
742  for(unsigned int i=0; i<current_length; ++i) {
743  current_=i;
744  sh_.clear();
745  assert(sublengths.size()==0); // not yet implemented
746  std::vector<unsigned int> my_permute_blocks;
747
748  // Determine the location of the values.
749  for(unsigned int k=0; k<value_permute.size(); ++k) {
750  for(unsigned int m=0; m<originals[i].size(); ++m) {
751  if(originals[i][m]==value_permute[k]) {
752  my_permute_blocks.push_back(m); // FIXME: non-unit block length?
753  break;
754  }
755  }
756  }
757
758  // std::cout << "handling sublengths" << std::endl;
759  if(sublengths_scattered.size()>0) {
760  // Re-order my_permute_blocks in such a way that the objects which sit
761  // in one sublength_scattered range are consecutive. This does not make
762  // any difference for the sign.
763  sh_.sublengths.clear();
764  std::vector<unsigned int> reordered_permute_blocks;
765  for(unsigned int m=0; m<sublengths_scattered.size(); ++m) {
766  int overlap=0;
767  for(unsigned int mm=0; mm<sublengths_scattered[m].size(); ++mm) {
768  // std::cout << "trying to find " << sublengths_scattered[m][mm] << " " << std::flush;
769  std::vector<unsigned int>::iterator it=my_permute_blocks.begin();
770  while(it!=my_permute_blocks.end()) {
771  if((*it)==sublengths_scattered[m][mm]) {
772  // std::cout << " found " << std::endl;
773  reordered_permute_blocks.push_back(*it);
774  my_permute_blocks.erase(it);
775  ++overlap;
776  break;
777  }
778  ++it;
779  }
780  // std::cout << std::endl;
781  }
782  if(overlap>0)
783  sh_.sublengths.push_back(overlap);
784  }
785  std::vector<unsigned int>::iterator it=my_permute_blocks.begin();
786  while(it!=my_permute_blocks.end()) {
787  reordered_permute_blocks.push_back(*it);
788  // std::cout << "adding one" << std::endl;
789  sh_.sublengths.push_back(1);
790  ++it;
791  }
792  my_permute_blocks=reordered_permute_blocks;
793  // std::cout << "handled sublengths" << std::endl;
794  }
795
796  // Put to-be-permuted data in originals.
797  for(unsigned int k=0; k<my_permute_blocks.size(); ++k) {
798  for(unsigned int kk=0; kk<block_length; ++kk) {
799  sh_.original.push_back(originals[i][my_permute_blocks[k]+kk]);
800  }
801  }
802
803  combin::range_vector_t subprob_input_asym;
804  sh_.current_multiplicity=1;
805  if(input_asym.size()>0) {
806  // Make a proper input_asym which refers to object locations
807  // in the permute blocks array, rather than in the original
808  // array.
809  for(unsigned int k=0; k<input_asym.size(); ++k) {
810  range_t newrange;
811  for(unsigned int m=0; m<input_asym[k].size(); ++m) {
812  // search in my_permute_blocks
813  for(unsigned int kk=0; kk<my_permute_blocks.size(); ++kk)
814  if(my_permute_blocks[kk]==input_asym[k][m]) {
815  newrange.push_back(kk);
816  break;
817  }
818  }
819  if(newrange.size()>1) {
820  subprob_input_asym.push_back(newrange);
821  sh_.current_multiplicity*=fact(newrange.size());
822  }
823  }
824  }
825  if(sh_.sublengths.size()==0)
826  sh_.set_unit_sublengths();
827  sh_.current_multiplicity*=combin::vector_prod_fact(sh_.sublengths);
828
829  // debugging
830  // std::cout << "my_permute_blocks: ";
831  // for(unsigned int ii=0; ii<my_permute_blocks.size(); ++ii)
832  // std::cout << my_permute_blocks[ii] << " ";
833  // std::cout << std::endl;
834  // std::cout << "sublengths: ";
835  // for(unsigned int ii=0; ii<sh_.sublengths.size(); ++ii)
836  // std::cout << sh_.sublengths[ii] << " ";
837  // std::cout << std::endl;
838
839  // Debugging output:
840  // std::cout << sh_.current_multiplicity << " asym: ";
841  // if(subprob_input_asym.size()>0) {
842  // for(unsigned int k=0; k<subprob_input_asym[0].size(); ++k)
843  // std::cout << subprob_input_asym[0][k] << " ";
844  // std::cout << std::endl;
845  // std::cout << subprob_input_asym.size() << std::endl;
846  // }
847  // else std::cout << "no asym" << std::endl;
848
849  permute_blocks=my_permute_blocks;
850  sh_.block_length=block_length;
851  sh_.input_asym=subprob_input_asym;
852  sh_.permute(start, end);
853
854  // Since we do not divide by the number of permutations, we need
855  // to adjust the multiplicity of the original.
856  multiplicity[i]*=sh_.current_multiplicity;
857  permute_blocks.clear(); // restore just in case
858  // for(unsigned int m=0; m<originals.size(); ++m) {
859  // for(unsigned int mm=0; mm<originals[m].size(); ++mm)
860  // std::cout << originals[m][mm] << " ";
861  // std::cout << std::endl;
862  // }
863  // break;
864  }
865  }
866  }
867  else { // permute by location
868  assert(value_permute.size()==0);
869  assert(permute_blocks.size()>0);
870  // When permuting by location, we have to apply the permutation
871  // algorithm separately to each and every permutation which is
872  // already stored in the symmetriser.
873  for(unsigned int i=0; i<current_length; ++i) {
874  current_=i;
875  sh_.clear();
876  for(unsigned int k=0; k<permute_blocks.size(); ++k) {
877  for(unsigned int kk=0; kk<block_length; ++kk) {
878  sh_.original.push_back(originals[i][permute_blocks[k]+kk]);
879  }
880  // sh_.sublengths.push_back(1);
881  }
882  assert(sublengths.size()==0); // not yet implemented
883  // sh_.sublengths=sublengths;
884  if(sh_.sublengths.size()==0)
885  sh_.set_unit_sublengths();
886  sh_.block_length=block_length;
887  sh_.input_asym=input_asym;
888  sh_.permute(start, end);
889  }
890  }
891
892  if(start!=-1) { // if start is not the first, have to erase the first
893  originals.erase(originals.begin());
894  multiplicity.erase(multiplicity.begin());
895  }
896  }
897
898  template<class T>
899  const std::vector<T>& symmetriser<T>::operator[](unsigned int i) const
900  {
901  assert(i<originals.size());
902  return originals[i];
903  }
904
905  template<class T>
906  unsigned int symmetriser<T>::size() const
907  {
908  return originals.size();
909  }
910
911  template<class T>
912  range_t symmetriser<T>::values_to_locations(const std::vector<T>& values) const
913  {
914  range_t ret;
915  for(unsigned int i=0; i<values.size(); ++i) {
916  // std::cout << "finding " << values[i] << std::endl;
917  for(unsigned int j=0; j<value_permute.size(); ++j) {
918  // std::cout << value_permute[j] << " ";
919  if(value_permute[i]==value_permute[j]) {
920  // std::cout << "found" << std::endl;
921  ret.push_back(j);
922  break;
923  }
924  // std::cout << std::endl;
925  }
926  }
927  return ret;
928  }
929
930  template<class T>
932  : current_multiplicity(1), first_one(true), owner_(tt)
933  {
934  }
935
936  template<class T>
938  {
939  first_one=true;
941  }
942
943  template<class T>
944  void symm_val_helper<T>::vector_generated(const std::vector<unsigned int>& vec)
945  {
946  ++this->vector_generated_called_;
947  if(first_one) {
948  first_one=false;
949  }
950  else {
951  if((this->start_==-1 || this->vector_generated_called_ >= this->start_) &&
952  (this->end_==-1 || this->vector_generated_called_ < this->end_)) {
953
954  // Since we permuted by value, we can do this permutation in one
955  // shot on all previously generated sets.
956  for(unsigned int i=0; i<owner_.current_; ++i) {
957  // owner_.multiplicity[i] *= current_multiplicity;
958  owner_.originals.push_back(owner_.originals[i]);
959
960  // Take care of the multiplicity & sign.
961  int multiplicity=owner_.multiplicity[i] * current_multiplicity;
962  if(owner_.permutation_sign==-1)
963  multiplicity*=ordersign(vec.begin(), vec.end());
964  owner_.multiplicity.push_back(multiplicity); //sign==1?true:false);
965
966  // We now have to find the permuted objects in the larger
967  // "original" set, and re-order these appropriately.
968  unsigned int loc=owner_.originals.size()-1;
969  for(unsigned int j=0; j<vec.size(); ++j) {
970  for(unsigned int k=0; k<owner_.originals[i].size(); ++k) {
971  if(owner_.originals[i][k]==this->original[j]) {
972  owner_.originals[loc][k]=this->original[vec[j]];
973  break;
974  }
975  }
976  }
977  }
978  }
979  }
980  }
981
982  template<class T>
984  : current_multiplicity(1), first_one(true), owner_(tt)
985  {
986  }
987
988  template<class T>
990  {
991  first_one=true;
993  }
994
995  template<class T>
996  int symmetriser<T>::signature(unsigned int i) const
997  {
998  assert(i<multiplicity.size());
999  return multiplicity[i]; //?1:-1;
1000  }
1001
1002  template<class T>
1003  void symmetriser<T>::set_multiplicity(unsigned int i, int val)
1004  {
1005  assert(i<multiplicity.size());
1006  multiplicity[i]=val;
1007  }
1008
1009  template<class T>
1010  void symm_helper<T>::vector_generated(const std::vector<unsigned int>& vec)
1011  {
1012  ++this->vector_generated_called_;
1013  if(first_one) {
1014  first_one=false;
1015  }
1016  else {
1017  if((this->start_==-1 || this->vector_generated_called_ >= this->start_) &&
1018  (this->end_==-1 || this->vector_generated_called_ < this->end_)) {
1019
1020  // std::cout << "produced ";
1021  // for(unsigned int m=0; m<vec.size(); ++m)
1022  // std::cout << vec[m] << " ";
1023  // std::cout << std::endl;
1024
1025  // owner_.multiplicity[owner_.current_] *= current_multiplicity;
1026  owner_.originals.push_back(owner_.originals[owner_.current_]);
1027  unsigned int siz=owner_.originals.size()-1;
1028
1029  // Take care of the permutation sign.
1030  int multiplicity=owner_.multiplicity[owner_.current_] * current_multiplicity;
1031  if(owner_.permutation_sign==-1)
1032  multiplicity*=ordersign(vec.begin(), vec.end());
1033  owner_.multiplicity.push_back(multiplicity);
1034
1035  for(unsigned int k=0; k<owner_.permute_blocks.size(); ++k) {
1036  for(unsigned int kk=0; kk<owner_.block_length; ++kk) {
1037  assert(owner_.permute_blocks[k]+kk<owner_.originals[0].size());
1038  owner_.originals[siz][owner_.permute_blocks[k]+kk]=
1039  owner_.originals[owner_.current_][owner_.permute_blocks[vec[k]]+kk];
1040  }
1041  }
1042  }
1043  }
1044  }
1045
1046  template<class T>
1047  std::ostream& operator<<(std::ostream& str, const symmetriser<T>& sym)
1048  {
1049  for(unsigned int i=0; i<sym.size(); ++i) {
1050  for(unsigned int j=0; j<sym[i].size(); ++j) {
1051  str << sym[i][j] << " ";
1052  }
1053  str << " ";
1055  str << std::setw(2) << sym.signature(i) << std::endl;
1056  }
1057  return str;
1058  }
1059
1060  }
1061
1062
1063
Definition: Combinatorics.hh:57
unsigned int block_length
Definition: Combinatorics.hh:73
std::vector< unsigned int > sublengths
Definition: Combinatorics.hh:74
virtual void vector_generated(const std::vector< unsigned int > &)=0
bool multiple_pick
Definition: Combinatorics.hh:77
bool final_weight_constraints_check() const
Definition: Combinatorics.hh:522
void set_unit_sublengths()
Definition: Combinatorics.hh:450
long start_
Definition: Combinatorics.hh:88
virtual void clear()
Definition: Combinatorics.hh:387
std::vector< weights_t > weights
Definition: Combinatorics.hh:78
bool is_allowed_by_weight_constraints(unsigned int i)
Definition: Combinatorics.hh:510
std::vector< unsigned int > temparr
Definition: Combinatorics.hh:86
std::vector< weight_cond > weight_conditions
Definition: Combinatorics.hh:80
unsigned int multiplier(const std::vector< T > &) const
Definition: Combinatorics.hh:472
combinations_base()
Definition: Combinatorics.hh:297
std::vector< int > current_weight
Definition: Combinatorics.hh:89
unsigned int total_permutations() const
Definition: Combinatorics.hh:444
std::vector< int > max_weights
Definition: Combinatorics.hh:79
long vector_generated_called_
Definition: Combinatorics.hh:88
unsigned int sub_problem_blocksize
Definition: Combinatorics.hh:81
virtual void clear_results()
Definition: Combinatorics.hh:402
virtual ~combinations_base()
Definition: Combinatorics.hh:321
long end_
Definition: Combinatorics.hh:88
void nextstep(unsigned int current, unsigned int fromalgehad, unsigned int groupindex, std::vector< bool > algehad)
Definition: Combinatorics.hh:558
void update_weights(unsigned int i)
Definition: Combinatorics.hh:542
unsigned int sum_of_sublengths() const
Definition: Combinatorics.hh:435
range_vector_t input_asym
Definition: Combinatorics.hh:75
void restore_weights(unsigned int i)
Definition: Combinatorics.hh:550
virtual bool entry_accepted(unsigned int current) const
Definition: Combinatorics.hh:345
std::vector< T > original
Definition: Combinatorics.hh:76
void permute(long start=-1, long end=-1)
Definition: Combinatorics.hh:351
weight_cond
Definition: Combinatorics.hh:71
@ weight_equals
Definition: Combinatorics.hh:71
@ weight_greater
Definition: Combinatorics.hh:71
@ weight_less
Definition: Combinatorics.hh:71
Definition: Combinatorics.hh:101
permuted_sets_t::const_iterator const_iterator
Definition: Combinatorics.hh:104
int ordersign(unsigned int) const
Definition: Combinatorics.hh:458
virtual ~combinations()
Definition: Combinatorics.hh:326
unsigned int size() const
Definition: Combinatorics.hh:429
const std::vector< T > & operator[](unsigned int) const
Definition: Combinatorics.hh:422
virtual void clear()
Definition: Combinatorics.hh:408
combinations()
Definition: Combinatorics.hh:309
virtual void vector_generated(const std::vector< unsigned int > &)
Definition: Combinatorics.hh:331
std::vector< std::vector< T > > permuted_sets_t
Definition: Combinatorics.hh:103
virtual void clear_results()
Definition: Combinatorics.hh:415
permuted_sets_t storage
Definition: Combinatorics.hh:121
combinations(const std::vector< T > &)
Definition: Combinatorics.hh:315
unsigned int multiplier(unsigned int) const
Definition: Combinatorics.hh:466
Definition: Combinatorics.hh:128
symm_helper(symmetriser< T > &)
Definition: Combinatorics.hh:983
symmetriser< T > & owner_
Definition: Combinatorics.hh:136
bool first_one
Definition: Combinatorics.hh:135
int current_multiplicity
Definition: Combinatorics.hh:133
virtual void vector_generated(const std::vector< unsigned int > &)
Definition: Combinatorics.hh:1010
virtual void clear()
Definition: Combinatorics.hh:989
Definition: Combinatorics.hh:141
int current_multiplicity
Definition: Combinatorics.hh:146
symm_val_helper(symmetriser< T > &)
Definition: Combinatorics.hh:931
bool first_one
Definition: Combinatorics.hh:148
virtual void clear()
Definition: Combinatorics.hh:937
virtual void vector_generated(const std::vector< unsigned int > &)
Definition: Combinatorics.hh:944
symmetriser< T > & owner_
Definition: Combinatorics.hh:149
Definition: Combinatorics.hh:154
unsigned int size() const
Definition: Combinatorics.hh:906
void clear()
Definition: Combinatorics.hh:626
std::vector< std::vector< T > > originals
Definition: Combinatorics.hh:188
void apply_symmetry(long start=-1, long end=-1)
Definition: Combinatorics.hh:695
symm_helper< T > sh_
Definition: Combinatorics.hh:185
unsigned int current_
Definition: Combinatorics.hh:187
range_t values_to_locations(const std::vector< T > &values) const
Convert vectors of values to vectors of locations in the original (mainly useful to create input_asym...
Definition: Combinatorics.hh:912
std::vector< unsigned int > sublengths
Definition: Combinatorics.hh:165
int permutation_sign
Definition: Combinatorics.hh:163
range_vector_t sublengths_scattered
Definition: Combinatorics.hh:167
std::vector< unsigned int > permute_blocks
Definition: Combinatorics.hh:161
std::vector< T > original
Definition: Combinatorics.hh:159
symmetriser()
Definition: Combinatorics.hh:620
range_vector_t input_asym
Definition: Combinatorics.hh:166
std::vector< int > multiplicity
Definition: Combinatorics.hh:189
void collect()
Collect equal entries, and adjust the multiplier field accordingly.
Definition: Combinatorics.hh:641
std::vector< T > value_permute
Definition: Combinatorics.hh:162
void remove_multiplicity_zero()
Definition: Combinatorics.hh:679
const std::vector< T > & operator[](unsigned int) const
Definition: Combinatorics.hh:899
void set_multiplicity(unsigned int pos, int val)
Definition: Combinatorics.hh:1003
unsigned int block_length
Definition: Combinatorics.hh:160
int signature(unsigned int) const
Definition: Combinatorics.hh:996
symm_val_helper< T > svh_
Definition: Combinatorics.hh:186
Definition: Combinatorics.hh:37
std::vector< range_t > range_vector_t
Definition: Combinatorics.hh:40
int determine_intersection_ranges(const range_vector_t &prod, const range_vector_t &indv, range_vector_t &target)
Definition: Combinatorics.cc:32
unsigned long factorial(unsigned int x)
Definition: Combinatorics.cc:23
std::vector< unsigned int > range_t
Definition: Combinatorics.hh:39
unsigned long vector_prod_fact(const std::vector< unsigned int > &)
product of factorials of elements
Definition: Combinatorics.cc:72
std::vector< int > weights_t
Definition: Combinatorics.hh:41
unsigned long vector_prod(const std::vector< unsigned int > &)
product of elements
Definition: Combinatorics.cc:64
std::ostream & operator<<(std::ostream &str, const symmetriser< T > &sym)
Definition: Combinatorics.hh:1047
T fact(T x)
Definition: Combinatorics.hh:283
long vector_sum(const std::vector< int > &)
sum of elements
Definition: Combinatorics.cc:56
bool operator==(const std::vector< unsigned int > &, const std::vector< unsigned int > &)
Definition: Combinatorics.cc:80
long hash(const std::vector< unsigned int > &)
compute a hash value for a vector of unsigned ints
Definition: Combinatorics.cc:90
int ordersign(iterator1 b1, iterator1 e1, iterator2 b2, iterator2 e2, int stepsize=1)
Definition: Combinatorics.hh:224
end
Definition: nevaluate.py:22
start
Definition: nevaluate.py:20
int k
Definition: passing.cc:4