/*=======================================================================*\
 * bimonotone.hh : sorted enumeration for bimonotone functions
 * Copyright (C) 2003-2005 Michael Eisermann, all rights reserved
 * Institut Fourier, Université Grenoble I, France
 *
 * Please send bug reports, suggestions, and comments to
 *   Michael Eisermann <Michael.Eisermann@ujf-grenoble.fr>
 * For more information and updates see the web page
 *   http://www-fourier.ujf-grenoble.fr/~eiserm
 * If you have used my software in your research, please let me know.
 *
 * 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 2 of the License, or 
 * (at your option) any later version. In particular, permission to use, 
 * copy, modify, or redistribute this software for research and education 
 * is granted without fee, provided that the above copyright notice and 
 * this permission appear in all copies and supporting documentation. 
 *
 * This software is distributed in the hope that it will be useful,
 * but it is provided "as is" without any warranty, without even 
 * the implied warranty of fitness for a particular purpose. 
 * See the GNU General Public License for more details.
 *
 *=======================================================================
 * Description:
 *-------------
 * The templates of this file provide generic classes for sorted 
 * enumeration of a bimonotone function, that is, a function  f(a,b)
 * that satisfies  f(a,b) <= f(a',b')  whenever  a <= a'  and  b <= b' .
 * For example, a two-variable polynomial  f(a,b)  is bimonotone if it 
 * has only non-negative coefficients. 
 *
 * The class  BimStream  and  its variant  BimStreamDiff  allow you 
 * to enumerate all values  f(a,b)  in increasing order, which is
 * called an *enumeration stream*. The main point is the efficient 
 * implementation (in time and memory) of the two basic operations:
 * - initialisation to a given level,
 * - advancing in the enumeration stream.
 * 
 * For a documentation of these algorithms and a cost analysis see 
 * Michael Eisermann: Bimonotone enumeration (preprint June 2005)
 *
 *=======================================================================
 * Notation and assumptions:
 *-------------------------- 
 * We wish to implement bimonotone functions f(a,b), where the parameters
 * a and b are of type A and B, respectively, and the result is of type Z.
 * To obtain a neat interface, we specify the following requirements.
 *
 * In the sequel we will assume that A and B implement iterative types,
 * that is, we require incrementation/decrementation as well as comparison:
 *
 * A    operator ++ ( A  a );
 * A    operator -- ( A  a );
 * bool operator == ( A a1, A a2 );
 * bool operator != ( A a1, A a2 );
 *
 * The type Z implements an ordered set, that is, we require comparison 
 * but neither incrementation nor decrementation:
 *
 * bool operator == ( Z z1, Z z2 );
 * bool operator != ( Z z1, Z z2 );
 * bool operator <= ( Z z1, Z z2 );
 * bool operator <  ( Z z1, Z z2 );
 *
 * Optionally, the type D is a difference type for Z, that is, 
 * we require operators of the following form:
 *
 * D    operator +  ( D d1, D d2 );
 * D    operator -  ( D d1, D d2 );
 * bool operator >  ( D d1, D d2 );
 * Z    operator +  ( Z  z, D  d );  // e.g. by converting D to Z
 * Z    operator -  ( Z z1, Z z2 );  // will be converted to D
 *
 * This can be used for a function z = f(a,b) to handle increments as
 * f(a+1,b) = f(a,b) + diff1(a,b)  and  f(a,b+1) = f(a,b) + diff2(a,b).
 * (Of course we will assume that  z+0 = z  and  (z+d)+d' = z+(d+d').)
 *
 *=======================================================================
 * Classes implemented in this file:
 *----------------------------------
 *
 * class BimABtoZ <A,B,Z>
 *   abstract base class for bimonotone functions
 *
 * class BinStream <A,B,Z>
 *   abstract base class for binary enumeration streams
 *
 * class Confluence <A,B,Z>
 *   derived from BinStream <A,B,Z>
 *   implements a confluence of several enumeration streams
 *
 * class BimStream <MapABtoZ>
 *   derived from BinStream <A,B,Z>
 *   implements an enumeration stream for a bimonotone function
 * 
 * class BimStreamDiff <MapABtoZ>
 *   derived from BinStream<A,B,Z>
 *   implements an enumeration stream for a bimonotone function
 *   storing differentials instead of absolute values to save space
 *
 *=======================================================================
 * Implementation and change history (in reverse chronological order)
 *-----------------------------------------------------------------------
 * 2005/06/12, version 0.8, author: Michael Eisermann
 * - Implemented specialized heap operations to test against STL.
 *   (So far this has been done and tested only for BimStreamDiff.)
 *   Remark: When advancing the stream, about 70% of the time is spent 
 *   with pop() and about 10% with push(). This justifies optimization.
 *   Empirically, we can accelerate 10%, which is surprisingly little.
 *-----------------------------------------------------------------------
 * 2005/06/01, version 0.7, author: Michael Eisermann
 * - Implemented class BimABtoZ (base class for bimonotone functions)
 * - Revised, reorganized, and simplified existing templates
 * - Updated and enhanced comments
 *-----------------------------------------------------------------------
 * 2004/08/25, version 0.6, author: Michael Eisermann
 * - Revised and tested class BimStream
 * - Implemented class BimStreamDiff (differential model)
 *-----------------------------------------------------------------------
 * 2003/10/25, version 0.5, author: Michael Eisermann
 * - Implemented class Confluence (merging enumeration streams)
 *-----------------------------------------------------------------------
 * 2003/10/16, version 0.4, author: Michael Eisermann
 * - Some small revisions after extensive testing
 *-----------------------------------------------------------------------
 * 2003/06/01, version 0.3, author: Michael Eisermann
 * - Implemented length counters
 *-----------------------------------------------------------------------
 * 2003/05/16, version 0.2, author: Michael Eisermann
 * - Implemented class BinStream (abstract base class)
 * - Implemented symmetric enumeration (now obsolete)
 *-----------------------------------------------------------------------
 * 2003/03/21, version 0.1, author: Michael Eisermann
 * - Implemented class BimStream (bimonotone enumeration)
\*=======================================================================*/

#ifndef __BIMONOTONE_ENUMERATION__
#define __BIMONOTONE_ENUMERATION__

// Container classes and algorithms from the standard template library:
#include <functional> // generic functionals
#include <algorithm>  // generic algorithms
#include <vector>     // generic vector class
#include <list>       // generic bidirectional list
using namespace std;

// #define _heap      // use own heap implementation
// #undef _heap       // use STL heap algorithms

/*=======================================================================*\
 * class BimABtoZ
 * abstract base class for bimonotone functions.
 *-----------------------------------------------------------------------
 * Provides:
 * Abstract base class for bimonotone functions. The class BimABtoZ is
 * derived from STL's generic class binary_function<A,B,Z> and concrete
 * implementations will serve as function objects for enumeration.
 *-----------------------------------------------------------------------
 * Assumptions:
 * In order to ensure correctness of the enumeration algorithms 
 * implemented below, we have to impose some conditions on f : X -> Z.
 * 
 * Firstly, the domain X where f if defined must be bimonotone,
 * that is, we require X = { (a,b) | a >= h(b) and b >= g(a) }, 
 * where the function h is a monotone map from B to A,
 * and g is a monotone map from A to B, such that
 * h(g(a)) <= a  with equality only for  a = a_min ,
 * g(h(b)) <= b  with equality only for  b = b_min .
 *
 * This means, graphically speaking, that X is bounded by the graphs
 * (a,g(a)) and (h(b),b) of the two auxiliary functions g and h.
 * In particular, X has a smallest element (a_min,b_min);
 * all other points (a,b) satisfy a>=a_min and b>=b_min.
 * For example, the domain { (a,b) | a>=a_min, b>=b_min } 
 * is bimonotone, as is { (a,b) | a_min <= a <= b }.
 * 
 * Secondly, f must be bimonotone, which means  f(a,b) <= f(a',b') 
 * whenever  a <= a'  and  b <= b' . On a bimonotone domain X as above,
 * our function f is bimonotone if and only if f(a,b) <= f(a+1,b) and 
 * f(a,b) <= f(a,b+1) for all pairs (a,b), (a+1,b), (a,b+1) in X.
 *
 * Thirdly, if the differential model is employed, we must require 
 * f(a+1,b) = f(a,b) + diff1(a,b) and f(a,b+1) = f(a,b) + diff2(a,b).
\*=======================================================================*/
template <typename Atype, typename Btype, typename Ztype, typename Dtype= Ztype >
class BimABtoZ : public binary_function< Atype, Btype, Ztype >
{
public:
  // Turn the parameter typenames into public typenames 
  typedef Atype A;
  typedef Btype B;
  typedef Ztype Z;
  typedef Dtype D;

  // Evaluate the function f for parameters (a,b).
  // This function will be assumed to be bimonotone, that is
  // f(a1,b1) <= f(a2,b2) whenever a1 <= a2 and b1 <= b2.
  Z operator() ( A a, B b ) const { return Z(); };
  
  // The domain X of A x B where the function is defined.
  bool domain_contains( A a, B b ) const { return true; };
  A a_min() const { return A(); };
  B b_min() const { return B(); };

  // Define a synomyn for the difference type 
  typedef D difference_type;

  // First differential f(a+1,b) = f(a,b) + diff1(a,b)
  D diff1( A a, B b ) const { return D(); };
  
  // Second differential f(a,b+1) = f(a,b) + diff2(a,b)
  D diff2( A a, B b ) const { return D(); };
};


/*=======================================================================*\
 * class BinStream<A,B,Z>
 * implements a (purely virtual) binary enumeration stream
 *-----------------------------------------------------------------------
 * Provides:
 * Abstract base class defining a binary enumeration stream. It serves
 * as a common interface for derived classes, for example semimonotone 
 * or bimonotone enumeration.
\*=======================================================================*/
template <typename A, typename B, typename Z>
class BinStream
{
public:
  // Public data types
  struct  ZAB { Z value; A first; B second; };
  typedef unsigned long long int width_type;
  typedef unsigned long long int length_type;
  
protected:
  length_type theLength;  
  
public:
  // Public constructor
  BinStream() : theLength(0) {};
  BinStream( const BinStream& s )
    : theLength( s.theLength ) {};
  
  // Destructor should be virtual
  virtual ~BinStream() {};

  // Swap (more efficient than copy)
  virtual void swap( BinStream& s )
    { std::swap( theLength, s.theLength ); };

  // Number of values that have already been output
  length_type length(void) const { return theLength; };
  
  // Reset the enumeration (to a given level)
  virtual void reset( void ) = 0;
  virtual void reset( const Z& level ) = 0;
  
  // Advance in the enumeration
  virtual void advance(void) = 0;

  // Size of internal buffer
  virtual width_type width(void) const = 0;
  
  // Check for end of stream
  virtual bool empty(void) const = 0;
  
  // Read the next value and its arguments
  virtual const Z& value(void)  const = 0;
  virtual const A& first(void)  const = 0;
  virtual const B& second(void) const = 0;
  
  // Alternative forms using operator notation
  const Z& operator * () const { return value(); };
  bool operator ! () const { return empty(); };
  operator bool () const { return !empty(); };
  void operator ++ () { advance(); };

  BinStream& operator >> ( Z& z ) 
    { z= value(); advance(); return *this; };
  BinStream& operator >> ( ZAB& zab )
    { zab.value = value(); zab.first = first(); 
      zab.second= second(); advance(); return *this; };

  /*
  // Comparison of binary enumeration streams
  bool operator == ( const BinStream& s ) const { return value() == s.value(); };
  bool operator != ( const BinStream& s ) const { return value() != s.value(); };
  bool operator <= ( const BinStream& s ) const { return value() <= s.value(); };
  bool operator <  ( const BinStream& s ) const { return value() <  s.value(); };
  bool operator >= ( const BinStream& s ) const { return value() >= s.value(); };
  bool operator >  ( const BinStream& s ) const { return value() >  s.value(); };
  */
};


/*=======================================================================*\
 * class Confluence<A,B,Z>
 * implements a confluence of several enumeration streams
 *-----------------------------------------------------------------------
 * Requires: 
 * Binary enumeration streams e_1,e_2,...,e_n of type BinStream<A,B,Z>,
 * called the tributaries, to be merged into one single stream.
 *-----------------------------------------------------------------------
 * Provides:
 * A binary enumeration stream of type BinStream<A,B,Z> that enumerates
 * the values of its tributaries in increasing order.
\*=======================================================================*/
template <typename A, typename B, typename Z>
class Confluence : public BinStream<A,B,Z>
{
public:
  // Typenames inherited from abstract base class BinStream<A,B,Z>
  typedef typename BinStream<A,B,Z>::width_type  width_type;
  typedef typename BinStream<A,B,Z>::length_type length_type;

protected:
  // List of tributaries and priority queue of values
  typedef BinStream<A,B,Z>  Stype;
  typedef Stype *           Ptype;
  typedef vector<Ptype>     Pvector;
  typedef pair<Z,Ptype>     ZPpair; 
  typedef vector<ZPpair>    ZPqueue;
  
  // Variables to keep track of the tributaries and their current values
  Pvector tributaries;
  ZPqueue values;
  
  // Comparison operator used in the priority queue
  class greater
  {
  public:
    bool operator() ( const ZPpair& x, const ZPpair& y ) const
    { return ( y.first < x.first ); };
  };

  // Push: add a new element to the priority queue.
  virtual void push( const Z& z, const Ptype& p ) 
    { 
      values.push_back( ZPpair( z, p ) ); 
      push_heap( values.begin(), values.end(), greater() );
    };
  
  // Pop: remove the smallest element from the priority queue.
  virtual void pop( void ) 
    { 
      pop_heap( values.begin(), values.end(), greater() );
      values.pop_back();
    };

private:
  // No copy constructor nor copy operator is provided!
  // Notice: it would not suffice to copy the priority queue 
  // and the list of tributaries. We would have to copy each
  // tributary! This approach does not seem reasonable.
  Confluence( const Confluence& s ) {};
  virtual Confluence& operator = ( const Confluence& s );

public:
  // Public constructors (initialization with with 0, 1, 2, or 3 tributaries)
  Confluence() {};
  Confluence( Stype& tribut1 )
    { adjoin( tribut1 ); };
  Confluence( Stype& tribut1, Stype& tribut2 )
    { adjoin( tribut1 ); adjoin( tribut2 ); };
  Confluence( Stype& tribut1, Stype& tribut2, Stype& tribut3 )
    { adjoin( tribut1 ); adjoin( tribut2 ); adjoin( tribut3 ); };

  // Adjoin a tributary
  virtual void adjoin( Stype& tribut )
    {
      tributaries.push_back( &tribut );
      if ( ! tribut.empty() )
	push( tribut.value(), &tribut );
      theLength+= tribut.length();
    };
  
  // Destructor should be virtual
  virtual ~Confluence() {};

  // Swap (more efficient than copy)
  virtual void swap( Confluence& s )
    {
      std::swap( theLength, s.theLength );
      tributaries.swap( s.tributaries );
      values.swap( s.values );
    };
  
  // Reset the enumeration (to a given level)
  virtual void reset(void);
  virtual void reset( const Z& level );
  
  // Advance in the enumeration
  virtual void advance(void);

  // Size of internal buffers (i.e. the sum of tributaries's widt)
  virtual width_type width(void) const 
    { 
      width_type width= 0;
      for ( unsigned int i=0; i<tributaries.size(); ++i )
	width+= tributaries[i]->width();
      return width;
    };

  // Check for end of stream
  virtual bool empty(void) const 
    { return values.empty(); };
  
  // Read the smallest value and its arguments
  virtual const Z& value(void) const
    { return values.front().first; };
  virtual const A& first(void) const
    { return values.front().second->first(); };
  virtual const B& second(void) const
    { return values.front().second->second(); };
};

/*-----------------------------------------------------------------------*\
 * The method reset() resets all tributaties
\*-----------------------------------------------------------------------*/
template <typename A, typename B, typename Z>
void Confluence<A,B,Z>::
reset(void)
{
  theLength= 0; 
  values.clear();
  for ( unsigned int i=0; i<tributaries.size(); ++i )
    {
      tributaries[i]->reset();
      if ( ! tributaries[i]->empty() )
	push( tributaries[i]->value(), tributaries[i] );
    };
};

/*-----------------------------------------------------------------------*\
 * The method reset(level) resets all tributaries to a given level.
\*-----------------------------------------------------------------------*/
template <typename A, typename B, typename Z>
void Confluence<A,B,Z>::
reset( const Z& level )
{
  theLength= 0; 
  values.clear();
  for ( unsigned int i=0; i<tributaries.size(); ++i )
    {
      tributaries[i]->reset(level);
      if ( ! tributaries[i]->empty() )
	push( tributaries[i]->value(), tributaries[i] );
      theLength+= tributaries[i]->length(); 
    };
};

/*-----------------------------------------------------------------------*\
 * The method advance() implements the sorted enumeration algorithm.
\*-----------------------------------------------------------------------*/
template <typename A, typename B, typename Z>
void Confluence<A,B,Z>::
advance(void)
{
  // Security check
  if ( values.empty() ) return;
  
  // Fetch the first value and advance the corresponding stream
  Ptype tribut= values.front().second;
  pop();
  ++theLength;
  tribut->advance();
  if ( ! tribut->empty() )
    push( tribut->value(), tribut );
};


/*=======================================================================*\
 * class BimStream <MapABtoZ>
 * implements a bimonotone enumeration stream
 *-----------------------------------------------------------------------
 * Requires:
 * The class MapABtoZ must be derived from the abstract class BimABtoZ.
 * It behaves as a function object, implementing a bimonotone function f
 * mapping some bimonotone domain X of A x B to an ordered set Z.
 * For details see the the abstract base class BimABtoZ above.
 *-----------------------------------------------------------------------
 * Provides:
 * A stream enumerating all pairs (a,b) in X such that the values f(a,b) 
 * appear in increasing order.
\*=======================================================================*/
template <class MapABtoZ>
class BimStream : public BinStream< typename MapABtoZ::A,
		                    typename MapABtoZ::B,
		                    typename MapABtoZ::Z >
{
public:
  // Typenames provided by MapABtoZ
  typedef typename MapABtoZ::A A;
  typedef typename MapABtoZ::B B;
  typedef typename MapABtoZ::Z Z;

  // Typenames inherited from BinStream
  typedef typename BinStream<A,B,Z>::width_type  width_type;
  typedef typename BinStream<A,B,Z>::length_type length_type;

protected:
  // List of arguments and priority queue of values
  typedef pair<A,B>                 ABpair;
  typedef list<ABpair>              ABlist;
  typedef typename ABlist::iterator ABiter;
  typedef pair<Z,ABiter>            ZIpair;
  typedef vector<ZIpair>            ZIqueue;
  
protected:
  // Internal states of the enumeration stream
  MapABtoZ theFunction;
  ABlist   minima;
  ZIqueue  values;

protected:
  // Initial construction of the list and the priority queue.
  // (Assuming both to be empty.)
  virtual void init()
    {
      theLength= 0;
      A a0= theFunction.a_min();
      B b0= theFunction.b_min();
      minima.push_front( ABpair( a0, b0 ) );
      push( theFunction( a0, b0 ), minima.begin() );
    };

  // Comparison operator used in the priority queue
  class greater
  {
  public:
    bool operator() ( const ZIpair& x, const ZIpair& y ) const
    { return ( y.first < x.first ); };
  };
  
  // Push: add a new element to the priority queue.
  virtual void push( const Z& z, const ABiter& it ) 
    { 
      values.push_back( ZIpair( z, it ) ); 
      push_heap( values.begin(), values.end(), greater() );
    };
  
  // Pop: remove the smallest element from the priority queue.
  virtual void pop( void ) 
    { 
      pop_heap( values.begin(), values.end(), greater() );
      values.pop_back();
    };
      
public:
  // Public constructor using the data provided by the function f
  BimStream( const MapABtoZ& f= MapABtoZ() ) 
    : theFunction( f ) { init(); };
  
  // Copy constructor: duplication cannot be left to the default 
  // operator since the queue stores iterators of the list.
  BimStream( const BimStream& s )
    : theFunction( s.theFunction ), minima( s.minima )
    {
      typename ZIqueue::const_iterator jt= s.values.begin();
      for( ABiter it= minima.begin(); it != minima.end(); ++it, ++jt )
	push( jt->first, it );
    };
  
  // Destructor should be virtual
  virtual ~BimStream() {};
  
  // Copy operator: duplication as explained above
  virtual BimStream& operator = ( const BimStream& s );
  
  // Swap (more efficient than copy)
  virtual void swap( BimStream& s )
    {
      std::swap( theLength,   s.theLength   );
      std::swap( theFunction, s.theFunction );
      minima.swap( s.minima );
      values.swap( s.values );
    };

  // Reset the enumeration (to a given level)
  virtual void reset( void );
  virtual void reset( const Z& level );
  
  // Advance in the enumeration
  virtual void advance(void);

  // Size of internal buffer (i.e. of the priority queue)
  virtual width_type width(void) const 
    { return values.size(); };
  
  // Check for end of stream
  virtual bool empty(void) const 
    { return values.empty(); };
  
  // Read the smallest value and its arguments
  virtual const Z& value(void) const
    { return values.front().first; };
  virtual const A& first(void) const
    { return values.front().second->first; };
  virtual const B& second(void) const
    { return values.front().second->second; };
};

/*-----------------------------------------------------------------------*\
 * Tailor-made copy operator that duplicates the list of minima together
 * with the priority queue of values containing pointers to the list.
\*-----------------------------------------------------------------------*/
template <class MapABtoZ>
BimStream<MapABtoZ>& BimStream<MapABtoZ>::
operator = ( const BimStream& s )
{ 
  theLength   = s.theLength; 
  theFunction = s.theFunction;
  minima      = s.minima;
  values.clear();
  typename ZIqueue::const_iterator jt= s.values.begin();
  for ( ABiter it= minima.begin(); it != minima.end(); ++it, ++jt )
    push( jt->first, it );
  return *this; 
};

/*-----------------------------------------------------------------------*\
 * The method reset() resets the enumeration stream.
\*-----------------------------------------------------------------------*/
template <class MapABtoZ>
void BimStream<MapABtoZ>::
reset(void)
{
  minima.clear();
  values.clear();
  init();
};

/*=======================================================================*\
 * method reset(level)
 * resets the enumeration stream to a given level.
 *-----------------------------------------------------------------------
 * Requires:
 * In order to ensure correctness of this algorithm, we have to require
 * that the function f be bimonotone and defined on a bimonotone domain X.
 * For details see the the abstract base class BimABtoZ above.
 *-----------------------------------------------------------------------
 * Provides:
 * An ordered list of the minimal elements of { x in X | f(x) >= level },
 * together with a priority queue of their respective values f(x).
 * These data are stored in the member objects  minima  and  values.
 * Thus initialized the stream can enumerate all values >= level
 * by means of the sorted enumeration algorithm below.
\*=======================================================================*/
template <class MapABtoZ>
void BimStream<MapABtoZ>::
reset( const Z& level )
{
  minima.clear();
  values.clear();

  A a= theFunction.a_min();  // initial value for a
  B b= theFunction.b_min();  // initial value for b
  A a0= theFunction.a_min(); // minimal such that (a0,b) is in domain.
  theLength= 0;              // counts pairs (a,b) with f(a,b) < level
  length_type lineLength= 0; // counts pairs from (a0,b) to (a,b)

  // Find the end of the minima list
  for(;;)
    {
      // Increment a in order to arrive at  f(a,b) >= level
      while( theFunction.domain_contains(a,b) && theFunction(a,b) < level ) 
	{ ++a; ++lineLength; };
      if ( theFunction.domain_contains(a,b) ) break;

      // If we have run out of the domain, step back and increment b.
      --a; 
      ++b; 
      theLength+= lineLength;
      --lineLength;

      // If (a,b) is not in the domain either, then we have finished.
      if ( !theFunction.domain_contains(a,b) ) return;

      // Otherwise find the beginning of the line to update its length
      while( !theFunction.domain_contains(a0,b) )
	{ ++a0; --lineLength; };
    };

  // At this point we have found a point (a,b) in the domain
  // such that f(a,b) >= level and b is minimal with this property.
  for(;;)
    {
      // We do not know whether a is minimal, so let's minimize it.
      do { --a; --lineLength; }
      while( theFunction.domain_contains(a,b) && level <= theFunction(a,b) );
      ++a; ++lineLength;

      // At this point we know that (a,b) is a minimal element
      // with respect to the property f(a,b) >= level.
      minima.push_front( ABpair(a,b) );
      push( theFunction(a,b), minima.begin() );

      // Decrement a once, and start incrementing b...
      --a;
      do
	{
	  // Increment b in order to arrive at f(a,b) >= level.
	  ++b; theLength+= lineLength;

	  // If we have run out of the domain, then we have finished.
	  if ( !theFunction.domain_contains(a,b) ) return;

	  // Otherwise find the beginning of the line to update its length
	  while( !theFunction.domain_contains(a0,b) )
	    { ++a0; --lineLength; };
	}
      while( theFunction(a,b) < level );
      --lineLength;
    };
};

/*=======================================================================*\
 * method advance()
 * implements the sorted enumeration algorithm.
 *-----------------------------------------------------------------------
 * Requires:
 * In order to ensure correctness of this algorithm, we have to require
 * that the function f be bimonotone and defined on a bimonotone domain X.
 * For details see the the abstract base class BimABtoZ above.
 *-----------------------------------------------------------------------
 * Provides:
 * Deletes the smallest value f(a,b) from the priority queue and updates 
 * the list of minima and the queue of values. To do so, it suffices to
 * consider two possible new minima, (a+1,b) and (a,b+1), and to insert
 * them into the list and queue if necessary.
\*=======================================================================*/
template <class MapABtoZ>
void BimStream<MapABtoZ>::
advance(void)
{
  // Security check
  if ( values.empty() ) return;
  
  // Read the pair (a,b) that realizes the current minimum
  ABiter pos= values.front().second;
  A a( pos->first );
  B b( pos->second );
  
  // Should the first successor (a+1,b) be inserted?
  A a1(a); ++a1;
  ABiter succ(pos); ++succ;
  bool ins1= ( succ == minima.end() ? 
	       theFunction.domain_contains(a1,b) : 
	       ( a1 != succ->first ) );
  
  // Should the second successor (a,b+1) be inserted?
  B b1(b); ++b1;
  ABiter pred(pos); --pred;
  bool ins2 = ( pos == minima.begin() ?
		theFunction.domain_contains(a,b1) :
		( b1 != pred->second ) );
  
  // Update the list of minimal elements and the queue of minimal values
  pop();
  ++theLength;
  if ( ins1 )
    {
      if ( ins2 )
	{
	  pos->first= a1;
	  push( theFunction(a1,b), pos );
	  pos= minima.insert( pos, ABpair(a,b1) );
	  push( theFunction(a,b1), pos );
	}
      else
	{
	  pos->first= a1;
	  push( theFunction(a1,b), pos );
	}
    }
  else
    {
      if ( ins2 )
	{
	  pos->second= b1;
	  push( theFunction(a,b1), pos );
	}
      else
	{
	  minima.erase(pos);
	};
    };
};


/*=======================================================================*\
 * class BimStreamDiff <MapABtoZ>
 * implements a bimonotone enumeration stream using differentials
 * (storing incremental values rather than absolute values)
 *-----------------------------------------------------------------------
 * Requires:
 * The class MapABtoZ must be derived from the abstract class BimABtoZ.
 * It behaves as a function object, implementing a bimonotone function f
 * mapping some bimonotone domain X of A x B to an ordered set Z.
 * For details see the the abstract base class BimABtoZ above.
 *
 * For the differential model we additionally employ a type named D,
 * which can be added to Z, and member functions diff1() and diff(2) 
 * that calculate differentials  f(a+1,b) = f(a,b) + diff1(a,b)  and  
 * f(a,b+1) = f(a,b) + diff2(a,b). 
 *-----------------------------------------------------------------------
 * Provides:
 * A stream enumerating all pairs (a,b) in X such that the values f(a,b) 
 * appear in increasing order.
 *-----------------------------------------------------------------------
 * Notice:
 * Only the differentials are stored in a priority queue. This allows us 
 * to choose a type of smaller size, in order to save space and time.
 * Overflow is avoided by observing a maximal tolerance: when the first
 * differential (at the front of the priority queue) is larger than the
 * prescribed tolerance, then all differentials in the queue are reduced
 * and the offset is updated. To be efficient, one chooses the tolerance
 * as large as possible, but still small enough to avoid overflow.
\*=======================================================================*/
template <class MapABtoZ> 
class BimStreamDiff : public BinStream< typename MapABtoZ::A,
		                        typename MapABtoZ::B,
		                        typename MapABtoZ::Z >
{
public:
  // Typenames provided by MapABtoZ
  typedef typename MapABtoZ::A A;
  typedef typename MapABtoZ::B B;
  typedef typename MapABtoZ::Z Z;
  typedef typename MapABtoZ::D D;

  // Typenames inherited from BinStream
  typedef typename BinStream<A,B,Z>::width_type  width_type;
  typedef typename BinStream<A,B,Z>::length_type length_type;

protected:
  // List of arguments and priority queue of values
  typedef pair<A,B>                 ABpair;
  typedef list<ABpair>              ABlist;
  typedef typename ABlist::iterator ABiter;
  typedef pair<D,ABiter>            DIpair;
  typedef vector<DIpair>            DIqueue;
  
protected:
  // Internal states of the enumeration stream
  MapABtoZ theFunction;
  ABlist   minima;
  DIqueue  values;
  Z        minval;
  Z        offset;
  D        tolerance;

protected:
  // Initial construction of the list and the priority queue.
  // (Assuming both to be empty.)
  virtual void init()
    {
      theLength= 0; 
      A a0= theFunction.a_min();
      B b0= theFunction.b_min();
      minval= theFunction( a0, b0 );
      offset= minval;
      minima.push_front( ABpair( a0, b0 ) );
      push( D(0), minima.begin() );
    };

#ifndef _heap
  // Comparison operator used in the priority queue
  class greater
  {
  public:
    bool operator() ( const DIpair& x, const DIpair& y ) const
      { return ( y.first < x.first ); };
  };
#endif
  
  // Push: add a new element to the priority queue.
  virtual void push( const D& d, const ABiter& it ) 
    {
#ifndef _heap
      // push_heap: STL version (generic but slightly slower)
      values.push_back( DIpair( d, it ) ); 
      push_heap( values.begin(), values.end(), greater() );
      minval= offset + values.front().first;
      return;
#else
      // push_heap: specialized version for optimization & comparison
      typedef int Index;
      const Index size= values.size();
      if( size == 0 ) { values.push_back( DIpair(d,it) ); return; }
      
      values.resize( size + 1 );
      DIpair* const v= &(values.front());
      DIpair* i= v + size;
      Index predecessor= (size-1) >> 1;
      for(;;)
	{
	  DIpair* p= v + predecessor;
	  if ( p->first <= d ) break;
	  *i= *p; i= p;
	  if( predecessor == 0 ) break;
	  else predecessor= (predecessor-1) >> 1;
	};
      i->first= d;
      i->second= it;
      minval= offset + v->first;
      return;
#endif
    };
  
  // Pop: remove the smallest element from the priority queue.
  virtual void pop( void ) 
    {
#ifndef _heap
      // pop_heap: STL version (generic and slightly slower)
      pop_heap( values.begin(), values.end(), greater() );
      values.pop_back();
      if ( values.empty() ) return;
      D base= values.front().first;
      minval= offset + base;
      if ( base > tolerance )
	{
	  typename DIqueue::iterator it;
	  for ( it= values.begin(); it != values.end(); ++it )
	    (it->first)-= base; 
	  offset+= base;
	};
      return;
#else
      // pop_heap: specialized version for optimization & comparison
      typedef int Index;
      const Index size= values.size()-1;
      if ( size == 0 ) { values.clear(); return; }
      
      DIpair* const v= &(values.front());
      const DIpair* const j= v+size;
      DIpair* i= v;
      Index successor= 1;
      while( successor < size )
	{
	  DIpair* s= v + successor;
	  if( successor < size-1 )
	    if( s->first > (s+1)->first )
	      { ++successor; ++s; };
	  if( j->first > s->first )
	    {
	      *i= *s; i= s;
	      successor= 2 * successor + 1;
	    }
	  else break;
	};
      *i= *j;
      values.pop_back();

      // Reduce differentials if tolerance is exceeded.
      const D base= v->first;
      minval= offset + base;
      if ( base > tolerance )
	{
	  for( Index i=0; i<size; ++i )
	    (v+i)->first -= base; 
	  offset+= base;
	};
      return;
#endif      
    };
  
public:
  // Public constructor (requiring tolerance only)
  BimStreamDiff( const D& t )
    : tolerance(t) { init(); };

  // Public constructor (requiring function and tolerance)
  BimStreamDiff( const MapABtoZ& f, const D& t )
    : theFunction(f), tolerance(t) { init(); };

  // Copy constructor: duplication cannot be left to the default 
  // operator since the queue stores iterators of the list.
  BimStreamDiff( const BimStreamDiff& s )
    : theFunction( s.theFunction ),
      minima( s.minima ), minval( s.minval ),
      offset( s.offset ), tolerance( s.tolerance )
    {
      typename DIqueue::const_iterator jt= s.values.begin();
      for( ABiter it= minima.begin(); it != minima.end(); ++it, ++jt )
	push( jt->first, it );
    };
  
  // Destructor should be virtual
  virtual ~BimStreamDiff() {};
  
  // Copy operator: duplication as explained above
  virtual BimStreamDiff& operator = ( const BimStreamDiff& s );

  // Swap (more efficient than copy)
  virtual void swap( BimStreamDiff& s )
    {
      std::swap( theLength,   s.theLength   );
      std::swap( theFunction, s.theFunction );
      std::swap( tolerance,   s.tolerance   );
      std::swap( minval,      s.minval      );
      std::swap( offset,      s.offset      );
      minima.swap( s.minima );
      values.swap( s.values );
    };

  // Functions to manipulate the tolerated difference
  virtual void set_tolerance( const D& t ) { tolerance= t; };
  virtual D get_tolerance(void) const      { return tolerance; };

  // Reset the enumeration (to a given level)
  virtual void reset( void );
  virtual void reset( const Z& level );
  
  // Advance in the enumeration
  virtual void advance(void);

  // Size of internal buffer (i.e. of the priority queue)
  virtual width_type width(void) const 
    { return values.size(); };
  
  // Check for end of stream
  virtual bool empty(void) const 
    { return values.empty(); };
  
  // Read the smallest value and its arguments
  virtual const Z& value(void) const
    { return minval; };
  virtual const A& first(void) const
    { return values.front().second->first; };
  virtual const B& second(void) const
    { return values.front().second->second; };
};

/*-----------------------------------------------------------------------*\
 * Tailor-made copy operator that duplicates the list of minima together
 * with the priority queue of values containing pointers to the list.
\*-----------------------------------------------------------------------*/
template <class MapABtoZ>
BimStreamDiff<MapABtoZ>& BimStreamDiff<MapABtoZ>::
operator = ( const BimStreamDiff& s )
{ 
  theLength   = s.theLength; 
  theFunction = s.theFunction;
  minima      = s.minima;
  minval      = s.minval;
  offset      = s.offset;
  tolerance   = s.tolerance;
  values.clear();
  typename DIqueue::const_iterator jt= s.values.begin();
  for ( ABiter it= minima.begin(); it != minima.end(); ++it, ++jt )
    push( jt->first, it );
  return *this; 
};

/*-----------------------------------------------------------------------*\
 * The method reset() resets the enumeration stream.
\*-----------------------------------------------------------------------*/
template <class MapABtoZ>
void BimStreamDiff<MapABtoZ>::
reset(void)
{
  minima.clear();
  values.clear();
  init();
};
  
/*=======================================================================*\
 * method reset(level)
 * resets the enumeration stream to a given level.
 *-----------------------------------------------------------------------
 * Requires:
 * In order to ensure correctness of this algorithm, we have to require
 * that the function f be bimonotone and defined on a bimonotone domain X.
 * For details see the the abstract base class BimABtoZ above.
 *-----------------------------------------------------------------------
 * Provides:
 * An ordered list of the minimal elements of { x in X | f(x) >= level },
 * together with a priority queue of their respective values f(x).
 * These data are stored in the member objects  minima  and  values.
 * (Notice that, in the differential model, we do not store the absolute
 * values, but rather their differentials, based on some given offset.)
 * Thus initialized the stream can enumerate all values >= level
 * by means of the sorted enumeration algorithm below.
\*=======================================================================*/
template <class MapABtoZ>
void BimStreamDiff<MapABtoZ>::
reset( const Z& level )
{
  minima.clear();
  values.clear();
  minval= level;
  offset= level;

  A a= theFunction.a_min();  // initial value for a
  B b= theFunction.b_min();  // initial value for b
  A a0= theFunction.a_min(); // minimal such that (a0,b) is in domain.
  theLength= 0;              // counts pairs (a,b) with f(a,b) < level
  length_type lineLength= 0; // counts pairs from (a0,b) to (a,b)

  // Find the end of the minima list
  for(;;)
    {
      // Increment a in order to arrive at  f(a,b) >= level
      while( theFunction.domain_contains(a,b) && theFunction(a,b) < level ) 
	{ ++a; ++lineLength; };
      if ( theFunction.domain_contains(a,b) ) break;

      // If we have run out of the domain, step back and increment b.
      --a; 
      ++b; 
      theLength+= lineLength;
      --lineLength;

      // If (a,b) is not in the domain either, then we have finished.
      if ( !theFunction.domain_contains(a,b) ) return;

      // Otherwise find the beginning of the line to update its length
      while( !theFunction.domain_contains(a0,b) )
	{ ++a0; --lineLength; };
    };

  // At this point we have found a point (a,b) in the domain
  // such that f(a,b) >= level and b is minimal with this property.
  for(;;)
    {
      // We do not know whether a is minimal, so let's minimize it.
      do { --a; --lineLength; }
      while( theFunction.domain_contains(a,b) && level <= theFunction(a,b) );
      ++a; ++lineLength;

      // At this point we know that (a,b) is a minimal element
      // with respect to the property f(a,b) >= level.
      minima.push_front( ABpair(a,b) );
      push( theFunction(a,b)-offset, minima.begin() );

      // Decrement a once, and start incrementing b...
      --a;
      do
	{
	  // Increment b in order to arrive at f(a,b) >= level.
	  ++b; theLength+= lineLength;

	  // If we have run out of the domain, then we have finished.
	  if ( !theFunction.domain_contains(a,b) ) return;

	  // Otherwise find the beginning of the line to update its length
	  while( !theFunction.domain_contains(a0,b) )
	    { ++a0; --lineLength; };
	}
      while( theFunction(a,b) < level );
      --lineLength;
    };
};


/*=======================================================================*\
 * method advance()
 * implements the sorted enumeration algorithm.
 *-----------------------------------------------------------------------
 * Requires:
 * In order to ensure correctness of this algorithm, we have to require
 * that the function f be bimonotone and defined on a bimonotone domain X.
 * For details see the the abstract base class BimABtoZ above.
 *-----------------------------------------------------------------------
 * Provides:
 * Deletes the smallest value f(a,b) from the priority queue and updates 
 * the list of minima and the queue of values. To do so, it suffices to
 * consider two possible new minima, (a+1,b) and (a,b+1), and to insert
 * them into the list and queue if necessary.
 * (Notice that, in the differential model, we do not store the absolute
 * values, but rather their differentials, based on some given offset.)
\*=======================================================================*/
template <class MapABtoZ>
void BimStreamDiff<MapABtoZ>::
advance(void)
{
  // Security check
  if ( values.empty() ) return;
  
  // Read the pair (a,b) that realizes the current minimum
  ABiter pos= values.front().second;
  A a( pos->first );
  B b( pos->second );
  
  // Should the first successor (a+1,b) be inserted?
  A a1(a); ++a1;
  ABiter succ(pos); ++succ;
  bool ins1= ( succ == minima.end() ? 
	       theFunction.domain_contains(a1,b) : 
	       ( a1 != succ->first ) );
  
  // Should the second successor (a,b+1) be inserted?
  B b1(b); ++b1;
  ABiter pred(pos); --pred;
  bool ins2 = ( pos == minima.begin() ?
		theFunction.domain_contains(a,b1) :
		( b1 != pred->second ) );
  
  // Update the list of minimal elements and the queue of minimal values
  if ( ins1 )
    {
      if ( ins2 )
	{
	  pos->first= a1;
	  push( values.front().first + theFunction.diff1(a,b), pos );
	  pos= minima.insert( pos, ABpair(a,b1) );
	  push( values.front().first + theFunction.diff2(a,b), pos );
	}
      else
	{
	  pos->first= a1;
	  push( values.front().first + theFunction.diff1(a,b), pos );
	}
    }
  else
    {
      if ( ins2 )
	{
	  pos->second= b1;
	  push( values.front().first + theFunction.diff2(a,b), pos );
	}
      else
	{
	  minima.erase(pos);
	};
    };
  
  // It is important to delete the old value only at the very end, because 
  // it serves as reference for incrementation by diff1 and diff2 above.
  pop();
  ++theLength;
};

/*=======================================================================*/

#endif /* __MONOTONE_ENUMERATION__ */

/*=======================================================================*\
 * end of file "bimonotone.hh"
\*=======================================================================*/