element_tree.hpp

Go to the documentation of this file.

#include <vector>
#include <set>
#include <iostream>
#include <map>
#include <algorithm>
#include <bitset>
#include <numeric>
#include <cmath>
#include <tr1/unordered_map>
#include <limits>

#include "common_tree.hpp"

#ifndef ELEMENT_TREE_HPP
#define ELEMENT_TREE_HPP
namespace moab {
//forward declarations

template< typename _Entity_handles, 
      typename _Box, 
      typename _Moab,
      typename _Parametrizer> class Element_tree;

//non-exported functionality
namespace  {
namespace _element_tree { 
template< typename Iterator>
struct Iterator_comparator{
typedef typename Iterator::value_type Value;
bool operator()(const Value & a, const Value & b){
    return a->second.second.to_ulong() < b->second.second.to_ulong();
}
}; //Iterator_comparator

template< typename Data>
struct Split_comparator {
  //we minimizes ||left| - |right|| + |middle|^2 
  double split_objective( const Data & a) const {
    if (a.second.sizes[ 2]==0 || a.second.sizes[ 0] == 0){
        return std::numeric_limits< std::size_t>::max();
    }
    const double total = a.second.sizes[ 0]+a.second.sizes[ 2];
    const int max = 2*(a.second.sizes[ 2]>a.second.sizes[ 0]);
    
    return (a.second.sizes[ max] -a.second.sizes[ 2*(1-(max==2))])/total;
  }
  bool operator()( const Data & a, const Data & b) const {
    return split_objective( a) < split_objective( b);
  }
}; //Split_comparator

template< typename Partition_data, typename Box>
void correct_bounding_box( const Partition_data & data, Box & box, 
               const int child){
    const int dim = data.dim;
    switch( child){
        case 0:
            box.max[ dim] = data.left_rightline;
            break;
        case 1:
            box.max[ dim] = data.right_line;
            box.min[ dim] = data.left_line;
            break;
        case 2:
            box.min[ dim] = data.right_leftline;
            break;
    }
    #ifdef ELEMENT_TREE_DEBUG
    print_vector( data.bounding_box.max);
    print_vector( data.bounding_box.min);
    print_vector( box.max);
    print_vector( box.min);
    #endif
}

template< typename Box>
struct _Partition_data{
    typedef _Partition_data< Box> Self;
    //default constructor
    _Partition_data():sizes(3,0),dim(0){}
    _Partition_data( const Self & f){ *this=f; }
    _Partition_data( const Box & _box, int _dim): sizes(3,0),
    bounding_box( _box), split((_box.max[ _dim] + _box.min[ _dim])/2.0), 
    left_line( split), right_line( split), dim( _dim){}
    _Partition_data& operator=( const Self & f){
        sizes = f.sizes;
        bounding_box = f.bounding_box;
        split = f.split;
        left_line = f.left_line;
        right_line = f.right_line;
        right_leftline=f.right_leftline;
        left_rightline=f.left_rightline;
        dim = f.dim;
        return *this;
    }
    std::vector< std::size_t> sizes;
    Box bounding_box;
    double split;
    double left_line;
    double right_line;
    double right_leftline;
    double left_rightline;
    int dim;
    std::size_t left()  const { return sizes[ 0]; }
    std::size_t middle()const { return sizes[ 1]; }
    std::size_t right() const { return sizes[ 2]; }
}; // Partition_data

template< typename _Entity_handles, typename _Entities>
class _Node{
    //public types:
    public:
    typedef _Entity_handles Entity_handles;
    typedef _Entities Entities;

    //private types:
    private:
    typedef _Node< _Entity_handles, _Entities> Self;

    //Constructors
    public:
    //Default constructor
    _Node(): children(3,-1), 
    left_( children[ 0]), middle_( children[ 1]),
    right_( children[ 2]),
    dim( -1), split( 0),
         left_line( 0), right_line( 0),
        entities( 0) {}

    //Copy constructor
    _Node( const Self & from):
    children( from.children), 
    left_( children[ 0]), middle_( children[ 1]),
    right_( children[ 2]),
    dim( from.dim), split( from.split),
    left_line( from.left_line), right_line( from.right_line), 
    entities( from.entities) {}

    public:
    template< typename Iterator>
    void assign_entities(const Iterator & begin, const Iterator & end){
        entities.reserve( std::distance( begin, end)); 
        for( Iterator i = begin; i != end; ++i){
            entities.push_back( std::make_pair((*i)->second.first, 
                                (*i)->first));
        }
    }

    // Functionality
    public: 
    bool leaf() const { return children[ 0] == -1 && 
                       children[ 1] == -1 && 
                       children[ 2] == -1; }
    Self& operator=( const Self & from){
        children=from.children;
        dim=from.dim;
        left_ = from.left_;
        middle_ = from.middle_;
        right_ = from.right_;
        split=from.split;
        left_line=from.left_line;
        right_line=from.right_line;
        entities=from.entities;
        return *this;
    }
    template< typename Box>
    Self& operator=( const _Partition_data< Box> & from){
        dim=from.dim;
        split=from.split;
        left_line=from.left_line;
        right_line=from.right_line;
        return *this;
    }
    //private data members:
    private:
    //indices of children
    std::vector< int> children;
    int & left_, middle_, right_;
    int dim; //split dimension
    double split; //split position
    double left_line;
    double right_line;
    Entities entities;

    //Element_tree can touch my privates.
    template< Entity_handles, typename B> friend class moab::Element_tree;
}; //class Node

} //namespace _element_tree
} // anon namespace 

template< typename _Entity_handles,
      typename _Box,
      typename _Moab, 
      typename _Parametrizer>
class Element_tree {

//public types
public:
    typedef  _Entity_handles Entity_handles;
    typedef  _Box Box;
    typedef  _Moab Moab;
    typedef _Parametrizer Parametrizer;
    typedef typename Entity_handles::value_type Entity_handle;
    
//private types
private: 
    typedef Element_tree< _Entity_handles, 
                  _Box, 
                  _Moab, 
                  _Parametrizer> Self; 
    typedef std::pair< Box, Entity_handle>  Leaf_element;
    typedef _element_tree::_Node< Entity_handles, std::vector< Leaf_element> > Node;
    //int is because we only need to store  
    #define MAX_ITERATIONS 2
    typedef  common_tree::_Element_data< Box, std::bitset<NUM_DIM*MAX_ITERATIONS*2> > 
                                Element_data;
    typedef std::vector< Node> Nodes;
    //TODO: we really want an unordered map here, make sure this is kosher..
    typedef std::tr1::unordered_map< Entity_handle, Element_data> 
                                Element_map;
    typedef typename std::vector< typename Element_map::iterator> 
                                Element_list;
    typedef _element_tree::_Partition_data< Box> Partition_data;
//public methods
public:
//Constructor
Element_tree( Entity_handles & _entities, Moab & _moab, Box & _bounding_box, 
          Parametrizer & _entity_contains): 
    entity_handles_( _entities), tree_(), moab( _moab), 
    bounding_box( _bounding_box), entity_contains( _entity_contains){
    tree_.reserve( _entities.size());
    Element_map element_map( _entities.size());
    Partition_data _data;
    common_tree::construct_element_map( entity_handles_, element_map, 
                        _data.bounding_box, moab);
    bounding_box = _data.bounding_box;
    _bounding_box = bounding_box;
    Element_list element_ordering( element_map.size());
    std::size_t index = 0;
    for(typename Element_map::iterator i = element_map.begin(); 
                      i != element_map.end(); ++i, ++index){
        element_ordering[ index] = i;
    }
    //We only build nonempty trees
    if( element_ordering.size()){ 
        //initially all bits are set
        std::bitset< 3> directions( 7);
        tree_.push_back( Node());
        int depth = 0;
        build_tree( element_ordering.begin(), 
                element_ordering.end(),
                0, directions, _data, depth);
        std::cout << "depth: " << depth << std::endl; 
    }
}

//Copy constructor
Element_tree( Self & s): entity_handles_( s.entity_handles_), 
             tree_( s.tree_), moab( s.moab), 
             bounding_box( s.bounding_box){}

//private functionality
private:

template< typename Iterator, typename Split_data>
void compute_split( Iterator & begin, Iterator & end, 
             Split_data & split_data, bool iteration=false){
    typedef typename Iterator::value_type::value_type Map_value_type;
    typedef typename Map_value_type::second_type::second_type Bitset;
    //we will update the left/right line
    double & left_line = split_data.left_line;
    double & right_line = split_data.right_line;
    double & split = split_data.split;
    const int & dim = split_data.dim;
    #ifdef ELEMENT_TREE_DEBUG
    std::cout << std::endl; 
    std::cout << "-------------------" << std::endl; 
    std::cout << "dim: " << dim << " split: " << split << std::endl;
    std::cout << "bounding_box min: "; 
    print_vector( split_data.bounding_box.min); 
    std::cout << "bounding_box max: "; 
    print_vector( split_data.bounding_box.max);
    #endif
    //for each elt determine if left/middle/right
    for(Iterator i = begin; i != end; ++i){
        const Box & box =  (*i)->second.first;
        Bitset & bits =  (*i)->second.second;
        //will be 0 if on left, will be 1 if in the middle
        //and 2 if on the right;
        const bool on_left = (box.max[ dim] < split);
        const bool on_right = (box.min[ dim] > split);
        const bool in_middle = !on_left && !on_right;
        //set the corresponding bits in the bit vector
        // looks like: [x_1 = 00 | x_2 = 00 | .. | z_1 = 00 | z_2 = 00]
        // two bits, left = 00, middle = 01, right = 10
        const int index = 4*dim + 2*iteration;
        if( on_left){ split_data.sizes[ 0]++; }
        else if(in_middle){
            split_data.sizes[ 1]++;
            bits.set( index, 1);
            left_line  = std::min( left_line,  box.min[ dim]);
            right_line = std::max( right_line, box.max[ dim]);
        }else if( on_right){
            bits.set( index+1, 1);
            split_data.sizes[ 2]++;
        }
    }
    #ifdef ELEMENT_TREE_DEBUG
    std::size_t _count = std::accumulate( split_data.sizes.begin(), 
                          split_data.sizes.end(), 0);
    std::size_t total = std::distance( begin, end);
    if( total != _count ){
        std::cout << total << "vs. " <<  _count << std::endl;

    }
    std::cout <<  " left_line: " << left_line;
    std::cout <<  " right_line: " << right_line << std::endl;
    std::cout << "co/mputed partition size: ";
    print_vector( split_data.sizes);
    std::cout << "-------------------" << std::endl; 
    #endif
}

template< typename Split_data>
bool update_split_line( Split_data & data) const{
    const int max = 2*(data.sizes[ 2]>data.sizes[ 0]);
    const int min = 2*(1-(max==2));
    bool one_side_empty = data.sizes[ max]==0 || data.sizes[ min]==0;
    double balance_ratio = data.sizes[ max] - data.sizes[ min];
    //if ( !one_side_empty && balance_ratio < .05*total){ return false; }
    if( !one_side_empty){
        //if we have some imbalance on left/right 
        //try to fix the situation 
        balance_ratio /= data.sizes[ max];
        data.split += (max-1)*balance_ratio*(data.split/2.0);
    }else{
        //if the (left) side is empty move the split line just past the 
        //extent of the (left) line of the middle box.
        //if middle encompasses everything then wiggle 
        //the split line a bit and hope for the best..
        const double left_distance = 
                    std::abs(data.left_line-data.split);    
            const double right_distance = 
                    std::abs(data.right_line-data.split);
        if( (data.sizes[ 0] == 0) && (data.sizes[ 2] != 0)){
            data.split += right_distance;
        }else if (data.sizes[ 2]==0 && data.sizes[ 0] != 0){
            data.split -= left_distance;
        }else{
            data.split *=1.05;
        }
    }
    data.left_line = data.right_line = data.split;
    data.sizes.assign( data.sizes.size(), 0);
    return true;
}

template< typename Iterator, 
      typename Split_data, 
      typename Directions>
void determine_split( Iterator & begin, 
              Iterator & end, 
              Split_data & data, 
              const Directions & directions){ 
    typedef typename Iterator::value_type Pair;
    typedef typename Pair::value_type Map_value_type;
    typedef typename Map_value_type::second_type::second_type Bitset;
    typedef typename Map_value_type::second_type::first_type Box;
    typedef typename std::map< std::size_t, Split_data> Splits;
    typedef typename Splits::value_type Split;  
    typedef _element_tree::Split_comparator< Split> Comparator;
    Splits splits;
    for (std::size_t dir = 0; dir < directions.size(); ++dir){
        if( directions.test( dir)){
            Split_data split_data( data.bounding_box, dir);
            compute_split( begin, end, split_data);
            splits.insert( std::make_pair(2*dir, split_data));
            if( update_split_line( split_data)){
                compute_split( begin, end, split_data, true);
                splits.insert( std::make_pair( 2*dir+1,
                                   split_data) );
            }
        }
    }
    Split best = *std::min_element( splits.begin(), splits.end(), 
                    Comparator());
    #ifdef ELEMENT_TREE_DEBUG
    std::cout << "best: " << Comparator().split_objective( best) << " ";
    print_vector( best.second.sizes);
    #endif 
    const int dir = best.first/2;
    const int iter = best.first%2;
    double & left_rightline = best.second.left_rightline= 
                        best.second.bounding_box.min[ dir];
    double & right_leftline = best.second.right_leftline = 
                    best.second.bounding_box.max[ dir];
    Bitset mask( 0);
    mask.flip( 4*dir+2*iter).flip( 4*dir+2*iter+1);
    for(Iterator i = begin; i != end; ++i){
        Bitset & bits = (*i)->second.second;
        const Box & box =  (*i)->second.first;
        //replace 12 bits with just two.
        bits &= mask;
        bits >>= 4*dir+2*iter;
        //if box is labeled left/right but properly contained 
        //in the middle, move the element into the middle.
        //we can shrink the size of left/right
        switch( bits.to_ulong()){
            case 0:
                if ( box.max[ dir] > best.second.left_line){ 
                    left_rightline = std::max( 
                        left_rightline, box.max[ dir]);
                }
                break;
            case 2:
                if ( box.min[ dir] < best.second.right_line){ 
                    right_leftline = std::min( 
                        right_leftline, box.max[ dir]);
                }
                break;
        }
    }
    data = best.second;
}

//define here for now.
#define ELEMENTS_PER_LEAF 30
#define MAX_DEPTH 30
#define EPSILON 1e-1
template< typename Iterator, typename Node_index, 
      typename Directions, typename Partition_data>
void build_tree( Iterator begin, Iterator end, 
         const Node_index node, 
         const Directions & directions, 
         Partition_data & _data,
         int & depth, 
         const bool is_middle = false){
    std::size_t number_elements = std::distance(begin, end);
    if ( depth < MAX_DEPTH && 
        number_elements > ELEMENTS_PER_LEAF && 
        (!is_middle || directions.any())){
        determine_split( begin, end,  _data, directions);
        //count_sort( begin, end, _data);
        std::sort( begin, end,
             _element_tree::Iterator_comparator< Iterator>());
        //update the tree
        tree_[ node] = _data;
        Iterator middle_begin( begin+_data.left());
        Iterator middle_end( middle_begin+_data.middle());
        std::vector< int> depths(3, depth);
        //left subtree
        if( _data.left()>0){
            Partition_data data( _data);
            tree_.push_back( Node());
            tree_[ node].children[ 0] = tree_.size()-1;
            correct_bounding_box( _data, data.bounding_box, 0);
            Directions new_directions( directions);
            const bool axis_is_very_small = 
             (data.bounding_box.max[ _data.dim] - 
                data.bounding_box.min[ _data.dim] < EPSILON);
            new_directions.set( _data.dim, axis_is_very_small); 
            build_tree( begin, middle_begin, 
                    tree_[ node].children[ 0], 
                    new_directions, data, ++depths[ 0],
                    is_middle);
        }
        //middle subtree
        if( _data.middle()>0){
            Partition_data data( _data);
            tree_.push_back( Node());
            tree_[ node].children[ 1] = tree_.size()-1;
            correct_bounding_box( _data, data.bounding_box, 1);
            //force the middle subtree to split
            //in a different direction from this one
            Directions new_directions( directions);
            new_directions.flip( tree_[node].dim);
            bool axis_is_very_small = 
             (data.bounding_box.max[ _data.dim] - 
                data.bounding_box.min[ _data.dim] < EPSILON);
            new_directions.set( _data.dim, axis_is_very_small); 
            build_tree( middle_begin, middle_end,
                    tree_[ node].children[ 1], 
                    new_directions, data, 
                    ++depths[ 1], true);
        }
        //right subtree
        if( _data.right()>0){
            Partition_data data( _data);
            tree_.push_back( Node());
            tree_[ node].children[ 2] = tree_.size()-1;
            correct_bounding_box( _data, data.bounding_box, 2);
            Directions new_directions( directions);
            const bool axis_is_very_small = 
             (data.bounding_box.max[ _data.dim] - 
                data.bounding_box.min[ _data.dim] < EPSILON);
            new_directions.set( _data.dim, axis_is_very_small); 

            build_tree( middle_end, end, tree_[ node].children[ 2],
                    directions, data, ++depths[ 2], is_middle);
        }
        depth = *std::max_element(depths.begin(), depths.end());
    }
    if( tree_[ node].leaf()){
        common_tree::assign_entities(tree_[ node].entities, begin, end);
    }
}

template< typename Vector, typename Node_index, typename Result>
Result& _find_point( const Vector & point, 
                 const Node_index & index,
             Result & result) const{
    typedef typename Node::Entities::const_iterator Entity_iterator;
    typedef typename std::pair< bool, Vector> Return_type;
    const Node & node = tree_[ index];
    if( node.leaf()){
        //check each node
        for( Entity_iterator i = node.entities.begin(); 
                     i != node.entities.end(); ++i){
            if( common_tree::box_contains_point( i->first, point)){
                Return_type r = entity_contains( moab, 
                                     i->second, 
                                 point);
                if( r.first){ result = 
                    std::make_pair( i->second, r.second);
                }
                return result;
            }
        }
        return Result(0, point);
    }
    if( point[ node.dim] < node.left_line){
        return _find_point( point, node.left_, result);
    }else if( point[ node.dim] > node.right_line){
        return _find_point( point, node.right_, result);
    } else {
        Entity_handle middle =  _find_point( point, 
                            node.middle_, result);
        if( middle != 0){ return result; }
        if( point[ node.dim] < node.split){ 
            return _find_point( point, node.left_, result); 
        }
        return  _find_point( point, node.right_, result);
    }
}

//public functionality
public:
template< typename Vector, typename Result>
Result& find( const Vector & point, Result & result) const{
    typedef typename Vector::const_iterator Point_iterator;
    typedef typename Box::Pair Pair; 
    typedef typename Pair::first_type Box_iterator;
    return  _find_point( point, 0, result);
}
    

//public accessor methods
public:

//private data members  
private:
    const Entity_handles & entity_handles_;
    Nodes tree_;
    Moab & moab;
    Box bounding_box;
    Parametrizer entity_contains;

}; //class Element_tree

} // namespace moab

#endif //ELEMENT_TREE_HPP