test_isolation.hpp

/*
 * test_isolation.hpp
 *
 *  Created on: Mar 11, 2016
 *      Author: i-bird
 */
 
#ifndef SRC_TEST_ISOLATION_HPP_
#define SRC_TEST_ISOLATION_HPP_
 
/*
 * Test that are not test but more are isolated here
 *
 *
 */
 
BOOST_AUTO_TEST_SUITE( antoniol_test_isolation )
 
 
BOOST_AUTO_TEST_CASE( CartDecomposition_test_2D )
{
    //size_t balance_values_4p_64[] = {2.86,2.86,2.86,6.7,7.43,4.9,8.07,1.82,1.82,4.47,5.3};
 
    // Vcluster
    Vcluster & vcl = *global_v_cluster;
 
    // non-periodic boundary condition
    size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
 
    // Initialize the global VCluster
    init_global_v_cluster(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);
 
    CartDecomposition<2, float> dec(vcl);
 
    // Init DLB tool
    DLB dlb(vcl);
 
    // Physical domain
    Box<2, float> box( { 0.0, 0.0 }, { 10.0, 10.0 });
    size_t div[2];
 
    // Get the number of processor and calculate the number of sub-domain
    // for each processor (SUB_UNIT_FACTOR=64)
    size_t n_proc = vcl.getProcessingUnits();
    size_t n_sub = n_proc * SUB_UNIT_FACTOR;
 
    // Set the number of sub-domains on each dimension (in a scalable way)
    for (int i = 0; i < 2; i++)
    {
        div[i] = openfpm::math::round_big_2(pow(n_sub,1.0/2));
    }
 
    // Define ghost
    Ghost<2, float> g(0.01);
 
    // Decompose
    dec.setParameters(div, box, bc, g);
 
    // Set unbalance threshold
    dlb.setHeurisitc(DLB::Heuristic::UNBALANCE_THRLD);
    dlb.setThresholdLevel(DLB::ThresholdLevel::THRLD_MEDIUM);
 
    // Add weights to points
 
    // First create the center of the weights distribution, check it is coherent to the size of the domain
    Point<2, float> center( { 2.0, 2.0 });
 
    // Radius of the weights distribution
    float radius = 2.0;
 
    // Weight if the distribution (high)
    size_t weight_h = 5, weight_l = 1;
 
    setComputationCosts(dec, dec.getNSubSubDomains(), center, radius, weight_h, weight_l);
 
    dec.getDistribution().write("DLB_test_graph_0.vtk");
 
    dec.decompose();
 
    dec.getDistribution().write("DLB_test_graph_1.vtk");
 
    float stime = 0.0, etime = 10.0, tstep = 0.1;
 
    for(float t = stime, i = 1; t < etime; t = t + tstep, i++)
    {
 
        if(t < etime/2)
        {
            center.get(0) += tstep;
            center.get(1) += tstep;
        }
        else
        {
            center.get(0) -= tstep;
            center.get(1) -= tstep;
        }
 
        setComputationCosts(dec, dec.getNSubSubDomains(), center, radius, weight_h, weight_l);
 
        dlb.endIteration();
 
        if(dec.rebalance(dlb))
            dec.write("DLB_test_graph_cart_" + std::to_string(i+1) + "_");
 
        std::stringstream str;
        str << "DLB_test_graph_" << i + 1 << ".vtk";
        dec.getDistribution().write(str.str());
    }
 
    // For each calculated ghost box
    for (size_t i = 0; i < dec.getNIGhostBox(); i++)
    {
        Box<2,float> b = dec.getIGhostBox(i);
        size_t proc = dec.getIGhostBoxProcessor(i);
 
        // sample one point inside the box
        Point<2,float> p = b.rnd();
 
        // Check that ghost_processorsID return that processor number
        const openfpm::vector<size_t> & pr = dec.template ghost_processorID<CartDecomposition<2,float>::processor_id>(p);
 
        bool found = false;
 
        for (size_t j = 0; j < pr.size(); j++)
        {
            if (pr.get(j) == proc)
            {   found = true; break;}
        }
 
        if (found == false)
        {
            const openfpm::vector<size_t> pr2 = dec.template ghost_processorID<CartDecomposition<2,float>::processor_id>(p);
        }
 
        BOOST_REQUIRE_EQUAL(found,true);
    }
 
    // Check the consistency
 
    bool val = dec.check_consistency();
    BOOST_REQUIRE_EQUAL(val,true);
}
 
BOOST_AUTO_TEST_CASE( CartDecomposition_test_2D_sar)
{
    // Vcluster
    Vcluster & vcl = *global_v_cluster;
 
    // non-periodic boundary condition
    size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
 
    // Initialize the global VCluster
    init_global_v_cluster(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);
 
    CartDecomposition<2, float> dec(vcl);
 
    // Init DLB tool
    DLB dlb(vcl);
 
    // Physical domain
    Box<2, float> box( { 0.0, 0.0 }, { 10.0, 10.0 });
    size_t div[2];
 
    // Get the number of processor and calculate the number of sub-domain
    // for each processor (SUB_UNIT_FACTOR=64)
    size_t n_proc = vcl.getProcessingUnits();
    size_t n_sub = n_proc * SUB_UNIT_FACTOR;
 
    // Set the number of sub-domains on each dimension (in a scalable way)
    for (int i = 0; i < 2; i++)
    {
        div[i] = openfpm::math::round_big_2(pow(n_sub,1.0/2));
    }
 
    // Define ghost
    Ghost<2, float> g(0.01);
 
    // Decompose
    dec.setParameters(div, box, bc, g);
 
    // Set type of heuristic
    dlb.setHeurisitc(DLB::Heuristic::SAR_HEURISTIC);
 
    // Add weights to points
 
    // First create the center of the weights distribution, check it is coherent to the size of the domain
    Point<2, float> center( { 2.0, 2.0 });
 
    // Radius of the weights distribution
    float radius = 2.0;
 
    // Weight if the distribution (high)
    size_t weight_h = 5, weight_l = 1;
 
    size_t n_v = pow(div[0], 2);
 
    setComputationCosts(dec, n_v, center, radius, weight_h, weight_l);
 
    dec.decompose();
 
    dec.getDistribution().write("DLB_test_graph_0.vtk");
 
    float stime = 0.0, etime = 10.0, tstep = 0.1;
 
    dlb.setSimulationStartTime(0);
    dlb.setSimulationEndTime(10);
 
    for(float t = stime, i = 1; t < etime; t = t + tstep, i++)
    {
        dlb.startIteration();
 
        if(t < etime/2)
        {
            center.get(0) += tstep;
            center.get(1) += tstep;
        }
        else
        {
            center.get(0) -= tstep;
            center.get(1) -= tstep;
        }
 
        setComputationCosts(dec, n_v, center, radius, weight_h, weight_l);
 
        sleep((n_v/dec.getProcessorLoad())/vcl.getProcessingUnits());
 
        dlb.endIteration();
 
        if(dec.rebalance(dlb))
        {
            dec.write("DLB_test_graph_cart_" + std::to_string(i) + "_");
        }
 
        std::stringstream str;
        str << "DLB_test_graph_" << i << ".vtk";
        dec.getDistribution().write(str.str());
    }
 
    // For each calculated ghost box
    for (size_t i = 0; i < dec.getNIGhostBox(); i++)
    {
        Box<2,float> b = dec.getIGhostBox(i);
        size_t proc = dec.getIGhostBoxProcessor(i);
 
        // sample one point inside the box
        Point<2,float> p = b.rnd();
 
        // Check that ghost_processorsID return that processor number
        const openfpm::vector<size_t> & pr = dec.template ghost_processorID<CartDecomposition<2,float>::processor_id>(p);
 
        bool found = false;
 
        for (size_t j = 0; j < pr.size(); j++)
        {
            if (pr.get(j) == proc)
            {   found = true; break;}
        }
 
        if (found == false)
        {
            const openfpm::vector<size_t> pr2 = dec.template ghost_processorID<CartDecomposition<2,float>::processor_id>(p);
        }
 
        BOOST_REQUIRE_EQUAL(found,true);
    }
 
    // Check the consistency
 
    bool val = dec.check_consistency();
    BOOST_REQUIRE_EQUAL(val,true);
}
 
BOOST_AUTO_TEST_CASE( CartDecomposition_test_3D)
{
    // Vcluster
    Vcluster & vcl = *global_v_cluster;
 
    // non-periodic boundary condition
    size_t bc[3] = { NON_PERIODIC, NON_PERIODIC, NON_PERIODIC };
 
    // Initialize the global VCluster
    init_global_v_cluster(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);
 
    CartDecomposition<3, float> dec(vcl);
 
    // Init DLB tool
    DLB dlb(vcl);
 
    // Physical domain
    Box<3, float> box( { 0.0, 0.0, 0.0 }, { 10.0, 10.0, 10.0 });
    size_t div[3];
 
    // Get the number of processor and calculate the number of sub-domain
    // for each processor (SUB_UNIT_FACTOR=64)
    size_t n_proc = vcl.getProcessingUnits();
    size_t n_sub = n_proc * SUB_UNIT_FACTOR;
 
    // Set the number of sub-domains on each dimension (in a scalable way)
    for (int i = 0; i < 3; i++)
    {
        div[i] = openfpm::math::round_big_2(pow(n_sub,1.0/3));
    }
 
    // Define ghost
    Ghost<3, float> g(0.01);
 
    // Decompose
    dec.setParameters(div, box, bc, g);
 
    // Set unbalance threshold
    dlb.setHeurisitc(DLB::Heuristic::UNBALANCE_THRLD);
    dlb.setThresholdLevel(DLB::ThresholdLevel::THRLD_MEDIUM);
 
    // Add weights to points
 
    // First create the center of the weights distribution, check it is coherent to the size of the domain
    Point<3, float> center( { 2.0, 2.0, 2.0 });
 
    // Radius of the weights distribution
    float radius = 2.0;
 
    // Weight if the distribution (high)
    size_t weight_h = 5, weight_l = 1;
 
    size_t n_v = pow(div[0], 3);
 
    setComputationCosts3D(dec, n_v, center, radius, weight_h, weight_l);
 
    dec.decompose();
 
    dec.getDistribution().write("DLB_test_graph_0.vtk");
 
    float stime = 0.0, etime = 10.0, tstep = 0.1;
 
    for(float t = stime, i = 1; t < etime; t = t + tstep, i++)
    {
 
        if(t < etime/2)
        {
            center.get(0) += tstep;
            center.get(1) += tstep;
            center.get(2) += tstep;
        }
        else
        {
            center.get(0) -= tstep;
            center.get(1) -= tstep;
            center.get(2) -= tstep;
        }
 
        setComputationCosts3D(dec, n_v, center, radius, weight_h, weight_l);
 
        dlb.endIteration();
 
        dec.rebalance(dlb);
 
        std::stringstream str;
        str << "DLB_test_graph_" << i << ".vtk";
        dec.getDistribution().write(str.str());
    }
 
    // For each calculated ghost box
    for (size_t i = 0; i < dec.getNIGhostBox(); i++)
    {
        Box<3,float> b = dec.getIGhostBox(i);
        size_t proc = dec.getIGhostBoxProcessor(i);
 
        // sample one point inside the box
        Point<3,float> p = b.rnd();
 
        // Check that ghost_processorsID return that processor number
        const openfpm::vector<size_t> & pr = dec.template ghost_processorID<CartDecomposition<3,float>::processor_id>(p);
 
        bool found = false;
 
        for (size_t j = 0; j < pr.size(); j++)
        {
            if (pr.get(j) == proc)
            {   found = true; break;}
        }
 
        if (found == false)
        {
            const openfpm::vector<size_t> pr2 = dec.template ghost_processorID<CartDecomposition<3,float>::processor_id>(p);
        }
 
        BOOST_REQUIRE_EQUAL(found,true);
    }
 
    // Check the consistency
 
    bool val = dec.check_consistency();
    BOOST_REQUIRE_EQUAL(val,true);
}
 
 
BOOST_AUTO_TEST_CASE( dist_map_graph_use_cartesian)
{
    Vcluster & vcl = *global_v_cluster;
 
//  if(vcl.getProcessingUnits() != 3)
//  return;
 
    // non-periodic boundary condition
    size_t bc[3] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
 
    CartDecomposition<3, float> dec(vcl);
 
    // Physical domain
    Box<3, float> box( { 0.0, 0.0, 0.0 }, { 1.0, 1.0, 1.0 });
    size_t div[3] = {8,8,8};
 
    // Grid size and info
    size_t gsz[3] = {8,8,8};
    grid_sm<3,void> g_sm(gsz);
 
    // Define ghost
    Ghost<3, float> g(0.01);
 
    // Decompose
    dec.setParameters(div, box, bc, g);
    dec.decompose();
 
    grid_dist_id_iterator_dec<CartDecomposition<3,float>> it_dec(dec,gsz);
 
    size_t cnt = 0;
    while (it_dec.isNext())
    {
        cnt++;
        ++it_dec;
    }
 
    openfpm::vector<size_t> v_cnt(vcl.getProcessingUnits());
 
    // Sent and receive the size of each subgraph
    vcl.allGather(cnt, v_cnt);
    vcl.execute();
 
    cnt = 0;
    for (long int i = 0; i <= ((long int)vcl.getProcessUnitID()) - 1 ; ++i)
        cnt += v_cnt.get(i);
 
    // count the points
 
    DistGraph_CSR<aggregate<size_t[3]>, aggregate<size_t>> dg;
 
    grid_dist_id_iterator_dec<CartDecomposition<3,float>> it_dec2(dec,gsz);
    while (it_dec2.isNext())
    {
        auto key = it_dec2.get();
 
        aggregate<size_t[3]> v;
 
        v.template get<0>()[0] = key.get(0);
        v.template get<0>()[1] = key.get(1);
        v.template get<0>()[2] = key.get(2);
 
        size_t gid = g_sm.LinId(key);
        dg.add_vertex(v, gid, cnt);
 
        cnt++;
        ++it_dec2;
    }
 
    dg.init();
 
    // we ask for some random vertex
 
    std::default_random_engine rg;
    std::uniform_int_distribution<size_t> d(0,g_sm.size()-1);
 
    openfpm::vector<size_t> v_req;
 
/*  for (size_t i = 0 ; i < 16 ; i++)
    {
        size_t v = d(rg);*/
 
        if (vcl.getProcessUnitID() == 0)
        {dg.reqVertex(450);}
 
/*      dg.reqVertex(v);
    }*/
 
    dg.sync();
 
    if (vcl.getProcessUnitID() == 0)
    {
        grid_key_dx<3> key;
        // get the position information
        key.set_d(0,dg.getVertex(450).template get<0>()[0]);
        key.set_d(1,dg.getVertex(450).template get<0>()[1]);
        key.set_d(2,dg.getVertex(450).template get<0>()[2]);
 
        size_t lin_id = g_sm.LinId(key);
 
    //  BOOST_REQUIRE_EQUAL(lin_id,v_req.get(i));
 
        std::cout << "Error: " << "   " << lin_id << "    " << key.to_string() << "\n";
    }
 
/*  for (size_t i = 0 ; i < 16 ; i++)
    {
        grid_key_dx<3> key;
        // get the position information
        key.set_d(0,dg.getVertex(v_req.get(i)).template get<0>()[0]);
        key.set_d(1,dg.getVertex(v_req.get(i)).template get<0>()[1]);
        key.set_d(2,dg.getVertex(v_req.get(i)).template get<0>()[2]);
 
        size_t lin_id = g_sm.LinId(key);
 
    //  BOOST_REQUIRE_EQUAL(lin_id,v_req.get(i));
 
        std::cout << "Error: " << i << "   " << lin_id << "  " << v_req.get(i) << "\n";
    }*/
 
/*  if (vcl.getProcessUnitID() == 0)
        std::cout << "Error: " << i << "   " << lin_id << "  " << v_req.get(i) << "\n";*/
}
 
BOOST_AUTO_TEST_CASE( Parmetis_distribution_test_random_walk )
{
    typedef Point<3,float> s;
 
    Vcluster & v_cl = *global_v_cluster;
 
    // set the seed
    // create the random generator engine
    std::srand(v_cl.getProcessUnitID());
    std::default_random_engine eg;
    std::uniform_real_distribution<float> ud(0.0f, 1.0f);
 
    size_t nsz[] = { 0, 32, 4 };
    nsz[0] = 65536 * v_cl.getProcessingUnits();
 
    print_test_v( "Testing 3D random walk vector k<=",nsz[0]);
 
    // 3D test
    for (size_t i = 0; i < 3; i++ )
    {
        size_t k = nsz[i];
 
        BOOST_TEST_CHECKPOINT( "Testing 3D random walk k=" << k );
 
        Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
 
        // Grid info
        grid_sm<3, void> info( { GS_SIZE, GS_SIZE, GS_SIZE });
 
        // Boundary conditions
        size_t bc[3] = { NON_PERIODIC,NON_PERIODIC,NON_PERIODIC };
 
        // factor
        float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f);
 
        // ghost
        Ghost<3,float> ghost(0.01 / factor);
 
        // Distributed vector
        vector_dist<3,float, Point_test<float>, CartDecomposition<3, float, HeapMemory, ParMetisDistribution<3, float>>> vd(k,box,bc,ghost);
 
        auto it = vd.getIterator();
 
        while (it.isNext())
        {
            auto key = it.get();
 
            vd.template getPos<s::x>(key)[0] = ud(eg);
            vd.template getPos<s::x>(key)[1] = ud(eg);
            vd.template getPos<s::x>(key)[2] = ud(eg);
 
            ++it;
        }
 
        vd.map();
 
        vd.addComputationCosts();
 
        vd.getDecomposition().getDistribution().write("parmetis_random_walk_" + std::to_string(0) + ".vtk");
 
        // 10 step random walk
 
        for (size_t j = 0; j < 10; j++)
        {
            auto it = vd.getDomainIterator();
 
            while (it.isNext())
            {
                auto key = it.get();
 
                vd.template getPos<s::x>(key)[0] += 0.01 * ud(eg);
                vd.template getPos<s::x>(key)[1] += 0.01 * ud(eg);
                vd.template getPos<s::x>(key)[2] += 0.01 * ud(eg);
 
                ++it;
            }
 
            vd.map();
 
 
            //vd.ghost_get<>();
            //vd.getDomainIterator;
 
 
            vd.addComputationCosts();
 
            vd.getDecomposition().rebalance(10);
 
            vd.map();
 
            vd.getDecomposition().getDistribution().write("parmetis_random_walk_" + std::to_string(j+1) + ".vtk");
 
            size_t l = vd.size_local();
            v_cl.sum(l);
            v_cl.execute();
 
            // Count the local particles and check that the total number is consistent
            size_t cnt = total_n_part_lc(vd,bc);
 
            //BOOST_REQUIRE_EQUAL((size_t)k,cnt);
        }
    }
}
 
BOOST_AUTO_TEST_CASE( Parmetis_distribution_test_random_walk_2D )
{
 
    //Particle: position, type of poistion, type of animal (0 rabbit, 1 fox), dead or alive (0 or 1), time the fox stays alive without eating
    typedef Point<2,float> s;
 
    Vcluster & v_cl = *global_v_cluster;
 
    size_t JUST_EATEN = 5;
 
    // set the seed
    // create the random generator engine
    std::srand(v_cl.getProcessUnitID());
    std::default_random_engine eg;
    std::uniform_real_distribution<float> ud(0.0f, 0.5f);
 
    size_t k = 100000;
 
    print_test_v( "Testing 2D random walk vector k<=",k);
 
    BOOST_TEST_CHECKPOINT( "Testing 2D random walk k=" << k );
 
    Box<2,float> box({0.0,0.0},{1.0,1.0});
 
    // Grid info
    grid_sm<2, void> info( { GS_SIZE, GS_SIZE });
 
    // Boundary conditions
    size_t bc[2] = { PERIODIC, PERIODIC };
 
    // factor
    float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f);
 
    // ghost
    Ghost<2,float> ghost(0.01 / factor);
 
    // Distributed vector
    vector_dist<2,float, Point_test<float>, CartDecomposition<2, float, HeapMemory, ParMetisDistribution<2, float>>> vd(k,box,bc,ghost);
 
    // Init DLB tool
    DLB dlb(v_cl);
 
    // Set unbalance threshold
    dlb.setHeurisitc(DLB::Heuristic::UNBALANCE_THRLD);
    dlb.setThresholdLevel(DLB::ThresholdLevel::THRLD_MEDIUM);
 
    auto it = vd.getIterator();
 
    size_t c = 0;
    while (it.isNext())
    {
        auto key = it.get();
        if(c % 5)
        {
            vd.template getPos<s::x>(key)[0] = ud(eg);
            vd.template getPos<s::x>(key)[1] = ud(eg);
        }else{
            vd.template getPos<s::x>(key)[0] = ud(eg)*2;
            vd.template getPos<s::x>(key)[1] = ud(eg)*2;
        }
        ++it;
        ++c;
    }
 
    vd.map();
 
    vd.addComputationCosts();
 
    vd.getDecomposition().rebalance(dlb);
 
    vd.map();
 
    vd.getDecomposition().write("dec_init");
    vd.getDecomposition().getDistribution().write("parmetis_random_walk_" + std::to_string(0) + ".vtk");
    vd.write("particles_", 0, NO_GHOST);
 
    // 10 step random walk
    for (size_t j = 0; j < 50; j++)
    {
        std::cout << "Iteration " << (j+1) << "\n";
 
        auto it = vd.getDomainIterator();
 
        while (it.isNext())
        {
            auto key = it.get();
 
            vd.template getPos<s::x>(key)[0] += 0.01 * ud(eg);
            vd.template getPos<s::x>(key)[1] += 0.01 * ud(eg);
 
            ++it;
        }
 
        vd.map();
 
        vd.addComputationCosts();
 
        vd.getDecomposition().rebalance(dlb);
 
        vd.map();
 
        vd.getDecomposition().getDistribution().write("parmetis_random_walk_" + std::to_string(j+1) + ".vtk");
        vd.write("particles_", j+1, NO_GHOST);
        vd.getDecomposition().write("dec_");
 
        size_t l = vd.size_local();
        v_cl.sum(l);
        v_cl.execute();
 
        // Count the local particles and check that the total number is consistent
        //size_t cnt = total_n_part_lc(vd,bc);
 
        //BOOST_REQUIRE_EQUAL((size_t)k,cnt);
    }
 
}
 
 
BOOST_AUTO_TEST_CASE( Parmetis_distribution_test_prey_and_predators )
{
    Vcluster & v_cl = *global_v_cluster;
 
    //time the animal stays alive without eating or reproducing
    size_t TIME_A = 5;
 
    size_t PREDATOR = 1, PREY = 0;
    size_t ALIVE = 1, DEAD = 0;
 
    // Predators reproducing probability
    float PRED_REPR = 0.1;
 
    // Predators eating probability
    float PRED_EAT = 0.2;
 
    // Prey reproducing probability
    float PREY_REPR = 0.1;
 
    // set the seed
    // create the random generator engine
    std::srand(v_cl.getProcessUnitID());
    std::default_random_engine eg;
    std::uniform_real_distribution<float> ud(0.0f, 1.0f);
 
    size_t k = 50000;
 
    print_test_v( "Testing 2D random walk vector k<=",k);
 
    BOOST_TEST_CHECKPOINT( "Testing 2D random walk k=" << k );
 
    Box<2,float> box({0.0,0.0},{1.0,1.0});
 
    // Grid info
    grid_sm<2, void> info( { GS_SIZE, GS_SIZE });
 
    // Boundary conditions
    size_t bc[2] = { PERIODIC, PERIODIC };
 
    // factor
    float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f);
 
    // interaction radius
    float r_cut = 0.01 / factor;
 
    // ghost
    Ghost<2,float> ghost(0.01 / factor);
 
    // Distributed vector
    vector_dist<2,float, animal, CartDecomposition<2, float, HeapMemory, ParMetisDistribution<2, float>>> vd(k,box,bc,ghost);
 
    // Init DLB tool
    DLB dlb(v_cl);
 
    // Set unbalance threshold
    dlb.setHeurisitc(DLB::Heuristic::UNBALANCE_THRLD);
    dlb.setThresholdLevel(DLB::ThresholdLevel::THRLD_MEDIUM);
 
    auto it = vd.getIterator();
 
    size_t c = 0;
    while (it.isNext())
    {
        auto key = it.get();
        if(c % 3)
        {
            vd.template getPos<animal::pos>(key)[0] = ud(eg);
            vd.template getPos<animal::pos>(key)[1] = ud(eg);
            vd.template getProp<animal::genre>(key) = 0; //prey
            vd.template getProp<animal::status>(key) = 1; //alive
            vd.template getProp<animal::time_a>(key) = TIME_A; //alive
        }else{
            vd.template getPos<animal::pos>(key)[0] = ud(eg);
            vd.template getPos<animal::pos>(key)[1] = ud(eg);
            vd.template getProp<animal::genre>(key) = 1; //predator
            vd.template getProp<animal::status>(key) = 1; //alive
            vd.template getProp<animal::time_a>(key) = TIME_A; //alive
        }
        ++it;
        ++c;
    }
 
    vd.map();
 
    vd.addComputationCosts();
 
    vd.getDecomposition().rebalance(dlb);
 
    vd.map();
 
    vd.getDecomposition().write("dec_init");
    vd.getDecomposition().getDistribution().write("parmetis_random_walk_" + std::to_string(0) + ".vtk");
    vd.write("particles_", 0, NO_GHOST);
 
    // 10 step random walk
    for (size_t j = 0; j < 50; j++)
    {
        std::cout << "Iteration " << (j+1) << "\n";
 
        auto it = vd.getDomainIterator();
 
        while (it.isNext())
        {
            auto key = it.get();
 
            vd.template getPos<animal::pos>(key)[0] += 0.01 * ud(eg);
            vd.template getPos<animal::pos>(key)[1] += 0.01 * ud(eg);
 
            ++it;
        }
 
        vd.map();
 
 
        vd.ghost_get<0>();
 
        openfpm::vector<size_t> deads;
 
        // get the cell list with a cutoff radius
 
        bool error = false;
 
        auto NN = vd.getCellList(0.01 / factor);
 
        // iterate across the domain particle
 
        auto it2 = vd.getDomainIterator();
 
        while (it2.isNext())
        {
            auto p = it2.get();
 
            Point<2,float> xp = vd.getPos<0>(p);
 
            size_t gp = vd.getProp<animal::genre>(p);
            size_t sp = vd.getProp<animal::status>(p);
 
            auto Np = NN.getIterator(NN.getCell(vd.getPos<0>(p)));
 
            while (Np.isNext())
            {
                auto q = Np.get();
 
                size_t gq = vd.getProp<animal::genre>(q);
                size_t sq = vd.getProp<animal::status>(q);
 
                // repulsive
 
                Point<2,float> xq = vd.getPos<0>(q);
                Point<2,float> f = (xp - xq);
 
                float distance = f.norm();
 
                //if p is a fox and q a rabit and they are both alive then the fox eats the rabbit
                if (distance < 2*r_cut*sqrt(2) && sp == ALIVE)
                {
                    if(gp == PREDATOR && gq == PREY && sq == ALIVE)
                    {
                        vd.getProp<animal::status>(q) = DEAD;
                        vd.getProp<animal::time_a>(q) = TIME_A;
                    }
                    else if (gp == PREY && gq == PREY && sq != DEAD)
                    {
                        vd.add();
                        //vd.getLastPos<animal::pos>()[0] = ud(eg);
                        //vd.getLastPos<animal::pos>()[1] = ud(eg);
                        vd.getLastProp<animal::genre>() = 0;
                    }
                }
 
                ++Np;
            }
 
            if(vd.getProp<animal::status>(p) == DEAD)
            {
                deads.add(p.getKey());
            }
 
            ++it2;
        }
 
        vd.remove(deads, 0);
        deads.resize(0);
 
 
        vd.addComputationCosts();
 
        vd.getDecomposition().rebalance(dlb);
 
        vd.map();
 
        vd.getDecomposition().getDistribution().write("parmetis_random_walk_" + std::to_string(j+1) + ".vtk");
        vd.write("particles_", j+1, NO_GHOST);
        vd.getDecomposition().write("dec_");
 
        size_t l = vd.size_local();
        v_cl.sum(l);
        v_cl.execute();
 
        // Count the local particles and check that the total number is consistent
        //size_t cnt = total_n_part_lc(vd,bc);
 
        //BOOST_REQUIRE_EQUAL((size_t)k,cnt);
    }
 
}
 
BOOST_AUTO_TEST_SUITE_END()
 
#endif /* SRC_TEST_ISOLATION_HPP_ */