Combinatorics.hh

Go to the documentation of this file.

/* 

    Cadabra: a field-theory motivated computer algebra system.
    Copyright (C) 2001-2014  Kasper Peeters <kasper.peeters@phi-sci.com>

   This program is free software: you can redistribute it and/or
   modify it under the terms of the GNU General Public License as
   published by the Free Software Foundation, either version 3 of the
   License, or (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
*/

/*
                                  length vector
   normal combinations:     one element, value=total length.
    normal permutations:     n elements, each equal to 1.

*/

#pragma once

#include <vector>
#include <cassert>
#include <algorithm>
#include <iomanip>
#include <iostream>
#include <map>

namespace combin {

typedef std::vector<unsigned int>  range_t;
typedef std::vector<range_t>       range_vector_t;
typedef std::vector<int>           weights_t;
    
unsigned long factorial(unsigned int x);
long          vector_sum(const std::vector<int>&);
unsigned long vector_prod(const std::vector<unsigned int>&);
unsigned long vector_prod_fact(const std::vector<unsigned int>&);

bool operator==(const std::vector<unsigned int>&, const std::vector<unsigned int>&);

long hash(const std::vector<unsigned int>&);

template<class T> 
class combinations_base {
    public:
        combinations_base();
        combinations_base(const std::vector<T>&);
        virtual ~combinations_base();

        void         permute(long start=-1, long end=-1);
        virtual void clear();
        virtual void clear_results();
        unsigned int sum_of_sublengths() const;
        void         set_unit_sublengths();
        unsigned int multiplier(const std::vector<T>&) const;
        unsigned int total_permutations() const; // including the ones not stored

        enum weight_cond { weight_equals, weight_less, weight_greater };

        unsigned int               block_length;
        std::vector<unsigned int>  sublengths;
        range_vector_t             input_asym;
      std::vector<T>             original;
        bool                       multiple_pick;
        std::vector<weights_t>     weights;
        std::vector<int>           max_weights;
        std::vector<weight_cond>   weight_conditions;
        unsigned int               sub_problem_blocksize; // when non-zero, do permutations within

    protected:
        virtual void   vector_generated(const std::vector<unsigned int>&)=0;
        virtual bool   entry_accepted(unsigned int current) const;
        std::vector<unsigned int>  temparr;

        long                       start_, end_, vector_generated_called_;
        std::vector<int>           current_weight;
    private:
        bool is_allowed_by_weight_constraints(unsigned int i);
        bool final_weight_constraints_check() const;
        void update_weights(unsigned int i);
        void restore_weights(unsigned int i);
      void nextstep(unsigned int current, unsigned int fromalgehad, unsigned int groupindex, 
                          std::vector<bool> algehad);

};

template<class T>
class combinations : public combinations_base<T> {
   public:
        typedef typename std::vector<std::vector<T> >    permuted_sets_t;
      typedef typename permuted_sets_t::const_iterator const_iterator;

      combinations();
        combinations(const std::vector<T>&);
        virtual ~combinations();

        virtual void           clear();
        virtual void           clear_results();
        const std::vector<T>&  operator[](unsigned int) const;
        int                    ordersign(unsigned int) const;
        unsigned int           size() const;
        unsigned int           multiplier(unsigned int) const;
        
    protected:
        virtual void vector_generated(const std::vector<unsigned int>&);

   private:
        permuted_sets_t storage;
};

template<class T>
class symmetriser;

template<class T>
class symm_helper : public combinations_base<T> {
    public: 
        symm_helper(symmetriser<T>&);
        virtual void clear();

        int             current_multiplicity;
    protected:
        bool            first_one;
        symmetriser<T>& owner_;
        virtual void vector_generated(const std::vector<unsigned int>&);
};

template<class T>
class symm_val_helper : public combinations_base<T> {
    public: 
        symm_val_helper(symmetriser<T>&);
        virtual void clear();

        int             current_multiplicity;
    protected:
        bool            first_one;
        symmetriser<T>& owner_;
        virtual void vector_generated(const std::vector<unsigned int>&);
};

template<class T>
class symmetriser {
    public:
        symmetriser();
        void apply_symmetry(long start=-1, long end=-1);

        std::vector<T>            original; 
        unsigned int              block_length;
        std::vector<unsigned int> permute_blocks;   // offset in unit elements! (not in blocks)
        std::vector<T>            value_permute;
        int                       permutation_sign;

        std::vector<unsigned int> sublengths; // refers to position within permute_blocks
        range_vector_t            input_asym; // as in combinations_base
        range_vector_t            sublengths_scattered; // sublengths, in original, but not connected.

        range_t                   values_to_locations(const std::vector<T>& values) const;

        const std::vector<T>&     operator[](unsigned int) const;
        int                       signature(unsigned int) const;
        void                      set_multiplicity(unsigned int pos, int val);
        unsigned int              size() const;
        void                      clear();
        void                      collect();
        void                      remove_multiplicity_zero();

        friend class symm_helper<T>;
        friend class symm_val_helper<T>;
    private:
        symm_helper<T>               sh_;
        symm_val_helper<T>           svh_;
        unsigned int                 current_;
        std::vector<std::vector<T> > originals;
        std::vector<int>             multiplicity;
};

int determine_intersection_ranges(const range_vector_t& prod,
                                             const range_vector_t& indv,
                                             range_vector_t& target);

template<class iterator1, class iterator2>
int ordersign(iterator1 b1, iterator1 e1, iterator2 b2, iterator2 e2, int stepsize=1);

template<class iterator1>
int ordersign(iterator1 b1, iterator1 e1);

template<class T>
T fact(T x);

template<class T>
std::ostream& operator<<(std::ostream& str, const symmetriser<T>& sym);


/* 
I assume PI consists of the integers 1 to N.
It can be done with O(N) comparisons and transpositions of integers 
in the list.

sign:= 1;
for i from 1 to N do
  while PI[i] <> i do
    interchange PI[i] and PI[PI[i]];
    sign:= -sign
  od 
od
*/

template<class iterator1, class iterator2>
int ordersign(iterator1 b1, iterator1 e1, iterator2 b2, iterator2 e2, int stepsize) 
    {
    int sign=1;
    std::vector<bool> crossedoff(std::distance(b1,e1),false);
    while(b1!=e1) {
        int otherpos=0;
        iterator2 it=b2;
        while(it!=e2) {
            if( (*it)==(*b1) && crossedoff[otherpos]==false) {
                crossedoff[otherpos]=true;
                break;
                }
            else {
                if(!crossedoff[otherpos])
                    sign=-sign;
                }
            it+=stepsize;
            ++otherpos;
            }
        b1+=stepsize;
        }
    return sign;
    }

//template<class iterator1, class iterator2, class comparator>
//int ordersign(iterator1 b1, iterator1 e1, iterator2 b2, iterator2 e2, comparator cmp, int stepsize) 
//  {
//  int sign=1;
//  std::vector<bool> crossedoff(std::distance(b1,e1),false);
//  while(b1!=e1) {
//      int otherpos=0;
//      iterator2 it=b2;
//      while(it!=e2) {
//          if(cmp((*it), (*b1)) && crossedoff[otherpos]==false) {
//              crossedoff[otherpos]=true;
//              break;
//              }
//          else {
//              if(!crossedoff[otherpos])
//                  sign=-sign;
//              }
//          it+=stepsize;
//          ++otherpos;
//          }
//      b1+=stepsize;
//      }
//  return sign;
//  }

template<class iterator1>
int ordersign(iterator1 b1, iterator1 e1)
    {
    std::vector<unsigned int> fil;
    for(int k=0; k<distance(b1,e1); ++k)
        fil.push_back(k);
    return ordersign(fil.begin(), fil.end(), b1, e1);
    }

template<class T>
T fact(T x) 
    {
    T ret=1;
    assert(x>=0);
    while(x!=0) { ret*=x--; }
    return ret;
    }


// Implementations

template<class T>
combinations_base<T>::combinations_base()
    : block_length(1), multiple_pick(false), sub_problem_blocksize(0)
   { 
   }

template<class T>
combinations_base<T>::combinations_base(const std::vector<T>& oa)
    : block_length(1), original(oa), multiple_pick(false), sub_problem_blocksize(0)
   { 
   }

template<class T>
combinations<T>::combinations()
    : combinations_base<T>()
    {
    }

template<class T>
combinations<T>::combinations(const std::vector<T>& oa)
    : combinations_base<T>(oa)
    {
    }

template<class T>
combinations_base<T>::~combinations_base()
    {
    }

template<class T>
combinations<T>::~combinations()
    {
    }

template<class T>
void combinations<T>::vector_generated(const std::vector<unsigned int>& toadd) 
    {
    ++this->vector_generated_called_;
    if((this->start_==-1 || this->vector_generated_called_ >= this->start_) && 
        (this->end_==-1   || this->vector_generated_called_ < this->end_)) {
        std::vector<T> newone(toadd.size()*this->block_length);
        for(unsigned int i=0; i<toadd.size(); ++i) 
            for(unsigned int bl=0; bl<this->block_length; ++bl) 
                newone[i*this->block_length+bl]=this->original[toadd[i]*this->block_length+bl];
        storage.push_back(newone);
        }
    }

template<class T>
bool combinations_base<T>::entry_accepted(unsigned int) const
    {
    return true;
    }

template<class T>
void combinations_base<T>::permute(long start, long end)
    {
    start_=start;
    end_=end;
    vector_generated_called_=-1;

    // Initialise weight handling.
    current_weight.clear();
    current_weight.resize(weights.size(), 0);
    for(unsigned int i=0; i<weights.size(); ++i) 
        assert(weights[i].size() == original.size()/block_length);
    if(weights.size()>0) {
        if(weight_conditions.size()==0)
            weight_conditions.resize(weights.size(), weight_equals);
        else assert(weight_conditions.size()==weights.size());
        }
    else assert(weight_conditions.size()==0);

    // Sublength handling.
    assert(sublengths.size()!=0);
    unsigned int len=sum_of_sublengths();

    // Consistency checks.
    assert(original.size()%block_length==0);
    if(!multiple_pick)
        assert(len*block_length<=original.size());

    for(unsigned int i=0; i<this->input_asym.size(); ++i)
        std::sort(this->input_asym[i].begin(), this->input_asym[i].end());

   temparr=std::vector<unsigned int>(len/* *block_length*/);
   std::vector<bool> algehad(original.size()/block_length,false);
   nextstep(0,0,0,algehad);
    }

template<class T>
void combinations_base<T>::clear()
    {
    block_length=1;
    sublengths.clear();
    this->input_asym.clear();
    original.clear();
    weights.clear();
    max_weights.clear();
    weight_conditions.clear();
    sub_problem_blocksize=0;
    temparr.clear();
    current_weight.clear();
    }

template<class T>
void combinations_base<T>::clear_results()
    {
    temparr.clear();
    }

template<class T>
void combinations<T>::clear()
    {
   storage.clear();
    combinations_base<T>::clear();
    }

template<class T>
void combinations<T>::clear_results()
    {
    storage.clear();
    combinations_base<T>::clear_results();
    }

template<class T>
const std::vector<T>& combinations<T>::operator[](unsigned int i) const
    {
    assert(i<storage.size());
    return storage[i];
    }

template<class T>
unsigned int combinations<T>::size() const
    {
    return storage.size();
    }

template<class T>
unsigned int combinations_base<T>::sum_of_sublengths() const
    {
    unsigned int ret=0;
    for(unsigned int i=0; i<sublengths.size(); ++i)
        ret+=sublengths[i];
    return ret;
    }

template<class T>
unsigned int combinations_base<T>::total_permutations() const
    {
    return vector_generated_called_+1;
    }

template<class T>
void combinations_base<T>::set_unit_sublengths() 
    {
    sublengths.clear();
    for(unsigned int i=0; i<original.size()/block_length; ++i)
        sublengths.push_back(1);
    }

template<class T>
int combinations<T>::ordersign(unsigned int num) const
    {
    assert(num<storage.size());
    return combin::ordersign(storage[0].begin(), storage[0].end(),
                             storage[num].begin(), storage[num].end(), this->block_length);
    }

template<class T>
unsigned int combinations<T>::multiplier(unsigned int num) const
    {
    return combinations_base<T>::multiplier(this->storage[num]);
    }

template<class T>
unsigned int combinations_base<T>::multiplier(const std::vector<T>& stor) const
    {
    unsigned long numerator=1;
    for(unsigned int i=0; i<this->input_asym.size(); ++i) 
        numerator*=fact(this->input_asym[i].size());

    unsigned long denominator=1;
    for(unsigned int i=0; i<this->input_asym.size(); ++i) {
        // for each input asym, and for each output asym, count
        // the number of overlap elements.
        unsigned int current=0;
        for(unsigned int k=0; k<this->sublengths.size(); ++k) {
            if(this->sublengths[k]>1) {
                unsigned int overlap=0;
                for(unsigned int slc=0; slc<this->sublengths[k]; ++slc) {
                    for(unsigned int j=0; j<this->input_asym[i].size(); ++j) {
                        unsigned int index=0;
                        while(!(stor[current]==this->original[index]))
                            ++index;
                        if(index==this->input_asym[i][j])
                            ++overlap;
                            }
                    ++current;
                    }
                if(overlap>0)
                    denominator*=fact(overlap);
                // FIXME: for each overlap thus found, divide out by a factor
                // due to the fact that output asym ranges can overlap.
// well, that's not right either.
                }
            else ++current;
            }
        }
    
    return numerator/denominator;
    }

template<class T>
bool combinations_base<T>::is_allowed_by_weight_constraints(unsigned int i)
    {
    if(weights.size()==0) return true;
    for(unsigned int cn=0; cn<current_weight.size(); ++cn) {
        if(weight_conditions[cn]==weight_less)
            if(current_weight[cn]+weights[cn][i] >= max_weights[cn])
                return false;
        }
    return true;
    }

template<class T>
bool combinations_base<T>::final_weight_constraints_check() const
    {
    for(unsigned int cn=0; cn<current_weight.size(); ++cn) {
        switch(weight_conditions[cn]) {
            case weight_equals:
                if(current_weight[cn]!=max_weights[cn])
                    return false;
                break;
            case weight_less:
                break;
            case weight_greater:
                if(current_weight[cn]<=max_weights[cn])
                    return false;
                break;
            }
        }
    return true;
    }

template<class T>
void combinations_base<T>::update_weights(unsigned int i)
    {
    if(weights.size()==0) return;
    for(unsigned int cn=0; cn<current_weight.size(); ++cn)
        current_weight[cn]+=weights[cn][i];
    }

template<class T>
void combinations_base<T>::restore_weights(unsigned int i)
    {
    if(weights.size()==0) return;
    for(unsigned int cn=0; cn<current_weight.size(); ++cn)
        current_weight[cn]-=weights[cn][i];
    }

template<class T>
void combinations_base<T>::nextstep(unsigned int current, unsigned int lowest_in_group, unsigned int groupindex, 
                                         std::vector<bool> algehad)
   {
    unsigned int grouplen=0;
    for(unsigned int i=0; i<=groupindex; ++i)
        grouplen+=sublengths[i];
   if(current==grouplen) { // group is filled
        ++groupindex;
        if(groupindex==sublengths.size()) {
            if(final_weight_constraints_check())
                vector_generated(temparr);
            return;
            }
        lowest_in_group=0;
      }
    
    unsigned int starti=0, endi=original.size()/block_length;
    if(sub_problem_blocksize>0) {
        starti=current-current%sub_problem_blocksize;
        endi=starti+sub_problem_blocksize;
        }
   for(unsigned int i=starti; i<endi; i++) {
        if(!algehad[i] || multiple_pick) {
            bool discard=false;
            if(is_allowed_by_weight_constraints(i)) {
                // handle input_asym
                for(unsigned k=0; k<this->input_asym.size(); ++k) {
                    for(unsigned int kk=0; kk<this->input_asym[k].size(); ++kk) {
                        if(i==this->input_asym[k][kk]) {
                            unsigned int k2=kk;
                            while(k2!=0) {
                                --k2;
                                if(!algehad[this->input_asym[k][k2]]) {
//                                  std::cout << "discarding " << std::endl;
                                    discard=true;
                                    break;
                                    }
                                }
                            }
                        }
                    if(discard) break;
                    }
                }
            else discard=true;
            if(!discard)
                if(i+1>lowest_in_group) {
                    algehad[i]=true;
                    update_weights(i);
                    temparr[current]=i;
//                  for(unsigned bl=0; bl<block_length; ++bl) 
//                      temparr[current*block_length+bl]=original[i*block_length+bl];
                    if(entry_accepted(current)) {
                        nextstep(current+1, i, groupindex, algehad);
                        }
                    algehad[i]=false;
                    restore_weights(i);
                    }
            }
      }
   }

template<class T>
symmetriser<T>::symmetriser()
    : block_length(1), permutation_sign(1), sh_(*this), svh_(*this)
    {
    }

template<class T>
void symmetriser<T>::clear()
    {
    original.clear();
    block_length=1;
    permute_blocks.clear();
    value_permute.clear();
    permutation_sign=1;
    sublengths.clear();
    input_asym.clear();
    sublengths_scattered.clear();
    originals.clear();
    multiplicity.clear();
    }

template<class T>
void symmetriser<T>::collect()
    {
    std::cout << "collecting" << std::endl;
    // Fill the hash map: entries which are equal have to sit in the same
    // bin, but there may be other entries in that bin which still have to
    // be separated.
    std::multimap<long, unsigned int> hashmap;
    for(unsigned int i=0; i<originals.size(); ++i)
        hashmap.insert(std::pair<long, unsigned int>(hash(originals[i]), i));

    // Collect equal vectors.
    std::multimap<long, unsigned int>::iterator it=hashmap.begin(), thisbin1, thisbin2, tmpit;
    while(it!=hashmap.end()) {
        long current_hash=it->first;
        thisbin1=it;
        while(thisbin1!=hashmap.end() && thisbin1->first==current_hash) {
            thisbin2=thisbin1;
            ++thisbin2;
            while(thisbin2!=hashmap.end() && thisbin2->first==current_hash) {
                if(originals[(*thisbin1).second]==originals[(*thisbin2).second]) {
                    multiplicity[(*thisbin1).second]+=multiplicity[(*thisbin2).second];
                    multiplicity[(*thisbin2).second]=0;
                    tmpit=thisbin2;
                    ++tmpit;
                    hashmap.erase(thisbin2);
                    thisbin2=tmpit;
                    }
                else ++thisbin2;
                }
            ++thisbin1;
            }
        it=thisbin1;
        }

    remove_multiplicity_zero();
    }

template<class T>
void symmetriser<T>::remove_multiplicity_zero()
    {
    std::vector<std::vector<T> > new_originals;
    std::vector<int>             new_multiplicity;
    for(unsigned int k=0; k<originals.size(); ++k) {
        if(multiplicity[k]!=0) {
            new_originals.push_back(originals[k]);
            new_multiplicity.push_back(multiplicity[k]);
            }
        }
    originals=new_originals;
    multiplicity=new_multiplicity;
    }


template<class T>
void symmetriser<T>::apply_symmetry(long start, long end)
    {
    unsigned int current_length=originals.size();
    if(current_length==0) {
        originals.push_back(original);
        multiplicity.push_back(1);
        current_length=1;
        }

    // Some options are mutually exclusive.
    assert(permute_blocks.size()>0 || value_permute.size()>0);
    assert(sublengths.size()==0 || sublengths_scattered.size()==0);

    if(permute_blocks.size()==0) { // permute by value
        assert(value_permute.size()!=0);

        if(input_asym.size()==0 && sublengths_scattered.size()==0) {
            // When permuting by value, we can do the permutation once,
            // and then figure out (see vector_generated of symm_val_helper),
            // for each permutation which is already stored in the symmetriser,
            // how the objects are moved.
            
            current_=current_length;
            svh_.clear();
            svh_.original=value_permute;
            svh_.input_asym.clear();
            svh_.sublengths=sublengths;
            svh_.current_multiplicity=combin::vector_prod_fact(sublengths);
            if(svh_.sublengths.size()==0)
                svh_.set_unit_sublengths();

            svh_.permute(start, end);
            // Since we do not divide by the number of permutations, we need
            // to adjust the multiplicity of all the originals.
            //      for(unsigned int i=0; i<current_; ++i)
//              multiplicity[i] *= svh_.current_multiplicity;
            }
        else {
            // However, when there is input_asym or sublength_scattered
            // are present, we cannot just do the permutation on the
            // values and then put them into all existing sets, since the
            // overlap of input_asym with the objects to be permuted will
            // be different for every set.  Therefore, we have to apply
            // the permutation algorithm separately to each and every set
            // which is already stored in the symmetriser.  We convert
            // the problem to a permute-by-location problem.

            for(unsigned int i=0; i<current_length; ++i) {
                current_=i;
                sh_.clear();
                assert(sublengths.size()==0); // not yet implemented
                std::vector<unsigned int> my_permute_blocks;

                // Determine the location of the values.
                for(unsigned int k=0; k<value_permute.size(); ++k) {
                    for(unsigned int m=0; m<originals[i].size(); ++m) {
                        if(originals[i][m]==value_permute[k]) {
                            my_permute_blocks.push_back(m); // FIXME: non-unit block length?
                            break;
                            }
                        }
                    }

//              std::cout << "handling sublengths" << std::endl;
                if(sublengths_scattered.size()>0) {
                    // Re-order my_permute_blocks in such a way that the objects which sit
                    // in one sublength_scattered range are consecutive. This does not make
                    // any difference for the sign.
                    sh_.sublengths.clear();
                    std::vector<unsigned int> reordered_permute_blocks;
                    for(unsigned int m=0; m<sublengths_scattered.size(); ++m) {
                        int overlap=0;
                        for(unsigned int mm=0; mm<sublengths_scattered[m].size(); ++mm) {
//                          std::cout << "trying to find " << sublengths_scattered[m][mm] << " " << std::flush;
                            std::vector<unsigned int>::iterator it=my_permute_blocks.begin();
                            while(it!=my_permute_blocks.end()) {
                                if((*it)==sublengths_scattered[m][mm]) {
//                                  std::cout << " found " << std::endl;
                                    reordered_permute_blocks.push_back(*it);
                                    my_permute_blocks.erase(it);
                                    ++overlap;
                                    break;
                                    }
                                ++it;
                                }
//                          std::cout << std::endl;
                            }
                        if(overlap>0)
                            sh_.sublengths.push_back(overlap);
                        }
                    std::vector<unsigned int>::iterator it=my_permute_blocks.begin();
                    while(it!=my_permute_blocks.end()) {
                        reordered_permute_blocks.push_back(*it);
//                      std::cout << "adding one" << std::endl;
                        sh_.sublengths.push_back(1);
                        ++it;
                        }
                    my_permute_blocks=reordered_permute_blocks;
//                  std::cout << "handled sublengths" << std::endl;
                    }
                
                // Put to-be-permuted data in originals.
                for(unsigned int k=0; k<my_permute_blocks.size(); ++k) {
                    for(unsigned int kk=0; kk<block_length; ++kk) {
                        sh_.original.push_back(originals[i][my_permute_blocks[k]+kk]);
                        }
                    }

                combin::range_vector_t subprob_input_asym;
                sh_.current_multiplicity=1;
                if(input_asym.size()>0) {
                    // Make a proper input_asym which refers to object locations
                    // in the permute blocks array, rather than in the original 
                    // array.
                    for(unsigned int k=0; k<input_asym.size(); ++k) {
                        range_t newrange;
                        for(unsigned int m=0; m<input_asym[k].size(); ++m) {
                            // search in my_permute_blocks
                            for(unsigned int kk=0; kk<my_permute_blocks.size(); ++kk)
                                if(my_permute_blocks[kk]==input_asym[k][m]) {
                                    newrange.push_back(kk);
                                    break;
                                    }
                            }
                        if(newrange.size()>1) {
                            subprob_input_asym.push_back(newrange);
                            sh_.current_multiplicity*=fact(newrange.size());
                            }
                        }
                    }
                if(sh_.sublengths.size()==0)
                    sh_.set_unit_sublengths();
                sh_.current_multiplicity*=combin::vector_prod_fact(sh_.sublengths);

                // debugging
//              std::cout << "my_permute_blocks: ";
//              for(unsigned int ii=0; ii<my_permute_blocks.size(); ++ii)
//                  std::cout << my_permute_blocks[ii] << " ";
//              std::cout << std::endl;
//              std::cout << "sublengths: ";
//              for(unsigned int ii=0; ii<sh_.sublengths.size(); ++ii)
//                  std::cout << sh_.sublengths[ii] << " ";
//              std::cout << std::endl;

                // Debugging output:
//              std::cout << sh_.current_multiplicity << " asym: ";
//              if(subprob_input_asym.size()>0) {
//                  for(unsigned int k=0; k<subprob_input_asym[0].size(); ++k)
//                      std::cout << subprob_input_asym[0][k] << " ";
//                  std::cout << std::endl;
//                  std::cout << subprob_input_asym.size() << std::endl;
//                  }
//              else std::cout << "no asym" << std::endl;
                
                permute_blocks=my_permute_blocks;
                sh_.block_length=block_length;
                sh_.input_asym=subprob_input_asym;
                sh_.permute(start, end);
                
                // Since we do not divide by the number of permutations, we need
                // to adjust the multiplicity of the original.
                multiplicity[i]*=sh_.current_multiplicity;
                permute_blocks.clear(); // restore just in case
//              for(unsigned int m=0; m<originals.size(); ++m) {
//                  for(unsigned int mm=0; mm<originals[m].size(); ++mm)
//                      std::cout << originals[m][mm] << " ";
//                  std::cout << std::endl;
//                  }
//              break;
                }
            }
        }
    else {                         // permute by location
        assert(value_permute.size()==0);
        assert(permute_blocks.size()>0);
        // When permuting by location, we have to apply the permutation
        // algorithm separately to each and every permutation which is
        // already stored in the symmetriser.
        for(unsigned int i=0; i<current_length; ++i) {
            current_=i;
            sh_.clear();
            for(unsigned int k=0; k<permute_blocks.size(); ++k) {
                for(unsigned int kk=0; kk<block_length; ++kk) {
                    sh_.original.push_back(originals[i][permute_blocks[k]+kk]);
                    }
//              sh_.sublengths.push_back(1);
                }
            assert(sublengths.size()==0); // not yet implemented
            // sh_.sublengths=sublengths;
            if(sh_.sublengths.size()==0)
                sh_.set_unit_sublengths();
            sh_.block_length=block_length;
            sh_.input_asym=input_asym;
            sh_.permute(start, end);
            }
        }
    
    if(start!=-1) { // if start is not the first, have to erase the first
        originals.erase(originals.begin());
        multiplicity.erase(multiplicity.begin());
        }
    }

template<class T>
const std::vector<T>& symmetriser<T>::operator[](unsigned int i) const
    {
    assert(i<originals.size());
    return originals[i];
    }

template<class T>
unsigned int symmetriser<T>::size() const
    {
    return originals.size();
    }

template<class T>
range_t symmetriser<T>::values_to_locations(const std::vector<T>& values) const
    {
    range_t ret;
    for(unsigned int i=0; i<values.size(); ++i) {
//      std::cout << "finding " << values[i] << std::endl;
        for(unsigned int j=0; j<value_permute.size(); ++j) {
//          std::cout << value_permute[j] << " ";
            if(value_permute[i]==value_permute[j]) {
//              std::cout << "found" << std::endl;
                ret.push_back(j);
                break;
                }
//          std::cout << std::endl;
            }
        }
    return ret;
    }

template<class T>
symm_val_helper<T>::symm_val_helper(symmetriser<T>& tt)
    : current_multiplicity(1), first_one(true), owner_(tt)
    {
    }

template<class T>
void symm_val_helper<T>::clear() 
    {
    first_one=true;
    combinations_base<T>::clear();
    }

template<class T>
void symm_val_helper<T>::vector_generated(const std::vector<unsigned int>& vec)
    {
    ++this->vector_generated_called_;
    if(first_one) {
        first_one=false;
        }
    else {
        if((this->start_==-1 || this->vector_generated_called_ >= this->start_) && 
            (this->end_==-1   || this->vector_generated_called_ < this->end_)) {
            
            // Since we permuted by value, we can do this permutation in one
            // shot on all previously generated sets.
            for(unsigned int i=0; i<owner_.current_; ++i) {
//              owner_.multiplicity[i] *= current_multiplicity;
                owner_.originals.push_back(owner_.originals[i]);
                
                // Take care of the multiplicity & sign.
                int multiplicity=owner_.multiplicity[i] * current_multiplicity;
                if(owner_.permutation_sign==-1)
                    multiplicity*=ordersign(vec.begin(), vec.end());
                owner_.multiplicity.push_back(multiplicity); //sign==1?true:false);
                
                // We now have to find the permuted objects in the larger
                // "original" set, and re-order these appropriately.
                unsigned int loc=owner_.originals.size()-1;
                for(unsigned int j=0; j<vec.size(); ++j) {
                    for(unsigned int k=0; k<owner_.originals[i].size(); ++k) {
                        if(owner_.originals[i][k]==this->original[j]) {
                            owner_.originals[loc][k]=this->original[vec[j]];
                            break;
                            }
                        }
                    }
                }
            }
        }
    }
                    
template<class T>
symm_helper<T>::symm_helper(symmetriser<T>& tt)
    : current_multiplicity(1), first_one(true), owner_(tt)
    {
    }

template<class T>
void symm_helper<T>::clear()
    {
    first_one=true;
    combinations_base<T>::clear();
    }

template<class T>
int symmetriser<T>::signature(unsigned int i) const
    {
    assert(i<multiplicity.size());
    return multiplicity[i]; //?1:-1;
    }

template<class T>
void symmetriser<T>::set_multiplicity(unsigned int i, int val) 
    {
    assert(i<multiplicity.size());
    multiplicity[i]=val;
    }

template<class T>
void symm_helper<T>::vector_generated(const std::vector<unsigned int>& vec)
    {
    ++this->vector_generated_called_;
    if(first_one) {
        first_one=false;
        }
    else {
        if((this->start_==-1 || this->vector_generated_called_ >= this->start_) && 
            (this->end_==-1   || this->vector_generated_called_ < this->end_)) {

//          std::cout << "produced ";
//          for(unsigned int m=0; m<vec.size(); ++m)
            //          std::cout << vec[m] << " ";
//          std::cout << std::endl;
            
//          owner_.multiplicity[owner_.current_] *= current_multiplicity;
            owner_.originals.push_back(owner_.originals[owner_.current_]);
            unsigned int siz=owner_.originals.size()-1;

            // Take care of the permutation sign.
            int multiplicity=owner_.multiplicity[owner_.current_] * current_multiplicity;
            if(owner_.permutation_sign==-1)
                multiplicity*=ordersign(vec.begin(), vec.end());
            owner_.multiplicity.push_back(multiplicity);
            
            for(unsigned int k=0; k<owner_.permute_blocks.size(); ++k) {
                for(unsigned int kk=0; kk<owner_.block_length; ++kk) {
                    assert(owner_.permute_blocks[k]+kk<owner_.originals[0].size());
                    owner_.originals[siz][owner_.permute_blocks[k]+kk]=
                        owner_.originals[owner_.current_][owner_.permute_blocks[vec[k]]+kk];
                    }
                }
            }
        }
    }

template<class T>
std::ostream& operator<<(std::ostream& str, const symmetriser<T>& sym)
    {
    for(unsigned int i=0; i<sym.size(); ++i) {
        for(unsigned int j=0; j<sym[i].size(); ++j) {
            str << sym[i][j] << " ";
            }
        str << " ";
        str.setf(std::ios::right, std::ios::adjustfield);
        str << std::setw(2) << sym.signature(i) << std::endl;
        }
    return str;
    }

}