unittests.cpp

//
//  Florian Probst
//  E-Mail: probstf@informatik.uni-freiburg.de / derbalvald@gmail.com
//

#include "unittests.h"

#include <execution>
#include <iostream>
#include "grid.h"
#include "particles.h"
#include <chrono>
#include "lbvh.h"

/*
int testuniformgrid()
{
    int h_init = 2;
    int dims = 3;
    int num_layers = 3;
    int num_particles = 27;

    Eigen::Array3Xd pos = Eigen::Array3Xd::Zero(dims, num_particles);
    Eigen::Array3Xd vel = Eigen::Array3Xd::Zero(dims, num_particles);
    Eigen::ArrayXd mass = Eigen::ArrayXd::Ones(num_particles) * 4.0;

    Particles particles = Particles(dims, pos, vel, mass, h_init);

    int index = 0;

    for (int x = 0; x < num_layers; x++)
    {
        for (int y = 0; y < num_layers; y++)
        {
            for (int z = 0; z < num_layers; z++)
            {
                if (index >= num_particles)
                {
                    break;
                }

                particles.positions.col(index) << x * h_init, y * h_init, z * h_init;

                double random1 = (static_cast<double>(rand()) / (RAND_MAX)) * 2.0;
                double random2 = (static_cast<double>(rand()) / (RAND_MAX)) * 2.0;
                double random3 = (static_cast<double>(rand()) / (RAND_MAX)) * 2.0;

                particles.accelerations.col(index) << random1, random2, random3;
                particles.velocities.col(index) << 0.0, 0.0, 0.0;
                particles.is_boundary(index) = false;
                index++;

            }
        }
    }

    grid u_grid(particles);
    u_grid.build_grid(particles);
    u_grid.compute_hashes(particles);
    u_grid.hash_sort(particles);
    u_grid.compute_cellStartEnd(particles);
    u_grid.get_neighbors(particles);

    for (auto particle : particles.neighbors)
    {
        std::cout << particle.size() << std::endl;
    }

    std::cout << particles.neighbors.size() << std::endl;

    return 0;
}


int testspatialhashlist()
{
    int h_init = 2;
    int dims = 3;
    int num_layers = 3;
    int num_particles = 27;
    int num_all_particles = num_particles;

    std::cout << "TEST1" << std::endl;

    Eigen::Array3Xd pos = Eigen::Array3Xd::Zero(dims, num_all_particles);
    Eigen::Array3Xd vel = Eigen::Array3Xd::Zero(dims, num_all_particles);
    Eigen::ArrayXd mass = Eigen::ArrayXd::Ones(num_all_particles) * 4.0;

    std::cout << "TEST2" << std::endl;

    Particles particles = Particles(dims, pos, vel, mass, h_init);

    std::cout << "TEST3" << std::endl;

    int index = 0;

    for (int x = 0; x < num_layers; x++)
    {
        for (int y = 0; y < num_layers; y++)
        {
            for (int z = 0; z < num_layers; z++)
            {
                if (index >= num_particles)
                {
                    break;
                }

                particles.positions.col(index) << x * h_init, y * h_init, z * h_init;

                double random1 = (static_cast<double>(rand()) / (RAND_MAX)) * 2.0;
                double random2 = (static_cast<double>(rand()) / (RAND_MAX)) * 2.0;
                double random3 = (static_cast<double>(rand()) / (RAND_MAX)) * 2.0;

                particles.accelerations.col(index) << random1, random2, random3;
                particles.velocities.col(index) << 0.0, 0.0, 0.0;
                particles.is_boundary(index) = false;
                index++;

            }
        }
    }

    std::cout << "TEST4" << std::endl;

    spatial_hash_list hash_list(11);

    std::cout << "TEST5" << std::endl;

    hash_list.build_hash_list(particles);

    std::cout << "TEST6" << std::endl;

    hash_list.get_neighbors(particles);

    std::cout << "TEST7" << std::endl;

    for (auto particle : particles.neighbors)
    {
        std::cout << particle.size() << std::endl;
    }

    std::cout << particles.neighbors.size() << std::endl;

    return 0;
}*/

int test_base()
{
    int h_init = 2;
    int dims = 3;
    int num_layers = 3;
    int num_particles = 27;
    int num_all_particles = num_particles;

    double rest_density = 1.0;

    std::cout << "TEST_BASE" << std::endl;

    Eigen::Array3Xd pos = Eigen::Array3Xd::Zero(dims, num_all_particles);
    Eigen::Array3Xd vel = Eigen::Array3Xd::Zero(dims, num_all_particles);
    Eigen::ArrayXd mass = Eigen::ArrayXd::Ones(num_all_particles) * 4.0;

    std::cout << "TEST_BASE" << std::endl;

    Particles particles = Particles(dims, pos, vel, mass, h_init, rest_density);

    std::cout << "TEST_BASE" << std::endl;

    int index = 0;

    for (int x = 0; x < num_layers; x++)
    {
        for (int y = 0; y < num_layers; y++)
        {
            for (int z = 0; z < num_layers; z++)
            {
                if (index >= num_particles)
                {
                    break;
                }

                particles.positions.col(index) << x * h_init, y * h_init, z * h_init;

                double random1 = (static_cast<double>(rand()) / (RAND_MAX)) * 2.0;
                double random2 = (static_cast<double>(rand()) / (RAND_MAX)) * 2.0;
                double random3 = (static_cast<double>(rand()) / (RAND_MAX)) * 2.0;

                particles.accelerations.col(index) << random1, random2, random3;
                particles.velocities.col(index) << 0.0, 0.0, 0.0;
                particles.is_boundary(index) = false;
                index++;

            }
        }
    }

    // find neighbors by just checking distance between particles being less than 2*h

    for (int i = 0; i < particles.positions.cols(); i++)
    {
        for (int j = 0; j < particles.positions.cols(); j++)
        {
            Eigen::Vector3d dist = particles.positions.col(i) - particles.positions.col(j);

            if (dist.norm() < particles.h*2 + 1e-6)
            {
                particles.neighbors[i].push_back(j);
            }
        }
    }

    //
    std::cout << particles.neighbors.size() << std::endl;

    for (int i = 0; i < particles.neighbors.size(); i++)
    {
        std::cout << particles.neighbors[i].size() << std::endl;


    }

    return 0;
}



int main_2() {

    int h_init = 1;
    int dims = 3;
    int num_layers = 200;
    int num_particles = num_layers * num_layers * num_layers;
    int num_all_particles = num_particles;

    double rest_density = 1.0;

    std::cout << "TEST_BASE" << std::endl;

    Eigen::Array3Xd pos = Eigen::Array3Xd::Zero(dims, num_all_particles);
    Eigen::Array3Xd vel = Eigen::Array3Xd::Zero(dims, num_all_particles);
    Eigen::ArrayXd mass = Eigen::ArrayXd::Ones(num_all_particles) * 4.0;

    std::cout << "TEST_BASE" << std::endl;

    Particles particles = Particles(dims, pos, vel, mass, h_init, rest_density);

    std::cout << "TEST_BASE" << std::endl;

    int index = 0;

    for (int x = 0; x < num_layers; x++)
    {
        for (int y = 0; y < num_layers; y++)
        {
            for (int z = 0; z < num_layers; z++)
            {
                if (index >= num_particles)
                {
                    break;
                }

                particles.positions.col(index) << x * h_init, y * h_init, z * h_init;

                double random1 = (static_cast<double>(rand()) / (RAND_MAX)) * 2.0;
                double random2 = (static_cast<double>(rand()) / (RAND_MAX)) * 2.0;
                double random3 = (static_cast<double>(rand()) / (RAND_MAX)) * 2.0;

                particles.accelerations.col(index) << random1, random2, random3;
                particles.velocities.col(index) << 0.0, 0.0, 0.0;
                particles.is_boundary(index) = false;
                index++;

            }
        }
    }

    double radius = h_init;
    double cellSize = radius * 2;

    std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();

    UniformGrid ug;
    auto test = particles.positions.rowwise().minCoeff();
    std::cout << "min = " << test << "\n";
    Eigen::Vector3d minCorner = Eigen::Vector3d(test);
    // Eigen::Vector3d offset = Eigen::Vector3d(particles.h*2, particles.h*2, particles.h*2);
    // minCorner -= offset;
    std::cout << "minCorner = " << minCorner << "\n";
    Eigen::Vector3d maxCorner = particles.positions.rowwise().maxCoeff();
    // maxCorner += offset;
    std::cout << "maxCorner = " << maxCorner << "\n";
    ug.init(minCorner, maxCorner, cellSize);
    ug.build(particles.positions);
    auto test_length = (radius * 2) + 1e-6;
    // Query neighbors for all particles and gather simple statistics (OpenMP parallelized)
    {
        long long total = 0;
        int maxn = 0;
#ifdef _OPENMP
        #pragma omp parallel
        {
            long long local_total = 0;
            int local_max = 0;
            #pragma omp for schedule(dynamic)
            for (int i = 0; i < num_particles; ++i)
            {
                std::vector<int> neigh = ug.queryNeighbors(particles.positions, i, test_length);
                local_total += static_cast<long long>(neigh.size());
                if (static_cast<int>(neigh.size()) > local_max) local_max = static_cast<int>(neigh.size());
            }
            #pragma omp atomic
            total += local_total;
            #pragma omp critical
            if (local_max > maxn) maxn = local_max;
        }
#else
        for (int i = 0; i < num_particles; ++i)
        {
            std::vector<int> neigh = ug.queryNeighbors(particles.positions, i, test_length);
            total += static_cast<long long>(neigh.size());
            if (static_cast<int>(neigh.size()) > maxn) maxn = static_cast<int>(neigh.size());
        }
#endif
        std::cout << "UniformGrid totalNeighbors = " << total
                  << " avg = " << (double)total / num_particles
                  << " max = " << maxn << "\n";
    }

    std::chrono::steady_clock::time_point end1 = std::chrono::steady_clock::now();

    SpatialHashGrid sh(cellSize);
    sh.build(particles.positions);
    {
        long long total = 0;
        int maxn = 0;
#ifdef _OPENMP
        #pragma omp parallel
        {
            long long local_total = 0;
            int local_max = 0;
            #pragma omp for schedule(dynamic)
            for (int i = 0; i < num_particles; ++i)
            {
                std::vector<int> neigh = sh.queryNeighbors(particles.positions, i, test_length);
                local_total += static_cast<long long>(neigh.size());
                if (static_cast<int>(neigh.size()) > local_max) local_max = static_cast<int>(neigh.size());
            }
            #pragma omp atomic
            total += local_total;
            #pragma omp critical
            if (local_max > maxn) maxn = local_max;
        }
#else
        for (int i = 0; i < num_particles; ++i)
        {
            std::vector<int> neigh = sh.queryNeighbors(particles.positions, i, test_length);
            total += static_cast<long long>(neigh.size());
            if (static_cast<int>(neigh.size()) > maxn) maxn = static_cast<int>(neigh.size());
        }
#endif
        std::cout << "SpatialHashGrid totalNeighbors = " << total
                  << " avg = " << (double)total / num_particles
                  << " max = " << maxn << "\n";
    }

    std::chrono::steady_clock::time_point end2 = std::chrono::steady_clock::now();

    SortedHashGrid sg(cellSize);
    sg.build(particles.positions);
    {
        long long total = 0;
        int maxn = 0;
#ifdef _OPENMP
        #pragma omp parallel
        {
            long long local_total = 0;
            int local_max = 0;
            #pragma omp for schedule(dynamic)
            for (int i = 0; i < num_particles; ++i)
            {
                std::vector<int> neigh = sg.queryNeighbors(particles.positions, i, test_length);
                local_total += static_cast<long long>(neigh.size());
                if (static_cast<int>(neigh.size()) > local_max) local_max = static_cast<int>(neigh.size());
            }
            #pragma omp atomic
            total += local_total;
            #pragma omp critical
            if (local_max > maxn) maxn = local_max;
        }
#else
        for (int i = 0; i < num_particles; ++i)
        {
            std::vector<int> neigh = sg.queryNeighbors(particles.positions, i, test_length);
            total += static_cast<long long>(neigh.size());
            if (static_cast<int>(neigh.size()) > maxn) maxn = static_cast<int>(neigh.size());
        }
#endif
        std::cout << "SortedHashGrid totalNeighbors = " << total
                  << " avg = " << (double)total / num_particles
                  << " max = " << maxn << "\n";
    }

    std::chrono::steady_clock::time_point end3 = std::chrono::steady_clock::now();

    SpatialHashGrid_Ihmsen shi(cellSize, num_particles);
    shi.build(particles.positions);
    {
        long long total = 0;
        int maxn = 0;
#ifdef _OPENMP
        #pragma omp parallel
        {
            long long local_total = 0;
            int local_max = 0;
            #pragma omp for schedule(dynamic)
            for (int i = 0; i < num_particles; ++i)
            {
                std::vector<int> tmp;
                shi.queryNeighbors(particles.positions, i, test_length, tmp);
                local_total += static_cast<long long>(tmp.size());
                if (static_cast<int>(tmp.size()) > local_max) local_max = static_cast<int>(tmp.size());
            }
            #pragma omp atomic
            total += local_total;
            #pragma omp critical
            if (local_max > maxn) maxn = local_max;
        }
#else
        std::vector<int> tmp;
        for (int i = 0; i < num_particles; ++i)
        {
            shi.queryNeighbors(particles.positions, i, test_length, tmp);
            total += static_cast<long long>(tmp.size());
            if (static_cast<int>(tmp.size()) > maxn) maxn = static_cast<int>(tmp.size());
        }
#endif
        std::cout << "SpatialHashGrid_Ihmsen totalNeighbors = " << total
                  << " avg = " << (double)total / num_particles
                  << " max = " << maxn << "\n";
    }

    std::chrono::steady_clock::time_point end4 = std::chrono::steady_clock::now();

    SortedHashGrid_Ihmsen sgi(cellSize, num_particles);
    sgi.build(particles.positions);
    {
        long long total = 0;
        int maxn = 0;
#ifdef _OPENMP
        #pragma omp parallel
        {
            long long local_total = 0;
            int local_max = 0;
            #pragma omp for schedule(dynamic)
            for (int i = 0; i < num_particles; ++i)
            {
                std::vector<int> neigh = sgi.queryNeighbors(particles.positions, i, test_length);
                local_total += static_cast<long long>(neigh.size());
                if (static_cast<int>(neigh.size()) > local_max) local_max = static_cast<int>(neigh.size());
            }
            #pragma omp atomic
            total += local_total;
            #pragma omp critical
            if (local_max > maxn) maxn = local_max;
        }
#else
        for (int i = 0; i < num_particles; ++i)
        {
            std::vector<int> neigh = sgi.queryNeighbors(particles.positions, i, test_length);
            total += static_cast<long long>(neigh.size());
            if (static_cast<int>(neigh.size()) > maxn) maxn = static_cast<int>(neigh.size());
        }
#endif
        std::cout << "SortedHashGrid_Ihmsen totalNeighbors = " << total
                  << " avg = " << (double)total / num_particles
                  << " max = " << maxn << "\n";
    }

    std::chrono::steady_clock::time_point end5 = std::chrono::steady_clock::now();

    UniformGrid_Morton ug_m;
    minCorner = particles.positions.rowwise().minCoeff();
    maxCorner = particles.positions.rowwise().maxCoeff();
    ug_m.init(minCorner, maxCorner, cellSize);
    ug_m.build(particles.positions);
    {
        long long total = 0;
        int maxn = 0;
#ifdef _OPENMP
        #pragma omp parallel
        {
            long long local_total = 0;
            int local_max = 0;
            #pragma omp for schedule(dynamic)
            for (int i = 0; i < num_particles; ++i)
            {
                std::vector<int> neigh = ug_m.queryNeighbors(particles.positions, i, test_length);
                local_total += static_cast<long long>(neigh.size());
                if (static_cast<int>(neigh.size()) > local_max) local_max = static_cast<int>(neigh.size());
            }
            #pragma omp atomic
            total += local_total;
            #pragma omp critical
            if (local_max > maxn) maxn = local_max;
        }
#else
        for (int i = 0; i < num_particles; ++i)
        {
            std::vector<int> neigh = ug_m.queryNeighbors(particles.positions, i, test_length);
            total += static_cast<long long>(neigh.size());
            if (static_cast<int>(neigh.size()) > maxn) maxn = static_cast<int>(neigh.size());
        }
#endif
        std::cout << "UniformGrid_Morton totalNeighbors = " << total
                  << " avg = " << (double)total / num_particles
                  << " max = " << maxn << "\n";
    }

    std::chrono::steady_clock::time_point end6 = std::chrono::steady_clock::now();

    // LinearBVH bvh(cellSize);
    // bvh.build(particles.positions);
    // std::cout << "LinearBVH neighbors[0] = " << bvh.queryNeighbors(particles.positions, 0, test_length).size() << "\n";

    std::chrono::steady_clock::time_point end7 = std::chrono::steady_clock::now();

    // --- LBVH (from lbvh.h) ---
    LBVH lbvh;
    lbvh.build(
        num_particles,
        [&](int i) -> LBVH::Vec3
        {
            return LBVH::Vec3(
                static_cast<LBVH::Scalar>(particles.positions(0, i)),
                static_cast<LBVH::Scalar>(particles.positions(1, i)),
                static_cast<LBVH::Scalar>(particles.positions(2, i))
            );
        },
        [&](int primId, LBVH::Vec3& outMin, LBVH::Vec3& outMax)
        {
            // Points: AABB degenerates to the point itself (still works for sphere queries)
            outMin = LBVH::Vec3(
                static_cast<LBVH::Scalar>(particles.positions(0, primId)),
                static_cast<LBVH::Scalar>(particles.positions(1, primId)),
                static_cast<LBVH::Scalar>(particles.positions(2, primId))
            );
            outMax = outMin;
        }
    );

    {
        long long total = 0;
        int maxn = 0;
#ifdef _OPENMP
        #pragma omp parallel
        {
            long long local_total = 0;
            int local_max = 0;
            #pragma omp for schedule(dynamic)
            for (int i = 0; i < num_particles; ++i)
            {
                std::vector<int> lbvhNeighbors;
                const LBVH::Vec3 center(
                    static_cast<LBVH::Scalar>(particles.positions(0, i)),
                    static_cast<LBVH::Scalar>(particles.positions(1, i)),
                    static_cast<LBVH::Scalar>(particles.positions(2, i))
                );
                lbvh.queryNeighbors(
                    center,
                    static_cast<LBVH::Scalar>(test_length),
                    [&](int primId) -> LBVH::Vec3 {
                        return LBVH::Vec3(
                            static_cast<LBVH::Scalar>(particles.positions(0, primId)),
                            static_cast<LBVH::Scalar>(particles.positions(1, primId)),
                            static_cast<LBVH::Scalar>(particles.positions(2, primId))
                        );
                    },
                    lbvhNeighbors
                );
                local_total += static_cast<long long>(lbvhNeighbors.size());
                if (static_cast<int>(lbvhNeighbors.size()) > local_max) local_max = static_cast<int>(lbvhNeighbors.size());
            }
            #pragma omp atomic
            total += local_total;
            #pragma omp critical
            if (local_max > maxn) maxn = local_max;
        }
#else
        for (int i = 0; i < num_particles; ++i)
        {
            std::vector<int> lbvhNeighbors;
            const LBVH::Vec3 center(
                static_cast<LBVH::Scalar>(particles.positions(0, i)),
                static_cast<LBVH::Scalar>(particles.positions(1, i)),
                static_cast<LBVH::Scalar>(particles.positions(2, i))
            );
            lbvh.queryNeighbors(
                center,
                static_cast<LBVH::Scalar>(test_length),
                [&](int primId) -> LBVH::Vec3 {
                    return LBVH::Vec3(
                        static_cast<LBVH::Scalar>(particles.positions(0, primId)),
                        static_cast<LBVH::Scalar>(particles.positions(1, primId)),
                        static_cast<LBVH::Scalar>(particles.positions(2, primId))
                    );
                },
                lbvhNeighbors
            );
            total += static_cast<long long>(lbvhNeighbors.size());
            if (static_cast<int>(lbvhNeighbors.size()) > maxn) maxn = static_cast<int>(lbvhNeighbors.size());
        }
#endif
        std::cout << "LBVH totalNeighbors = " << total
                  << " avg = " << (double)total / num_particles
                  << " max = " << maxn << "\n";
    }

    std::chrono::steady_clock::time_point end8 = std::chrono::steady_clock::now();

    // write out the times it needed to build and query the datastructure:

    std::cout << "ug: " << std::chrono::duration_cast<std::chrono::milliseconds>(end1 - begin).count() << "\n"
              << "sphg: " << std::chrono::duration_cast<std::chrono::milliseconds>(end2 - end1).count() << "\n"
              << "sohg: " << std::chrono::duration_cast<std::chrono::milliseconds>(end3 - end2).count() << "\n"
              << "sphg_i: " << std::chrono::duration_cast<std::chrono::milliseconds>(end4 - end3).count() << "\n"
              << "sohg_i: " << std::chrono::duration_cast<std::chrono::milliseconds>(end5 - end4).count() << "\n"
              << "ug_m: " << std::chrono::duration_cast<std::chrono::milliseconds>(end6 - end5).count() << "\n"
              // << "bvh: " << std::chrono::duration_cast<std::chrono::milliseconds>(end7 - end6).count() << "\n"
              << "lbvh: " << std::chrono::duration_cast<std::chrono::milliseconds>(end8 - end7).count()
              << std::endl;

    return 0;
}


int morten_encode_decode_test()
{
    CellCoord coord1 = CellCoord(7, 2097150, 2097151);
    CellCoord coord2 = CellCoord(66666, 3, 1230);
    CellCoord coord3 = CellCoord(23, 234, 123);

    // Range for cell coord used to be 65535. its now much higher.
    CellCoord special = CellCoord(65535, 234, 123);

    auto morton1 = encode_morton(coord1.x, coord1.y, coord1.z);
    auto morton2 = encode_morton(coord2.x, coord2.y, coord2.z);
    auto morton3 = encode_morton(coord3.x, coord3.y, coord3.z);

    auto morton_special = encode_morton(special.x, special.y, special.z);

    CellCoord test_coord1 = decode_morton(morton1);
    CellCoord test_coord2 = decode_morton(morton2);
    CellCoord test_coord3 = decode_morton(morton3);

    CellCoord test_special = decode_morton(morton_special);

    std::cout << "special x: " << special.x << "\n";
    std::cout << "special y: " << special.y << "\n";
    std::cout << "special z: " << special.z << "\n";
    std::cout << "morton special: " << morton_special << "\n";
    std::cout << "test_special x: " << test_special.x << "\n";
    std::cout << "test_special y: " << test_special.y << "\n";
    std::cout << "test_special z: " << test_special.z << "\n";

    std::cout << morton1 << "\n";
    std::cout << morton2 << "\n";
    std::cout << morton3 << "\n";

    std::cout << coord1.x << "\n";
    std::cout << coord1.y << "\n";
    std::cout << coord1.z << "\n";

    std::cout << test_coord1.x << "\n";
    std::cout << test_coord1.y << "\n";
    std::cout << test_coord1.z << "\n";

    assert(coord1.x == test_coord1.x);
    assert(coord1.y == test_coord1.y);
    assert(coord1.z == test_coord1.z);
    assert(coord2.x == test_coord2.x);
    assert(coord2.y == test_coord2.y);
    assert(coord2.z == test_coord2.z);
    assert(coord3.x == test_coord3.x);
    assert(coord3.y == test_coord3.y);
    assert(coord3.z == test_coord3.z);

    return 0;
}



int main()
{
    std::cout << "testing morton encoding and decoding" << std::endl;

    morten_encode_decode_test();

    std::cout << "Basis No Acceleration Structure" << std::endl;

    test_base();

    std::cout << "Testing Acceleration Structures" << std::endl;

    main_2();

    // std::cout << "Uniform Grid" << std::endl;

    // testuniformgrid();

    // std::cout << "Spatial Hash list" << std::endl;

    // testspatialhashlist();

}