#ifndef MB_MESHSET_HPP
#define MB_MESHSET_HPP

#ifndef IS_BUILDING_MB
#error "MB_MeshSet.hpp isn't supposed to be included into an application"
#endif

#include "moab_export.h"

#include "moab/Interface.hpp"
#include "Internals.hpp"
#include "moab/Range.hpp"
#include "moab/CN.hpp"

#include <assert.h>
#include <vector>
#include <algorithm>
#include <iterator>

namespace moab {

class AEntityFactory;

/** \brief Class to implement entity set functionality 
  * \author Jason Kraftcheck <kraftche@cae.wisc.edu>
  */
class MOAB_EXPORT MeshSet
{
public:

  //! create an empty meshset
  inline MeshSet();
  inline MeshSet(unsigned flags);

  //! destructor
  inline ~MeshSet();
  
  inline ErrorCode set_flags( unsigned flags, EntityHandle my_handle, AEntityFactory* adjacencies );
    

    //! get all children pointed to by this meshset
  inline const EntityHandle* get_children( int& count_out ) const ;

    //! get all parents pointed to by this meshset
  inline const EntityHandle* get_parents( int& count_out ) const;

    //! return the number of children pointed to by this meshset
  inline int num_children() const ;
    
    //! return the number of parents pointed to by this meshset
  inline int num_parents() const;

    //! add a parent to this meshset; returns true if parent was added, 0 if it was
    //! already a parent of this meshset
  int add_parent(EntityHandle parent);
    
    //! add a child to this meshset; returns true if child was added, 0 if it was
    //! already a child of this meshset
  int add_child(EntityHandle child);
    
    //! remove a parent from this meshset; returns true if parent was removed, 0 if it was
    //! not a parent of this meshset
  int remove_parent(EntityHandle parent);
    
    //! remove a child from this meshset; returns true if child was removed, 0 if it was
    //! not a child of this meshset
  int remove_child(EntityHandle child);

  unsigned flags() const { return mFlags; }
  //! returns whether entities of meshsets know this meshset 
  int tracking()     const { return mFlags & MESHSET_TRACK_OWNER; }
  int set()          const { return mFlags & MESHSET_SET; }
  int ordered()      const { return mFlags & MESHSET_ORDERED; }
  int vector_based() const { return ordered(); }

    //! replace one entity with another in the set (contents and parent/child
    //! lists); returns whether it was replaced or not
  ErrorCode replace_entities(EntityHandle my_handle, 
                               const EntityHandle *old_entities, 
                               const EntityHandle *new_entities,
                               size_t num_entities,
                               AEntityFactory* mAdjFact);
  
    /** Clear *contents* of set (not parents or children) */
  inline ErrorCode clear( EntityHandle myhandle, AEntityFactory* adjacencies );
  
    /** Clear all set lists (contents, parents, and children) */
  inline ErrorCode clear_all( EntityHandle myhandle, AEntityFactory* adjacencies );

    /** Get contents data array.  NOTE: this may not contain what you expect if not vector_based */
  inline const EntityHandle* get_contents( size_t& count_out ) const;
    /** Get contents data array.  NOTE: this may not contain what you expect if not vector_based */
  inline EntityHandle* get_contents( size_t& count_out );

    /** Get entities contained in set */
  inline ErrorCode get_entities(std::vector<EntityHandle>& entities) const;

    /** Get entities contained in set */
  inline ErrorCode get_entities( Range& entities ) const;
  
    //! get all entities in this MeshSet with the specified type
  inline ErrorCode get_entities_by_type(EntityType entity_type, std::vector<EntityHandle> &entity_list) const;
  
  inline ErrorCode get_entities_by_type( EntityType type, Range& entity_list) const;
    
    //! return the number of entities with the given type contained in this meshset
  inline unsigned int num_entities_by_type(EntityType type) const;
      
  inline ErrorCode get_entities_by_dimension( int dimension, std::vector<EntityHandle> &entity_list) const;

  inline ErrorCode get_entities_by_dimension( int dimension, Range& entity_list) const;
  
    //! return the number of entities with the given type contained in this meshset
  inline unsigned int num_entities_by_dimension(int dimension) const;
      
  inline ErrorCode get_non_set_entities( Range& range ) const;

  /** Test of meshset contains some or all of passed entities
   *
   *\param entities Array of entities
   *\param num_entities Length of array of entities.
   *\param op - Interface::UNION     : Test if set contains any of the input entities
   *          - Interface::INTERSECT : Test if set contains all of the input entities
   */
  inline bool contains_entities(const EntityHandle *entities, int num_entities, const int op) const;

  
  //! subtract/intersect/unite meshset_2 from/with/into meshset_1; modifies meshset_1
  inline ErrorCode subtract(const MeshSet *meshset_2,
                               EntityHandle my_handle,
                               AEntityFactory* adjacencies);

  ErrorCode intersect(const MeshSet *meshset_2,
                                EntityHandle my_handle,
                                AEntityFactory* adjacencies);

  inline ErrorCode unite(const MeshSet *meshset_2,
                            EntityHandle my_handle,
                            AEntityFactory* adjacencies);

  //! add these entities to this meshset
  inline ErrorCode add_entities(const EntityHandle *entity_handles,
                                   const int num_entities,
                                   EntityHandle my_handle,
                                   AEntityFactory* adjacencies);
    
    //! add these entities to this meshset
  inline ErrorCode add_entities(const Range &entities,
                                   EntityHandle my_handle,
                                   AEntityFactory* adjacencies);
    
    //! add these entities to this meshset
  inline ErrorCode remove_entities(const Range& entities,
                                      EntityHandle my_handle,
                                      AEntityFactory* adjacencies);
    
   
    //! remove these entities from this meshset
  inline ErrorCode remove_entities(const EntityHandle *entities,
                                      const int num_entities,
                                      EntityHandle my_handle,
                                      AEntityFactory* adjacencies);

    //! return the number of entities contained in this meshset
  inline unsigned int num_entities() const;
  
  inline bool empty() const { return mContentCount == ZERO; }
  
  unsigned long get_memory_use() const;

protected:
  
  /** Convert for changing flag values */
  ErrorCode convert( unsigned flags, EntityHandle my_handle, AEntityFactory* adj );

  /** Add explicit adjacencies from all contained entities to this (i.e. convert to tracking) */
  ErrorCode create_adjacencies( EntityHandle myhandle, AEntityFactory* adjacencies );

  /** Remvoe explicit adjacencies from all contained entities to this (i.e. convert from tracking) */
  ErrorCode remove_adjacencies( EntityHandle myhandle, AEntityFactory* adjacencies );

  /** Insert vector of handles into MeshSet */
  ErrorCode insert_entity_vector( const EntityHandle* vect, size_t len, EntityHandle my_h, AEntityFactory* adj );
 
  /** Insert vector of handle range pairs into MeshSet */
  ErrorCode insert_entity_ranges( const EntityHandle* range_vect, size_t len, EntityHandle my_h, AEntityFactory* adj );

  /** Insert Range of handles into MeshSet */
  ErrorCode insert_entity_ranges( const Range& range, EntityHandle my_h, AEntityFactory* adj );

  /** Remove vector of handles from MeshSet */
  ErrorCode remove_entity_vector( const EntityHandle* vect, size_t len, EntityHandle my_h, AEntityFactory* adj );

  /** Remove vector of handle range pairs from MeshSet */
  ErrorCode remove_entity_ranges( const EntityHandle* range_vect, size_t len, EntityHandle my_h, AEntityFactory* adj );

  /** Remove Range of handles from MeshSet */
  ErrorCode remove_entity_ranges( const Range& range, EntityHandle my_h, AEntityFactory* adj );

public:  
    //! Possible values of mParentCount and mChildCount
  enum Count { ZERO=0, ONE=1, TWO=2, MANY=3 };
    //! If the number of entities is less than 3, store
    //! the handles directly in the hnd member.  Otherwise
    //! use the ptr member to hold the beginning and end
    //! of a dynamically allocated array.
  union CompactList {
    EntityHandle hnd[2];  //!< Two handles
    EntityHandle* ptr[2]; //!< begin and end pointers for array
  };

private:
  //!Meshset propery flags
  unsigned char mFlags;
  //! If less than MANY, the number of parents stored inline in
  //! parentMeshSets.hnd.  If MANY, then parentMeshSets.ptr contains
  //! array begin and end pointers for a dynamically allocated array
  //! of parent handles.
  unsigned mParentCount : 2;
  //! If less than MANY, the number of children stored inline in
  //! childMeshSets.hnd.  If MANY, then childMeshSets.ptr contains
  //! array begin and end pointers for a dynamically allocated array
  //! of child handles.
  unsigned mChildCount : 2;
  //! If less than MANY, the number of children stored inline in
  //! contentList.hnd.  If MANY, then contentList.ptr contains
  //! array begin and end pointers for a dynamically allocated array..
  unsigned mContentCount : 2;
  //! Storage for data lists
  CompactList parentMeshSets, childMeshSets, contentList;

public:
    /** get dimension of enity */
  static inline int DIM_FROM_HANDLE( EntityHandle h )
    { return CN::Dimension( TYPE_FROM_HANDLE( h ) ); }
  
    /** Get smallest possible handle with specified dimension (first handle for first type of dimension) */
  static inline EntityHandle FIRST_OF_DIM( int dim )
    { return FIRST_HANDLE( CN::TypeDimensionMap[dim].first ); }
  
    /** Get largest possible handle with specified dimension (largest handle for last type of dimension) */
  static inline EntityHandle LAST_OF_DIM( int dim )
    { return LAST_HANDLE( CN::TypeDimensionMap[dim].second ); }

    /** functor: test if handle is not of type */
  struct not_type_test {
      inline not_type_test( EntityType type ) : mType(type) {}
      inline bool operator()( EntityHandle handle )
        { return TYPE_FROM_HANDLE(handle) != mType; }
      EntityType mType;
  };

    /** functor: test if handle is of type */
  struct type_test {
      inline type_test( EntityType type ) : mType(type) {}
      inline bool operator()( EntityHandle handle )
        { return TYPE_FROM_HANDLE(handle) == mType; }
      EntityType mType;
  };

    /** functor: test if handle is not of dimension */
  struct not_dim_test {
      inline not_dim_test( int dimension ) : mDim(dimension) {}
      inline bool operator()( EntityHandle handle ) const
        { return DIM_FROM_HANDLE(handle) != mDim; }
      int mDim;
  };

    /** functor: test if handle is of dimension */
  struct dim_test {
      inline dim_test( int dimension ) : mDim(dimension) {}
      inline bool operator()( EntityHandle handle ) const
        { return DIM_FROM_HANDLE(handle) == mDim; }
      int mDim;
  };
  
    /** Iterate over range of handles.  That is, given [first_handle,last_handle],
     *  step through all contained values.
     */
  struct hdl_iter {
    EntityHandle h;
    hdl_iter( EntityHandle val ) : h(val) {}
    hdl_iter& operator++() { ++h; return *this; }
    hdl_iter& operator--() { --h; return *this; }
    hdl_iter operator++(int) { return hdl_iter(h++); }
    hdl_iter operator--(int) { return hdl_iter(h--); }
    hdl_iter& operator+=(size_t s) { h += s; return *this; }
    hdl_iter& operator-=(size_t s) { h -= s; return *this; }
    EntityHandle operator*() const { return h; }
    bool operator==(hdl_iter other) const { return h == other.h; }
    bool operator!=(hdl_iter other) const { return h != other.h; }
    bool operator< (hdl_iter other) const { return h <  other.h; }
    bool operator> (hdl_iter other) const { return h >  other.h; }
    bool operator<=(hdl_iter other) const { return h <= other.h; }
    bool operator>=(hdl_iter other) const { return h >= other.h; }
    
    struct iterator_category : public std::random_access_iterator_tag { };
    typedef EntityID difference_type;
    typedef EntityHandle value_type;
    typedef EntityHandle* pointer;
    typedef EntityHandle& reference;
  };
  
};

inline MeshSet::hdl_iter::difference_type
 operator-( const MeshSet::hdl_iter a, const MeshSet::hdl_iter b )
  { return (MeshSet::hdl_iter::difference_type)a.h 
         - (MeshSet::hdl_iter::difference_type)b.h; }


  //! create an empty meshset
MeshSet::MeshSet()
  : mFlags(0), mParentCount(ZERO), mChildCount(ZERO), mContentCount(ZERO)
{ }

  //! create an empty meshset
MeshSet::MeshSet(unsigned flg)
        : mFlags((unsigned char)flg), mParentCount(ZERO), mChildCount(ZERO), mContentCount(ZERO)
{ }

  //! destructor
MeshSet::~MeshSet() 
{
  if (mChildCount == MANY)
    free( childMeshSets.ptr[0] );
  if (mParentCount == MANY)
    free( parentMeshSets.ptr[0] );
  if (mContentCount == MANY)
    free( contentList.ptr[0] );
  mChildCount = mParentCount = mContentCount = ZERO;
}
  
ErrorCode MeshSet::set_flags( unsigned flg, EntityHandle my_handle, AEntityFactory* adjacencies ) 
{
  if(ZERO != mContentCount) {
    ErrorCode result = convert(flg, my_handle, adjacencies);
    if(MB_SUCCESS != result) return result;
  }
  mFlags = (unsigned char)flg;
  return MB_SUCCESS;
}
    

//! get all children pointed to by this meshset
const EntityHandle* MeshSet::get_children( int& count_out ) const 
{ 
  count_out = mChildCount;
  if (count_out < MANY)
    return childMeshSets.hnd;

  count_out = childMeshSets.ptr[1] - childMeshSets.ptr[0];
  return childMeshSets.ptr[0];
}

//! get all parents pointed to by this meshset
const EntityHandle* MeshSet::get_parents( int& count_out ) const
{ 
  count_out = mParentCount;
  if (count_out < MANY)
    return parentMeshSets.hnd;

  count_out = parentMeshSets.ptr[1] - parentMeshSets.ptr[0];
  return parentMeshSets.ptr[0];
}

//! return the number of children pointed to by this meshset
int MeshSet::num_children() const 
{
  if (mChildCount < MANY)
    return mChildCount;
  else
    return childMeshSets.ptr[1] - childMeshSets.ptr[0];
}

//! return the number of parents pointed to by this meshset
int MeshSet::num_parents() const
{
  if (mParentCount < MANY)
    return mParentCount;
  else
    return parentMeshSets.ptr[1] - parentMeshSets.ptr[0];
}
  
inline ErrorCode MeshSet::clear( EntityHandle myhandle, AEntityFactory* adjacencies )
{ 
  if (tracking())
    remove_adjacencies( myhandle, adjacencies );
  if (mContentCount == MANY)
    free( contentList.ptr[0] );
  mContentCount = ZERO;
  return MB_SUCCESS;
}
  
inline ErrorCode MeshSet::clear_all( EntityHandle myhandle, AEntityFactory* adjacencies )
{ 
  ErrorCode rval = clear( myhandle, adjacencies );
  if (mChildCount == MANY)
    free( childMeshSets.ptr[0] );
  mChildCount = ZERO;
  if (mParentCount == MANY)
    free( parentMeshSets.ptr[0] );
  mParentCount = ZERO;
  return rval;
}

inline const EntityHandle* MeshSet::get_contents( size_t& count_out ) const
{
  if (mContentCount == MANY) {
    count_out = contentList.ptr[1] - contentList.ptr[0];
    return contentList.ptr[0];
  }
  else {
    count_out = mContentCount;
    return contentList.hnd;
  }
}

inline EntityHandle* MeshSet::get_contents( size_t& count_out )
{
  if (mContentCount == MANY) {
    count_out = contentList.ptr[1] - contentList.ptr[0];
    return contentList.ptr[0];
  }
  else {
    count_out = mContentCount;
    return contentList.hnd;
  }
}

inline ErrorCode MeshSet::get_entities(std::vector<EntityHandle>& entities) const
{
  size_t count;
  const EntityHandle* ptr = get_contents( count );
  if (vector_based()) {
    size_t old_size = entities.size();
    entities.resize( count+old_size );
    std::copy( ptr, ptr+count, entities.begin()+old_size );
  }
  else {
    assert(count%2 == 0);
    for (size_t i = 0; i < count; i += 2) 
      std::copy( hdl_iter(ptr[i]), hdl_iter(ptr[i+1]+1), std::back_inserter(entities) );
  }
  return MB_SUCCESS;
} 

inline ErrorCode MeshSet::get_entities(Range& entities) const
{
  size_t count;
  const EntityHandle* ptr = get_contents( count );
  if (vector_based()) {
    std::copy( ptr, ptr+count, range_inserter(entities) );
  }
  else {
    assert(count%2 == 0);
    Range::iterator in = entities.begin();
    for (size_t i = 0; i < count; i += 2) 
      in = entities.insert( in, ptr[i], ptr[i+1] );
  }
  return MB_SUCCESS;
} 


  //! get all entities in this MeshSet with the specified type
inline ErrorCode MeshSet::get_entities_by_type(EntityType type,
                                               std::vector<EntityHandle> &entity_list
                                               ) const
{
  size_t count;
  const EntityHandle* ptr = get_contents( count );
  if (MBMAXTYPE == type) {
    return get_entities( entity_list);
  }
  else if (vector_based()) {
    std::remove_copy_if( ptr, ptr+count, 
                         std::back_inserter( entity_list ),
                         not_type_test(type) );
  }
  else {
    size_t idx = std::lower_bound( ptr, ptr+count, FIRST_HANDLE(type) ) - ptr;
    if (idx < count && TYPE_FROM_HANDLE(ptr[idx]) == type) {
      if (idx % 2) { // only part of first block is of type
        std::copy( hdl_iter(FIRST_HANDLE(type)), hdl_iter(ptr[idx]+1), std::back_inserter( entity_list ) );
        ++idx;
      }
      for (; idx < count; idx += 2) {
        if (TYPE_FROM_HANDLE(ptr[idx+1]) == type) // whole block is of type
          std::copy( hdl_iter(ptr[idx]), hdl_iter(ptr[idx+1]+1), std::back_inserter( entity_list ) );
        else {
          if (TYPE_FROM_HANDLE(ptr[idx]) == type) // part of last block is of type
            std::copy( hdl_iter(ptr[idx]), hdl_iter(LAST_HANDLE(type)), std::back_inserter( entity_list ) );
          break;
        }
      }
    }
  }

  return MB_SUCCESS;
}


inline ErrorCode MeshSet::get_entities_by_type( EntityType type,
                                                    Range& entity_list) const
{
  size_t count;
  const EntityHandle* ptr = get_contents( count );
  if (MBMAXTYPE == type) {
    return get_entities( entity_list );
  }
  else if (vector_based()) {
    std::remove_copy_if( ptr, ptr+count, 
                         range_inserter( entity_list ),
                         not_type_test(type) );
  }
  else {
    size_t idx = std::lower_bound( ptr, ptr+count, FIRST_HANDLE(type) ) - ptr;
    Range::iterator in = entity_list.begin();
    if (idx < count && TYPE_FROM_HANDLE(ptr[idx]) == type) {
      if (idx % 2) { // only part of first block is of type
        in = entity_list.insert( in, FIRST_HANDLE(type), ptr[idx] );
        ++idx;
      }
      for (; idx < count; idx += 2) {
        if (TYPE_FROM_HANDLE(ptr[idx+1]) == type) // whole block is of type
          in = entity_list.insert( in, ptr[idx], ptr[idx+1] );
        else {
          if (TYPE_FROM_HANDLE(ptr[idx]) == type) // part of last block is of type
            entity_list.insert( in, ptr[idx], LAST_HANDLE(type) );
          break;
        }
      }
    }
  }

  return MB_SUCCESS;
}

  //! return the number of entities with the given type contained in this meshset
inline unsigned int MeshSet::num_entities_by_type(EntityType type) const
{
  unsigned int result;
  size_t count;
  const EntityHandle* ptr = get_contents( count );
  if (MBMAXTYPE == type) {
    return num_entities( );
  }
  else if (vector_based()) {
    #ifndef __SUNPRO_CC
      result = std::count_if( ptr, ptr+count, type_test(type) );
    #else
      std::count_if( ptr, ptr+count, type_test(type), result );
    #endif
  }
  else {
    result = 0;
    size_t idx = std::lower_bound( ptr, ptr+count, FIRST_HANDLE(type) ) - ptr;
    if (idx < count && TYPE_FROM_HANDLE(ptr[idx]) == type) {
      if (idx % 2) { // only part of first block is of type
        result += ptr[idx] - FIRST_HANDLE(type) + 1;
        ++idx;
      }
      for (; idx < count; idx += 2) {
        if (TYPE_FROM_HANDLE(ptr[idx+1]) == type) // whole block is of type
          result += ptr[idx+1] - ptr[idx] + 1;
        else {
          if (TYPE_FROM_HANDLE(ptr[idx]) == type) // part of last block is of type
            result += LAST_HANDLE(type) - ptr[idx] + 1;
          break;
        }
      }
    }
  }

  return result;
}

inline ErrorCode MeshSet::get_entities_by_dimension( int dimension,
                                                         std::vector<EntityHandle> &entity_list
                                                         ) const
{
  size_t count;
  const EntityHandle* ptr = get_contents( count );
  if (vector_based()) {
    std::remove_copy_if( ptr, ptr+count, 
                         std::back_inserter( entity_list ),
                         not_dim_test(dimension) );
  }
  else {
    size_t idx = std::lower_bound( ptr, ptr+count, FIRST_OF_DIM(dimension) ) - ptr;
    if (idx < count && DIM_FROM_HANDLE(ptr[idx]) == dimension) {
      if (idx % 2) { // only part of first block is of type
        std::copy( hdl_iter(FIRST_OF_DIM(dimension)), hdl_iter(ptr[idx]+1), std::back_inserter( entity_list ) );
        ++idx;
      }
      for (; idx < count; idx += 2) {
        if (DIM_FROM_HANDLE(ptr[idx+1]) == dimension) // whole block is of type
          std::copy( hdl_iter(ptr[idx]), hdl_iter(ptr[idx+1]+1), std::back_inserter( entity_list ) );
        else {
          if (DIM_FROM_HANDLE(ptr[idx]) == dimension) // part of last block is of type
            std::copy( hdl_iter(ptr[idx]), hdl_iter(LAST_OF_DIM(dimension)), std::back_inserter( entity_list ) );
          break;
        }
      }
    }
  }

  return MB_SUCCESS;
}


inline ErrorCode MeshSet::get_entities_by_dimension( int dimension,
                                                         Range& entity_list) const
{
  size_t count;
  const EntityHandle* ptr = get_contents( count );
  if (vector_based()) {
    std::remove_copy_if( ptr, ptr+count, 
                         range_inserter( entity_list ),
                         not_dim_test(dimension) );
  }
  else {
    size_t idx = std::lower_bound( ptr, ptr+count, FIRST_OF_DIM(dimension) ) - ptr;
    Range::iterator in = entity_list.begin();
    if (idx < count && DIM_FROM_HANDLE(ptr[idx]) == dimension) {
      if (idx % 2) { // only part of first block is of type
        in = entity_list.insert( in, FIRST_OF_DIM(dimension), ptr[idx] );
        ++idx;
      }
      for (; idx < count; idx += 2) {
        if (DIM_FROM_HANDLE(ptr[idx+1]) == dimension) // whole block is of type
          in = entity_list.insert( in, ptr[idx], ptr[idx+1] );
        else {
          if (DIM_FROM_HANDLE(ptr[idx]) == dimension) // part of last block is of type
            entity_list.insert( in, ptr[idx], LAST_OF_DIM(dimension) );
          break;
        }
      }
    }
  }

  return MB_SUCCESS;
}

  //! return the number of entities with the given type contained in this meshset
inline unsigned int MeshSet::num_entities_by_dimension(int dimension) const
{
  unsigned int result;
  size_t count;
  const EntityHandle* ptr = get_contents( count );
  if (vector_based()) {
    #ifndef __SUNPRO_CC
      result = std::count_if( ptr, ptr+count, dim_test(dimension) );
    #else
      std::count_if( ptr, ptr+count, dim_test(dimension), result );
    #endif
  }
  else {
    result = 0;
    size_t idx = std::lower_bound( ptr, ptr+count, FIRST_OF_DIM(dimension) ) - ptr;
    if (idx < count && DIM_FROM_HANDLE(ptr[idx]) == dimension) {
      if (idx % 2) { // only part of first block is of type
        result += ptr[idx] - FIRST_OF_DIM(dimension) + 1;
        ++idx;
      }
      for (; idx < count; idx += 2) {
        if (DIM_FROM_HANDLE(ptr[idx+1]) == dimension) // whole block is of type
          result += ptr[idx+1] - ptr[idx] + 1;
        else {
          if (DIM_FROM_HANDLE(ptr[idx]) == dimension) // part of last block is of type
            result += LAST_OF_DIM(dimension) - ptr[idx] + 1;
          break;
        }
      }
    }
  }

  return result;
}

inline ErrorCode MeshSet::get_non_set_entities( Range& range ) const
{
  size_t count;
  const EntityHandle* ptr = get_contents( count );
  if (vector_based()) {
    std::remove_copy_if( ptr, ptr+count, 
                         range_inserter( range ),
                         type_test(MBENTITYSET) );
  }
  else {
    Range::iterator in = range.begin();
    for (size_t idx = 0; idx < count; idx += 2) {
      if (TYPE_FROM_HANDLE(ptr[idx+1]) != MBENTITYSET)
        in = range.insert( in, ptr[idx], ptr[idx+1] );
      else {
        if (TYPE_FROM_HANDLE(ptr[idx]) != MBENTITYSET)
          in = range.insert( in, ptr[idx], LAST_HANDLE( MBENTITYSET - 1 ) );
        break;
      }
    }
  }

  return MB_SUCCESS;
}

inline bool MeshSet::contains_entities(const EntityHandle *entities, 
                                         int num_ents,
                                         const int op) const
{
  size_t count;
  const EntityHandle* const ptr = get_contents( count );
  const EntityHandle* const end = ptr + count;
  size_t found_count = 0;
  if (vector_based()) {
    for (int i = 0; i < num_ents; ++i)
      if (std::find( ptr, end, entities[i] ) < end)
        ++found_count; 
  }
  else {
    assert(0 == count % 2);
    for (int i = 0; i < num_ents; ++i) {
      const unsigned long idx = std::lower_bound( ptr, end, entities[i] ) - ptr;
      if (idx < count && (idx%2 != 0 || ptr[idx] == entities[i]))
        ++found_count;
    }
  }

  return found_count >= ((Interface::INTERSECT == op) ? (unsigned)num_ents : 1u);
}



//! subtract/intersect/unite meshset_2 from/with/into meshset_1; modifies meshset_1
inline ErrorCode MeshSet::subtract(const MeshSet *meshset_2,
                                       EntityHandle my_handle,
                                       AEntityFactory* adjacencies)
{
  size_t count;
  const EntityHandle* const ptr = meshset_2->get_contents( count );
  if (meshset_2->vector_based()) 
    return remove_entity_vector( ptr, count, my_handle, adjacencies );
  else
    return remove_entity_ranges( ptr, count, my_handle, adjacencies );
}

inline ErrorCode MeshSet::unite(const MeshSet *meshset_2,
                                    EntityHandle my_handle,
                                    AEntityFactory* adjacencies)
{
  size_t count;
  const EntityHandle* const ptr = meshset_2->get_contents( count );
  if (meshset_2->vector_based()) 
    return insert_entity_vector( ptr, count, my_handle, adjacencies );
  else
    return insert_entity_ranges( ptr, count, my_handle, adjacencies );
}

//! add these entities to this meshset
inline ErrorCode MeshSet::add_entities(const EntityHandle *entity_handles,
                                           const int num_ents,
                                           EntityHandle my_handle,
                                           AEntityFactory* adjacencies)
{
  return insert_entity_vector( entity_handles, num_ents, my_handle, adjacencies );
}

  //! add these entities to this meshset
inline ErrorCode MeshSet::add_entities(const Range &entities,
                                           EntityHandle my_handle,
                                           AEntityFactory* adjacencies)
{
  return insert_entity_ranges( entities, my_handle, adjacencies );
}

  //! add these entities to this meshset
inline ErrorCode MeshSet::remove_entities(const Range& entities,
                                              EntityHandle my_handle,
                                              AEntityFactory* adjacencies)
{
  return remove_entity_ranges(  entities, my_handle, adjacencies );
}


  //! remove these entities from this meshset
inline ErrorCode MeshSet::remove_entities(const EntityHandle *entities,
                                              const int num_ents,
                                              EntityHandle my_handle,
                                              AEntityFactory* adjacencies)
{
  return remove_entity_vector( entities, num_ents, my_handle, adjacencies );
}

  //! return the number of entities contained in this meshset
unsigned int MeshSet::num_entities() const
{
  size_t count;
  const EntityHandle* list = get_contents( count );
  if (vector_based())
    return count;
  
  int result = 0;
  const EntityHandle *const end = list + count;
  for (; list < end; list += 2)
    result += list[1] - list[0] + 1;
  return result;
}
  
} // namespace moab

#endif
