//////////////////////////////////////////////////////////////////////////
//
// ComputeSubsets.hh (unifies ComputeCircuits and ComputeCocircuits)
//
//    produced: 2024/02/29 jr
//
/////////////////////////////////////////////////////////////////////////
#ifndef COMPUTESUBSETS_HH
#define COMPUTESUBSETS_HH

#include <iostream>

#include "Global.hh"
#include "Message.hh"
#include "CommandlineOptions.hh"
#include "LabelSet.hh"
#include "PointConfiguration.hh"
#include "Chirotope.hh"
#include "Symmetry.hh"
#include "Circuits.hh"
#include "Cocircuits.hh"
#include "Permutation.hh"
#include "Symmetry.hh"
#include "SymmetricSubsetGraphMaster.hh"
#include "Signal.hh"

#ifdef TOPCOM_CONTAINERS
#include "PlainHashSet.hh"
#include "HashSet.hh"
namespace topcom {
  typedef PlainHashSet<subset_type> subset_container_type;
};
#else
#include <unordered_set>
namespace topcom {
  typedef std::unordered_set<subset_type, Hash<subset_type> > subset_container_type;
};
#endif

namespace topcom {

  // we want to enumerate objects of type T;
  // the class T has a constructor from subsets with
  // cardinality card_offset compared to rank and a
  // function returning a signature subset for T that
  // uniquely determines an object of type T;
  // moreover, bool T.is_valid() and LabelSet T.support() are defined;
  // alternatively, there is a mode for
  // SymmetricSubsetGraphMaster for T
  // (currently used for circuits and cocircuits enumeration):
  template <class T, parameter_type card_offset>
  class ComputeSubsets {
  public:
    static constexpr int INPUT_CHIRO     =  0x0001;
    static constexpr int OUTPUT_SUBSETS  =  0x0008;
    static constexpr int PREPROCESS      =  0x0010;
  private:
    parameter_type                     _no;
    parameter_type                     _rank;
    PointConfiguration*                _pointsptr;
    Chirotope*                         _chiroptr;
    SymmetryGroup*                     _symmetriesptr;
    SwitchTable<subset_type,lexmin_mode>* _switch_tableptr;
  public:
    bool                  input_chiro;
    bool                  output_subset_objects;
    bool                  preprocess;
  public:
    inline ComputeSubsets() :
      _no(-1),
      _rank(-1),
      _pointsptr(0),
      _chiroptr(0),
      _symmetriesptr(0),
      _switch_tableptr(0),
      input_chiro(false),
      output_subset_objects(false),
      preprocess(false) {}
    inline ComputeSubsets(int flags) :
      _no(-1),
      _rank(-1),
      _pointsptr(0),
      _chiroptr(0),
      _symmetriesptr(0),
      _switch_tableptr(0),
      input_chiro(flags & INPUT_CHIRO),
      output_subset_objects(flags & OUTPUT_SUBSETS),
      preprocess(flags & PREPROCESS) {}
    inline ~ComputeSubsets() {
      if (_pointsptr) {
	delete _pointsptr;
      }
      if (_chiroptr) {
	delete _chiroptr;
      }
      if (_symmetriesptr) {
	delete _symmetriesptr;
      }
      if (_switch_tableptr) {
	delete _switch_tableptr;
      }
    }

    Message& write_header(Message&) const;
    std::istream& read_input(std::istream&);

    int run_from_chiro();
    int run_from_points();

    int  run();
  };

  template <class T, parameter_type card_offset>
  Message& ComputeSubsets<T, card_offset>::write_header(Message& msg) const {
    msg << std::endl;
    msg << "------------------------------------------------------------------\n";
    msg << " computing "
	<< T::type_name()
	<< "s of a point configuration up to symmetry \n";
    msg << " TOPCOM client: " << CommandlineOptions::client() << '\n';
    msg << "------------------------------------------------------------------\n";
    msg << std::endl;
    return msg;
  }

  template <class T, parameter_type card_offset>
  std::istream& ComputeSubsets<T, card_offset>::read_input(std::istream& ist) {
    if (input_chiro) {
      _chiroptr = new Chirotope();
      if (!_chiroptr->read_string(std::cin)) {
	MessageStreams::forced() << "error while reading chirotope - exiting" << std::endl;
	exit(1);
      }
      MessageStreams::verbose() << "read chirotope with " << _chiroptr->no()
				<< " elements in rank " << _chiroptr->rank() << std::endl;
    }
    else {
      _pointsptr = new PointConfiguration();
      if (!_pointsptr->read(ist)) {
	MessageStreams::forced() << "error while reading point configuration - exiting" << std::endl;
	exit(1);
      }
      _pointsptr->transform_to_full_rank();
      MessageStreams::verbose() << "read point configuration with " << _pointsptr->no()
				<< " points in rank " << _pointsptr->rank() << std::endl;
      if ((_pointsptr->no() < 2) || (_pointsptr->rank() < 2)) {
	MessageStreams::forced() << "no of points and rank must be at least two - exiting" << std::endl;
	exit(1);
      }
      if (_pointsptr->rank() > _pointsptr->no()) {
	MessageStreams::forced() << "rank must not be larger than no of points - exiting" << std::endl;
	exit(1);
      }
      _chiroptr = new Chirotope(*_pointsptr, preprocess);
    }
    _no = _chiroptr->no();
    _rank = _chiroptr->rank();
    _symmetriesptr = new SymmetryGroup(_no); // all configurations have a symmetry group
    if (!CommandlineOptions::ignore_symmetries()) {
      if (CommandlineOptions::use_switch_tables()) {
	
	// here, we can get away with generators only:
	if (_symmetriesptr->read_generators(ist)) {
	  MessageStreams::verbose() << "read " << _symmetriesptr->generators().size() << " symmetry generators" << std::endl;
	  MessageStreams::debug() << "symmetry generators:" << '\n'
				  << _symmetriesptr->generators() << std::endl;
	}
	else {
	  MessageStreams::verbose() << "no valid symmetry generators found." << std::endl;
	}
	if (input_chiro) {
	  
	  // here, we need the switch table right away (otherwise it is computed later):
	  MessageStreams::verbose() << "computing switch table ..." << std::endl;
	  _switch_tableptr = new SwitchTable<subset_type, lexmin_mode>(_symmetriesptr->n(), _symmetriesptr->generators());
	  MessageStreams::verbose() << "... done." << std::endl;
	}
      }
      else {
	
	// here, we need a complete symmetry group:
	if (_symmetriesptr->read(ist)) {
	  MessageStreams::verbose() << "read symmetry group with " << _symmetriesptr->generators().size()
				    << " generators of order " << _symmetriesptr->size() + 1 << std::endl;
	  MessageStreams::debug() << "symmetries:" << '\n'
				  << *_symmetriesptr << std::endl;
	}
	else {
	  MessageStreams::verbose() << "no valid symmetry generators found." << std::endl;
	}
      }
    }
    if (CommandlineOptions::preprocess_points()) {
      Symmetry transformation(_pointsptr->no());
      _pointsptr->preprocess(transformation);
      _symmetriesptr->transform(transformation);
    }
    return ist;
  }

  template <class T, parameter_type card_offset>
  int ComputeSubsets<T, card_offset>::run_from_chiro() {
      
    // counters:
    size_type iterations(0);
    size_type count(0);
    size_type total_count(0);
    size_type node_count(0);

    // effort estimation:
    size_type no_of_subsets(global::binomial(_no, _rank + card_offset));
    
    // use either as full storage or cache of already processed sets:
    const bool activate_cache(!CommandlineOptions::memopt() || (CommandlineOptions::localcache() > 0));
    subset_container_type support_sets;
    subset_container_type done_sets;
    if (CommandlineOptions::memopt()) {
      if (activate_cache) {
	done_sets = subset_container_type(std::min<size_type>(CommandlineOptions::localcache(), no_of_subsets));
      }
    }
    else {
      done_sets = subset_container_type(no_of_subsets);
    }
    
    if (output_subset_objects) {
      MessageStreams::result() << _no << ',' << _rank << ":\n{\n";
    }
    
    // use for storage of canonical representatives:
    subset_container_type canonical_support_sets;
    subset_container_type canonical_subsets;
    subset_container_type orbit;
    Permutation subsetperm(_no, _rank + card_offset);
    
    do {
      if (++iterations % CommandlineOptions::report_frequency() == 0) {
	MessageStreams::verbose() << count << " symmetry classes --- "
				  << total_count << " total subset_objects ("
				  << iterations << " out of "
				  << no_of_subsets << " " << _rank + 1 << "-sets investigated [cache load: "
				  << done_sets.size() << ", cache size: "
				  << done_sets.bucket_count() << "])"
				  << std::endl;
      }
      if (done_sets.size() == no_of_subsets) {
	break;
      }
      subset_type subset(subsetperm);
      MessageStreams::debug() << "investigating " << subset << " ..." << std::endl;
      if ((activate_cache) && (done_sets.find(subset) != done_sets.end())) {
	
	// subset already considered and stored/cached:
	continue;
      }
      bool old_subset = false;
      if (CommandlineOptions::use_switch_tables()) {
	if (_switch_tableptr->not_canonical(subset)) {
	  old_subset = true;
	}
      }
      else {
	for (SymmetryGroup::const_iterator iter = _symmetriesptr->begin();
	     iter != _symmetriesptr->end();
	     ++iter) {
	  const Symmetry& g(*iter);
	  
	  // if there is a symmetry mapping this dependent set to something lexicographically smaller,
	  // then we have had an equivalent dependent set already:
	  if (g.lex_decreases(subset)) {
	    MessageStreams::debug() << g.map(subset) << " has been investigated already." << std::endl;
	    old_subset = true;
	    break;
	  }	  
	}
      }
      if (old_subset) {
	
	// nothing new:
	continue;
      }
      
      // store/cache the dependent set as done:
      if (activate_cache) {
	if (CommandlineOptions::memopt()) {
	  if (CommandlineOptions::localcache() > 0) {
	    if (CommandlineOptions::localcache() > done_sets.size()) {
	      done_sets.insert(subset);
	    }
	  }
	}
	else {
	  done_sets.insert(subset);
	}
      }
      
      // store/cache all equivalent dependent sets:
      for (SymmetryGroup::const_iterator iter = _symmetriesptr->begin();
	   iter != _symmetriesptr->end();
	   ++iter) {
	// canonicalize the support set:
	const Symmetry& g(*iter);
	const subset_type equivalent_subset(g.map(subset));
	if (activate_cache) {
	  if (CommandlineOptions::memopt()) {
	    if (CommandlineOptions::localcache() > 0) {
	      if (CommandlineOptions::localcache() > done_sets.size()) {
		done_sets.insert(equivalent_subset);
	      }
	    }
	  }
	  else {
	    done_sets.insert(equivalent_subset);
	  }
	}
      }
      T subset_object(*_chiroptr, subset);
      if (subset_object.is_valid()) {
	MessageStreams::debug() << "computing subset_object support for " << subset << " ..." << std::endl;
	subset_type support_set(subset_object.support());
	bool old_support_set = false;

	// the first subset found is by definition the canonical subset:
	subset_type canonical_support_set = support_set;
	if (support_sets.find(support_set) != support_sets.end()) {
	  
	  // the subset_object support is in the storage/cache:
	  old_support_set = true;
	}
	else if (canonical_support_sets.find(support_set) != canonical_support_sets.end()) {
	  
	  // the subset_object support was found already:
	  old_support_set = true;
	}
	else {
	  for (SymmetryGroup::const_iterator iter = _symmetriesptr->begin();
	       iter != _symmetriesptr->end();
	       ++iter) {
	    const Symmetry& g(*iter);
	    const subset_type equivalent_support_set(g.map(support_set));
	    if (support_sets.find(equivalent_support_set) != support_sets.end()) {

	      // a symmetric image of the subset_object support is in the storage
	      old_support_set = true;
	      break;
	    }
	    else if (canonical_support_sets.find(equivalent_support_set) != canonical_support_sets.end()) {

	      // a symmetric image of the subset_object support was found already
	      old_support_set = true;
	      break;
	    }
	  }
	}
	if (!old_support_set) {
	  if (output_subset_objects) {
	    if (CommandlineOptions::debug()) {
	      MessageStreams::result() << "C[" << count << "] := " << subset_object << ";" << std::endl;
	    }
	    else {
	      MessageStreams::result() << subset_object << std::endl;
	    }
	  }
	  
	  // add the canonicalized support set to the set of already found canonical support sets:
	  canonical_support_sets.insert(canonical_support_set);

	  // insert into support sets:
	  if (activate_cache) {
	    if (CommandlineOptions::memopt()) {
	      if (CommandlineOptions::localcache() > 0) {
		if (CommandlineOptions::localcache() > support_sets.size()) {
		  support_sets.insert(support_set);
		}
	      }
	    }
	    else {
	      support_sets.insert(support_set);
	    }
	  }
	  
	  // we now need to count the number of support sets in the orbit:
	  if (CommandlineOptions::use_switch_tables()) {

	    // stabilizer formula:
	    total_count += _switch_tableptr->orbit_size(support_set);
	  }
	  else {
	    orbit.insert(support_set);
	    for (SymmetryGroup::const_iterator iter = _symmetriesptr->begin();
		 iter != _symmetriesptr->end();
		 ++iter) {
	      
	      // store all equivalent support sets:
	      const Symmetry& g(*iter);
	      const subset_type equivalent_support_set(g.map(support_set));
	      orbit.insert(equivalent_support_set);
	      
	      // store/cache all equivalent support sets as found:
	      if (activate_cache) {
		if (CommandlineOptions::memopt()) {
		  if (CommandlineOptions::localcache() > 0) {
		    if (CommandlineOptions::localcache() > support_sets.size()) {
		      support_sets.insert(equivalent_support_set);
		    }
		  }
		}
		else {
		  support_sets.insert(equivalent_support_set);
		}
	      }
	    }
	    total_count += orbit.size();
	  }
	  count += 1;
	  orbit.clear();
	}
      }
      MessageStreams::debug() << "... done." << std::endl;
    } while (subsetperm.lexnext());
    if (output_subset_objects) {
      MessageStreams::result() << "}" << std::endl;
    }
    MessageStreams::verbose() << count << " symmetry classes";
    if (!CommandlineOptions::skip_orbitcount()) {
      MessageStreams::verbose() << " | " << total_count << " total " << T::type_name() << "s";
    }
    MessageStreams::verbose() << "." << std::endl;
    if (!output_subset_objects) {
      MessageStreams::result() << count << std::endl;
    }
    return 0;
  }
  
  template <>
  int ComputeSubsets<Circuit, 1>::run_from_points() {
    SymmetricSubsetGraphMaster<circuits> ssgm(_pointsptr->no(),
					      _pointsptr->rank(),
					      *_pointsptr,
					      *_symmetriesptr,
					      output_subset_objects,
					      false);
    const Integer uniform_number = field::binomial(Integer(_pointsptr->no()), Integer(_pointsptr->rank() + 1));
    if (Signal::signal_received()) {
      MessageStreams::verbose().print_dumpseparator();
      MessageStreams::verbose() << "### intermediate results at checkpoint forced by signal:" << std::endl;
      MessageStreams::verbose().print_dumpseparator();
    }
    MessageStreams::verbose() << ssgm.symcount() << " symmetry classes";
    if (!CommandlineOptions::skip_orbitcount()) {
      MessageStreams::verbose() << " | " << ssgm.totalcount() << " total circuits";
    }
    MessageStreams::verbose() << " using " << ssgm.nodecount() << " nodes";
    MessageStreams::verbose() << " | " << ssgm.deadendcount() << " branching deadends";
    MessageStreams::verbose() << " | " << ssgm.earlydeadendcount() << " early detected deadends";
    MessageStreams::verbose() << "." << std::endl;
    MessageStreams::verbose() << "(total number is " << ssgm.totalcount() * 100.0 / uniform_number.get_d()
			      << "% of uniform number = " << uniform_number
			      << ")";
    MessageStreams::verbose() << std::endl;
    if (!output_subset_objects) {
      MessageStreams::result() << ssgm.symcount() << std::endl;
    }
    return Signal::exit_value();
  }

  template <>
  int ComputeSubsets<Cocircuit, -1>::run_from_points() {
    SymmetricSubsetGraphMaster<cocircuits> ssgm(_pointsptr->no(),
						_pointsptr->rank(),
						*_pointsptr,
						*_symmetriesptr,
						output_subset_objects,
						false);
    const Integer uniform_number = field::binomial(Integer(_pointsptr->no()), Integer(_pointsptr->rank() + 1));
    if (Signal::signal_received()) {
      MessageStreams::verbose().print_dumpseparator();
      MessageStreams::verbose() << "### intermediate results at checkpoint forced by signal:" << std::endl;
      MessageStreams::verbose().print_dumpseparator();
    }
    MessageStreams::verbose() << ssgm.symcount() << " symmetry classes";
    if (!CommandlineOptions::skip_orbitcount()) {
      MessageStreams::verbose() << " | " << ssgm.totalcount() << " total cocircuits";
    }
    MessageStreams::verbose() << " using " << ssgm.nodecount() << " nodes";
    MessageStreams::verbose() << " | " << ssgm.deadendcount() << " branching deadends";
    MessageStreams::verbose() << " | " << ssgm.earlydeadendcount() << " early detected deadends";
    MessageStreams::verbose() << "." << std::endl;
    MessageStreams::verbose() << "(total number is " << ssgm.totalcount() * 100.0 / uniform_number.get_d()
			      << "% of uniform number = " << uniform_number
			      << ")";
    MessageStreams::verbose() << std::endl;
    if (!output_subset_objects) {
      MessageStreams::result() << ssgm.symcount() << std::endl;
    }
    return Signal::exit_value();
  }
  
  template <class T, parameter_type card_offset>
  int ComputeSubsets<T, card_offset>::run() {
    
    // write header for computations concerning subsets:
    write_header(MessageStreams::forced());
 
    // read input from stdin or file (depending on commandline options):
    read_input(std::cin);

    // get the trivial case out of the way first:
    if (_no == _rank) {
      MessageStreams::verbose() << 0 << " subset_objects in total." << std::endl;
      if (!output_subset_objects) {
	MessageStreams::result() << 0 << std::endl;
      }
      return 0;
    }

    if (input_chiro) {

      // there is no other way currently than to check all rank + 1 subsets mod symmetry:
      return run_from_chiro();
    }
    else {
      
      // if we know the points, we can enumerate subset_objects by symmetric lexicographic subset reverse search:
      return run_from_points();
    }
  }

}; // namespace topcom

#endif

// eof ComputeSubsets.hh