tree_searching_perf.cpp

Go to the documentation of this file.

#include "moab/Core.hpp"
#include "moab/SpatialLocator.hpp"
#include "moab/Tree.hpp"
#include "moab/HomXform.hpp"
#include "moab/ScdInterface.hpp"
#include "moab/CartVect.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "moab/BVHTree.hpp"
#include "moab/ProgOptions.hpp"
#include "moab/CpuTimer.hpp"

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

#include <cstdlib>
#include <sstream>

using namespace moab;

ErrorCode test_locator( SpatialLocator& sl, int npoints, double& cpu_time, double& percent_outside );
ErrorCode create_hex_mesh( Interface& mb, Range& elems, int n, int dim );

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 npoints = 100, dim = 3;
    int dints = 1, dleafs = 1, ddeps = 1, csints = 0;

    ProgOptions po;
    po.addOpt< int >( "candidateplaneset,c",
                      "Candidate plane set (0=SUBDIVISION,1=SUBDIV_SNAP,2=VERTEX_MEDIAN,3=VERTEX_SAMPLE", &csints );
    po.addOpt< int >( "ints,i", "Number of doublings of intervals on each side of scd mesh", &dints );
    po.addOpt< int >( "leaf,l", "Number of doublings of maximum number of elements per leaf", &dleafs );
    po.addOpt< int >( "max_depth,m", "Number of 5-intervals on maximum depth of tree", &ddeps );
    po.addOpt< int >( "npoints,n", "Number of query points", &npoints );
    po.addOpt< int >( "dim,d", "Dimension of the mesh", &dim );
    //  po.addOpt<void>( "print,p", "Print tree details", &print_tree);
    po.parseCommandLine( argc, argv );

    std::vector< int > ints, deps, leafs;
    ints.push_back( 10 );
    for( int i = 1; i < dints; i++ )
        ints.push_back( 2 * ints[i - 1] );
    deps.push_back( 30 );
    for( int i = 1; i < ddeps; i++ )
        deps.push_back( deps[i - 1] - 5 );
    leafs.push_back( 6 );
    for( int i = 1; i < dleafs; i++ )
        leafs.push_back( 2 * leafs[i - 1] );

    ErrorCode rval = MB_SUCCESS;
    std::cout << "Tree_type"
              << " "
              << "Elems_per_leaf"
              << " "
              << "Tree_depth"
              << " "
              << "Ints_per_side"
              << " "
              << "N_elements"
              << " "
              << "search_time"
              << " "
              << "perc_outside"
              << " "
              << "initTime"
              << " "
              << "nodesVisited"
              << " "
              << "leavesVisited"
              << " "
              << "numTraversals"
              << " "
              << "leafObjectTests" << std::endl;

    // outermost iteration: # elements
    for( std::vector< int >::iterator int_it = ints.begin(); int_it != ints.end(); ++int_it )
    {
        Core mb;
        Range elems;
        rval = create_hex_mesh( mb, elems, *int_it, dim );
        if( MB_SUCCESS != rval ) return rval;

        // iteration: tree depth
        for( std::vector< int >::iterator dep_it = deps.begin(); dep_it != deps.end(); ++dep_it )
        {

            // iteration: tree max elems/leaf
            for( std::vector< int >::iterator leafs_it = leafs.begin(); leafs_it != leafs.end(); ++leafs_it )
            {

                // iteration: tree type
                for( int tree_tp = 0; tree_tp < 2; tree_tp++ )
                {
                    // create tree
                    Tree* tree;
                    if( 0 == tree_tp )
                        tree = new BVHTree( &mb );
                    else
                        tree = new AdaptiveKDTree( &mb );

                    std::ostringstream opts;
                    opts << "MAX_DEPTH=" << *dep_it << ";MAX_PER_LEAF=" << *leafs_it;
                    if( csints )
                    {
                        if( opts.str().length() > 0 ) opts << ";";
                        opts << "PLANE_SET=" << csints;
                    }
                    FileOptions fo( opts.str().c_str() );
                    rval = tree->parse_options( fo );
                    if( MB_SUCCESS != rval ) return rval;
                    SpatialLocator sl( &mb, elems, tree );

                    // call evaluation
                    double cpu_time, perc_outside;
                    rval = test_locator( sl, npoints, cpu_time, perc_outside );
                    if( MB_SUCCESS != rval ) return rval;

                    std::cout << ( tree_tp == 0 ? "BVH" : "KD" ) << " " << *leafs_it << " " << *dep_it << " " << *int_it
                              << " " << ( *int_it ) * ( *int_it ) * ( *int_it ) << " " << cpu_time << " "
                              << perc_outside << " ";

                    tree->tree_stats().output_all_stats();

                }  // tree_tp

            }  // max elems/leaf

        }  // max depth

    }  // # elements

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

    return 0;
}

ErrorCode test_locator( SpatialLocator& sl, int npoints, double& cpu_time, double& percent_outside )
{
    BoundBox box     = sl.local_box();
    CartVect box_del = box.bMax - box.bMin;

    std::vector< CartVect > test_pts( npoints ), test_res( npoints );
    std::vector< EntityHandle > ents( npoints );
    int* is_in = new int[npoints];

    double denom = 1.0 / (double)RAND_MAX;
    for( int i = 0; i < npoints; i++ )
    {
        // generate a small number of random point to test
        double rx = (double)rand() * denom, ry = (double)rand() * denom, rz = (double)rand() * denom;
        test_pts[i] = box.bMin + CartVect( rx * box_del[0], ry * box_del[1], rz * box_del[2] );
    }

    CpuTimer ct;

    // call spatial locator to locate points
    ErrorCode rval = sl.locate_points( test_pts[0].array(), npoints, &ents[0], test_res[0].array(), &is_in[0] );
    if( MB_SUCCESS != rval )
    {
        delete[] is_in;
        return rval;
    }

    cpu_time = ct.time_elapsed();

    int num_out     = std::count( is_in, is_in + npoints, false );
    percent_outside = ( (double)num_out ) / npoints;
    delete[] is_in;

    return rval;
}

ErrorCode create_hex_mesh( Interface& mb, Range& elems, int n, int dim )
{
    ScdInterface* scdi;
    ErrorCode rval = mb.query_interface( scdi );
    if( MB_SUCCESS != rval ) return rval;

    HomCoord high( n - 1, -1, -1 );
    if( dim > 1 ) high[1] = n - 1;
    if( dim > 2 ) high[2] = n - 1;
    ScdBox* new_box;
    rval = scdi->construct_box( HomCoord( 0, 0, 0 ), high, NULL, 0, new_box );
    if( MB_SUCCESS != rval ) return rval;
    rval = mb.release_interface( scdi );
    if( MB_SUCCESS != rval ) return rval;

    rval = mb.get_entities_by_dimension( 0, dim, elems );
    if( MB_SUCCESS != rval ) return rval;

    return rval;
}