kd_tree_test.cpp

Go to the documentation of this file.

#include "moab/Core.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "moab/Range.hpp"
#include "moab/CartVect.hpp"

#ifdef MOAB_HAVE_MPI
#include "moab_mpi.h"
#endif

#include <cmath>
#include <cassert>
#include <cfloat>
#include <cstdio>

#include "TestUtil.hpp"

using namespace moab;

const unsigned INTERVALS = 4;
const unsigned DEPTH     = 7;  // 3*log2(INTERVALS)+1
const char* TAG_NAME     = "TEST_DATA";

EntityHandle create_tree( AdaptiveKDTree& tool, unsigned depth, int intervals, Tag* tag_handle = 0 );
void validate_tree( AdaptiveKDTree& tool, EntityHandle root, int depth, double intervals );

void test_tree_create();
void test_leaf_merge();
void test_tree_readwrite();
void test_tree_delete();
void test_iterator_back();
void test_point_search();

int main( int argc, char** argv )
{
#ifdef MOAB_HAVE_MPI
    int fail = MPI_Init( &argc, &argv );
    if( fail ) return fail;
#else
    // silence the warning of parameters not used, in serial; there should be a smarter way :(
    argv[0] = argv[argc - argc];
#endif

    int err = RUN_TEST( test_tree_create );
    if( err )  // can't run other tests if can't create tree
        return 1;
    err += RUN_TEST( test_leaf_merge );
#ifdef MOAB_HAVE_HDF5
    err += RUN_TEST( test_tree_readwrite );
#endif
    err += RUN_TEST( test_tree_delete );
    err += RUN_TEST( test_iterator_back );
    err += RUN_TEST( test_point_search );

#ifdef MOAB_HAVE_MPI
    fail = MPI_Finalize();
    if( fail ) return fail;
#endif

    return err;
}

EntityHandle create_tree( AdaptiveKDTree& tool, unsigned depth, int intervals, Tag* tag_handle )
{
    // Create tree root
    ErrorCode err;
    EntityHandle root, leaf;
    err = tool.create_root( CartVect( 0.0 ).array(), CartVect( intervals ).array(), root );
    assert( !err );

    // Use iterator to create tree to fixed depth of DEPTH
    AdaptiveKDTreeIter iter;
    err = tool.get_tree_iterator( root, iter );CHECK_ERR( err );
    while( err == MB_SUCCESS )
    {
        if( iter.depth() < depth )
        {
            // bisect leaves along alternating axes
            AdaptiveKDTree::Plane split;
            split.norm  = iter.depth() % 3;  // alternate split axes;
            split.coord = 0.5 * ( iter.box_min()[split.norm] + iter.box_max()[split.norm] );
            err         = tool.split_leaf( iter, split );  // advances iter to first new leaf
            CHECK_ERR( err );
        }
        // if current leaf is at desired depth, advance to next one
        else
        {
            err = iter.step();
        }
    }
    CHECK( MB_ENTITY_NOT_FOUND == err );

    if( !tag_handle ) return root;

    // define a tag to use to store integer values on tree leaves
    err = tool.moab()->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, *tag_handle, MB_TAG_DENSE | MB_TAG_EXCL );CHECK_ERR( err );

    // iterate over tree setting data
    int counter = 0;
    for( err = tool.get_tree_iterator( root, iter ); !err; err = iter.step() )
    {
        // store integer value on leaf
        ++counter;
        leaf = iter.handle();
        err  = tool.moab()->tag_set_data( *tag_handle, &leaf, 1, &counter );CHECK_ERR( err );
    }

    return root;
}

void validate_tree( AdaptiveKDTree& tool, EntityHandle root, unsigned depth, int intervals, Tag data )
{
    ErrorCode err;
    const double VOL = 1.0;  // all leaves should be 1x1x1 boxes
    int val;

    // iterate over tree, verifying leaves
    AdaptiveKDTreeIter iter;
    int counter = 0;
    for( err = tool.get_tree_iterator( root, iter ); !err; err = iter.step() )
    {
        // store integer value on leaf
        ++counter;
        EntityHandle leaf = iter.handle();
        CHECK( leaf != 0 );
        CHECK_EQUAL( MBENTITYSET, TYPE_FROM_HANDLE( leaf ) );

        // check size of leaf
        const double* min = iter.box_min();
        const double* max = iter.box_max();
        double dims[]     = { max[0] - min[0], max[1] - min[1], max[2] - min[2] };
        double volume     = dims[0] * dims[1] * dims[2];
        CHECK_REAL_EQUAL( VOL, volume, DBL_EPSILON );

        // check depth of leaf
        CHECK_EQUAL( depth, iter.depth() );

        // check tag value on leaf
        err = tool.moab()->tag_get_data( data, &leaf, 1, &val );CHECK_ERR( err );
        CHECK_EQUAL( counter, val );
    }
    // check number of leaves
    const int num_leaves = intervals * intervals * intervals;
    CHECK_EQUAL( num_leaves, counter );
}

void test_tree_create()
{
    Tag tag;
    Core mb;
    AdaptiveKDTree tool( &mb );
    const EntityHandle root = create_tree( tool, DEPTH, INTERVALS, &tag );
    validate_tree( tool, root, DEPTH, INTERVALS, tag );
}

void test_leaf_merge()
{
    ErrorCode err;
    Core mb;
    AdaptiveKDTree tool( &mb );
    Tag data;
    const EntityHandle root = create_tree( tool, DEPTH, INTERVALS, &data );

    // reduce tree depth to DEPTH-1 by merging adjacent leaf pairs,
    // make new "leaf" have smaller of two data values on original pair
    AdaptiveKDTreeIter iter;
    for( err = tool.get_tree_iterator( root, iter ); !err; err = iter.step() )
    {
        // get data for first leaf
        int data1;
        EntityHandle leaf = iter.handle();
        err               = mb.tag_get_data( data, &leaf, 1, &data1 );CHECK_ERR( err );
        // tree traversal is always such that two leaves with same parent are consective
        err = iter.step();CHECK_ERR( err );
        // get data for sibling
        int data2;
        leaf = iter.handle();
        err  = mb.tag_get_data( data, &leaf, 1, &data2 );CHECK_ERR( err );
        // as we stored increasing values, these had better be increasing
        CHECK_EQUAL( 1, data2 - data1 );
        // merge leaf pair (iter can be at either one)
        err = tool.merge_leaf( iter );  // changes iter to be new "merged" leaf
        CHECK_ERR( err );
        // store smaller of two values on new leaf
        leaf = iter.handle();
        err  = mb.tag_set_data( data, &leaf, 1, &data1 );CHECK_ERR( err );
    }

    // Iterate over tree, verifying leaves and checking data
    // Initial leaves had volume of 1 : merged pairs of leaves so volume should now be 2.
    // Initial leaves were enumerated in order : merged pairs so new leaves should
    //   have data incrementing in steps of 2.
    int counter = 1;
    for( err = tool.get_tree_iterator( root, iter ); !err; err = iter.step() )
    {
        // store integer value on leaf
        int data1;
        EntityHandle leaf = iter.handle();
        err               = mb.tag_get_data( data, &leaf, 1, &data1 );CHECK_ERR( err );
        CHECK_EQUAL( counter, data1 );
        counter += 2;

        // check size of leaf
        const double* min = iter.box_min();
        const double* max = iter.box_max();
        double dims[]     = { max[0] - min[0], max[1] - min[1], max[2] - min[2] };
        double volume     = dims[0] * dims[1] * dims[2];
        CHECK_REAL_EQUAL( 2.0, volume, DBL_EPSILON );

        // check depth of leaf
        CHECK_EQUAL( DEPTH - 1, iter.depth() );
    }
}

void test_tree_readwrite()
{
    ErrorCode err;
    Tag tag;
    Core mb;
    AdaptiveKDTree tool( &mb );
    create_tree( tool, DEPTH, INTERVALS, &tag );

    // write to file
    err = mb.write_file( "tree.h5m" );CHECK_ERR( err );

    // clear everything
    mb.delete_mesh();

    // read tree from file
    err = mb.load_file( "tree.h5m" );
    remove( "tree.h5m" );CHECK_ERR( err );

    // get tag handle by name, because the handle may have changed
    err = mb.tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag );CHECK_ERR( err );

    // get root handle for tree
    Range range;
    err = tool.find_all_trees( range );
    assert( !err );
    assert( range.size() == 1 );
    EntityHandle root = range.front();  // first (only) handle

    validate_tree( tool, root, DEPTH, INTERVALS, tag );
}

void test_tree_delete()
{
    ErrorCode err;
    Core mb;
    AdaptiveKDTree tool( &mb );
    Tag data;
    create_tree( tool, DEPTH, INTERVALS, &data );

    err = tool.reset_tree();CHECK_ERR( err );

    Range ents;
    err = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &data, 0, 1, ents );CHECK_ERR( err );
    CHECK( ents.empty() );
}

void test_iterator_back()
{
    Core mb;
    AdaptiveKDTree tool( &mb );
    const EntityHandle root = create_tree( tool, DEPTH, INTERVALS );

    AdaptiveKDTreeIter iter;
    ErrorCode rval = tool.get_tree_iterator( root, iter );CHECK_ERR( rval );

    CartVect min( iter.box_min() );
    CartVect max( iter.box_max() );
    EntityHandle leaf = iter.handle();

    // going back from first location should fail.
    rval = iter.back();
    CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
    rval = tool.get_tree_iterator( root, iter );CHECK_ERR( rval );

    // make sure iterator is valid
    CHECK_REAL_EQUAL( min[0], iter.box_min()[0], DBL_EPSILON );
    CHECK_REAL_EQUAL( min[1], iter.box_min()[1], DBL_EPSILON );
    CHECK_REAL_EQUAL( min[2], iter.box_min()[2], DBL_EPSILON );
    CHECK_REAL_EQUAL( max[0], iter.box_max()[0], DBL_EPSILON );
    CHECK_REAL_EQUAL( max[1], iter.box_max()[1], DBL_EPSILON );
    CHECK_REAL_EQUAL( max[2], iter.box_max()[2], DBL_EPSILON );
    CHECK_EQUAL( leaf, iter.handle() );

    while( MB_SUCCESS == iter.step() )
    {
        // Get values at current iterator location
        CartVect next_min( iter.box_min() );
        CartVect next_max( iter.box_max() );
        EntityHandle next_leaf = iter.handle();

        // step back to previous location
        rval = iter.back();CHECK_ERR( rval );

        // check expected values for previous location
        CHECK_REAL_EQUAL( min[0], iter.box_min()[0], DBL_EPSILON );
        CHECK_REAL_EQUAL( min[1], iter.box_min()[1], DBL_EPSILON );
        CHECK_REAL_EQUAL( min[2], iter.box_min()[2], DBL_EPSILON );
        CHECK_REAL_EQUAL( max[0], iter.box_max()[0], DBL_EPSILON );
        CHECK_REAL_EQUAL( max[1], iter.box_max()[1], DBL_EPSILON );
        CHECK_REAL_EQUAL( max[2], iter.box_max()[2], DBL_EPSILON );
        CHECK_EQUAL( leaf, iter.handle() );

        // advance iterator to 'current' location
        rval = iter.step();CHECK_ERR( rval );

        // check that iterator values are correct
        CHECK_REAL_EQUAL( next_min[0], iter.box_min()[0], DBL_EPSILON );
        CHECK_REAL_EQUAL( next_min[1], iter.box_min()[1], DBL_EPSILON );
        CHECK_REAL_EQUAL( next_min[2], iter.box_min()[2], DBL_EPSILON );
        CHECK_REAL_EQUAL( next_max[0], iter.box_max()[0], DBL_EPSILON );
        CHECK_REAL_EQUAL( next_max[1], iter.box_max()[1], DBL_EPSILON );
        CHECK_REAL_EQUAL( next_max[2], iter.box_max()[2], DBL_EPSILON );
        CHECK_EQUAL( next_leaf, iter.handle() );

        // store values for next iteration
        min  = next_min;
        max  = next_max;
        leaf = next_leaf;
    }
}

void test_point_search()
{
    Core mb;
    AdaptiveKDTree tool( &mb );
    const EntityHandle root = create_tree( tool, DEPTH, INTERVALS );

    ErrorCode rval;
    EntityHandle leaf;
    AdaptiveKDTreeIter iter, iter2;

    // points to test
    CartVect left( 0.5 );
    CartVect right( CartVect( INTERVALS ) - left );

    // compare leaf search to iterator search
    rval = tool.point_search( left.array(), leaf );CHECK_ERR( rval );
    rval = tool.point_search( left.array(), iter );CHECK_ERR( rval );
    CHECK_EQUAL( leaf, iter.handle() );

    // iterator should be at 'first' leaf
    rval = tool.get_tree_iterator( root, iter2 );CHECK_ERR( rval );
    for( ;; )
    {
        CHECK_EQUAL( iter.handle(), iter2.handle() );
        CHECK_EQUAL( iter.depth(), iter2.depth() );
        CHECK_REAL_EQUAL( iter.box_min()[0], iter2.box_min()[0], DBL_EPSILON );
        CHECK_REAL_EQUAL( iter.box_min()[1], iter2.box_min()[1], DBL_EPSILON );
        CHECK_REAL_EQUAL( iter.box_min()[2], iter2.box_min()[2], DBL_EPSILON );
        CHECK_REAL_EQUAL( iter.box_max()[0], iter2.box_max()[0], DBL_EPSILON );
        CHECK_REAL_EQUAL( iter.box_max()[1], iter2.box_max()[1], DBL_EPSILON );
        CHECK_REAL_EQUAL( iter.box_max()[2], iter2.box_max()[2], DBL_EPSILON );

        rval = iter2.step();
        if( MB_ENTITY_NOT_FOUND == rval ) break;CHECK_ERR( rval );
        rval = iter.step();CHECK_ERR( rval );
    }

    // compare leaf search to iterator search
    rval = tool.point_search( right.array(), leaf, 0.0, 0.0, NULL, const_cast< EntityHandle* >( &root ) );CHECK_ERR( rval );
    rval = tool.point_search( right.array(), iter, 0.0, 0.0, NULL, const_cast< EntityHandle* >( &root ) );CHECK_ERR( rval );
    assert( iter.handle() == leaf );

    // iterator should be at 'last' leaf
    rval = tool.get_last_iterator( root, iter2 );CHECK_ERR( rval );
    for( ;; )
    {
        CHECK_EQUAL( iter.handle(), iter2.handle() );
        CHECK_EQUAL( iter.depth(), iter2.depth() );
        CHECK_REAL_EQUAL( iter.box_min()[0], iter2.box_min()[0], DBL_EPSILON );
        CHECK_REAL_EQUAL( iter.box_min()[1], iter2.box_min()[1], DBL_EPSILON );
        CHECK_REAL_EQUAL( iter.box_min()[2], iter2.box_min()[2], DBL_EPSILON );
        CHECK_REAL_EQUAL( iter.box_max()[0], iter2.box_max()[0], DBL_EPSILON );
        CHECK_REAL_EQUAL( iter.box_max()[1], iter2.box_max()[1], DBL_EPSILON );
        CHECK_REAL_EQUAL( iter.box_max()[2], iter2.box_max()[2], DBL_EPSILON );

        rval = iter2.back();
        if( MB_ENTITY_NOT_FOUND == rval ) break;CHECK_ERR( rval );
        rval = iter.back();CHECK_ERR( rval );
    }
}