Commit dd27ab76 authored by Sebastian J. Mielke's avatar Sebastian J. Mielke
Browse files

Adaptive CellList Nearest Neighbors works!

parent 1b3d1aea
......@@ -8,20 +8,6 @@
#include <string>
/*
* How does it work (normal cell lists)?
*
* I have a particle and want to know all its interactions with particles in its subdomain.
* Since I don't want to check against all other particles, I only check the proximity condition for those that are "nearby enough".
* We could do so by dividing the subdomain into many cells of equal size (maximal cutoff/radius of the particles as the length of each side).
* Once we accomplished that we only have to check for possible interactions with all the particles that are in neighbor cells of the cell our initially chosen particle was in.
* This is what the NNIterator does: given some start cell list all elements elements in this start cell could maybe interact with.
* There are multiple ways to define neighbor cells, I can check all neighbor cells (full) or only those that are "further ahaed" (sym).
* If we know that interactions are symmetric, this makes them unique (and saves us nearly half the computations!)
*
* First, let us screw all that fancy logic and implement it as naively as possible... so let us compute with all points! \o/
*/
// Stub implementation
template<unsigned int dim, typename T, unsigned int impl=BALANCED, typename ele_container=openfpm::vector<std::pair<Point<dim+1, T>, size_t>> >
......@@ -114,18 +100,41 @@ class AdaptiveCellList<dim,T,BALANCED,ele_container> // : public CellDecomposer_
}
}
size_t firstIndexOfLevel(size_t level) {
inline size_t firstIndexOfLevel(size_t level) {
//sum i from 0 to (target_level-1): 2^(dim*i) = (2^(dim+dim*(target_level-1)) - 1) / (2^dim - 1)
return ((1l << (dim*level)) - 1) / ((1l << dim) - 1);
}
public:
// Object type that the structure store
typedef T value_type;
inline size_t maxlevel() {return max_level;}
inline size_t size() {return all_eles.size();}
inline T getEdgeLengthForLevel(size_t level) {
return D_m / static_cast<T>(1l << level);
}
inline T getRadiusForLevel(size_t level) {
return getEdgeLengthForLevel(level+1) * static_cast<T>(1.5f); // the *1.5 should be useless. TODO: make sure that it is
}
/*! \warning Checks for < low bound and >= (!) high bound
*
*/
inline bool isInside(Point<dim+1,T> p) {
for(unsigned int i=0; i<dim; i++)
if(p.get(i) < sbox.getLow(i) || p.get(i) >= sbox.getHigh(i))
return false;
return true;
}
// For iterating over all points to compute all interactions:
inline decltype(all_eles.begin()) begin() {return all_eles.begin();}
inline decltype(all_eles.end()) end() {return all_eles.end();}
/*! Initialize the cell list
*
* \param sbox Domain where this cell list is living
......@@ -179,7 +188,7 @@ public:
/*! \brief Sorts elements according so they are coherent in cells
*
* TODO: use some sort of isvalid-bool that is unset in add(), set in construct() and checked in the iterator
* TODO: use some sort of isvalid-bool that is unset in add(), set in construct() and checked in the iterator-getter
*
* \note Has to be called after insertion and before usage!
*
......@@ -195,7 +204,11 @@ public:
//std::cout << "Min. radius: " << minradius << std::endl;
max_level = std::ceil(std::log2(D_m / minradius)) - 1;
// fun fact: more than about 64/dim levels are not possible with these indices (and thats already with 8 byte for a size_t). Perhaps set it as maximum? TODO.
// fun fact: more than about 64/dim levels are not even theoretically possible with these indices (and thats already with 8 byte for a size_t).
// even worse: when using float, we get inconsistencies starting at level 24.
// So we should think about setting a maximum here (and making sure that everything else honors this):
// TODO .
//std::cout << "Max. level: " << max_level << std::endl;
auto end_iter = all_eles.begin();
......@@ -337,7 +350,7 @@ public:
for(unsigned int d = 0; d < dim; d++)
pos.get(d) = sbox.getP1().get(d);
// cutoff radius
pos.get(dim) = (D_m / static_cast<T>(1l << (level+1))) * static_cast<T>(1.5f); // the *1.5 should be useless. TODO: make sure that it is
pos.get(dim) = getRadiusForLevel(level);
// using this offset as a pivot on the whole row gives me the position w.r.t to the highest dimension on the highest level
// then pivoting on this half of the row gives me the pos in the next highest dim and so on until
......@@ -363,18 +376,6 @@ public:
return pos;
}
/*! \brief Get an element in the cell
*
* \param ele_id element id
*
* \return The element value
*
*/
inline size_t& get(size_t ele_id)
{
return all_eles.get(ele_id).second;
}
/*! \brief Swap the memory
*
* \param cl Cell list with witch you swap the memory
......@@ -390,23 +391,23 @@ public:
* \param cell cell id
*
*/
template<unsigned int impl> inline AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,FULL,impl> getNNIterator(size_t cell)
template<unsigned int impl> inline AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,FULL,impl> getNNIterator(Point<dim+1,T> p)
{
AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,FULL,impl> cln(*this);
AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,FULL,impl> cln(p, *this);
return cln;
}
template<unsigned int impl> inline AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,SYM,impl> getNNIteratorSym(size_t cell)
template<unsigned int impl> inline AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,SYM,impl> getNNIteratorSym(Point<dim+1,T> p)
{
AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,SYM,impl> cln(*this);
AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,SYM,impl> cln(p, *this);
return cln;
}
template<unsigned int impl> inline AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,CRS,impl> getNNIteratorCross(size_t cell)
template<unsigned int impl> inline AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,CRS,impl> getNNIteratorCross(Point<dim+1,T> p)
{
AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,CRS,impl> cln(*this);
AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,CRS,impl> cln(p, *this);
return cln;
}
......
......@@ -8,6 +8,8 @@
#ifndef ADAPTIVECELLLISTNNITERATOR_HPP_
#define ADAPTIVECELLLISTNNITERATOR_HPP_
#include "util/mathutil.hpp"
#include <vector>
#include <iostream>
template<unsigned int dim, typename CellS, unsigned int NNc_size, unsigned int impl> class AdaptiveCellNNIterator
......@@ -15,17 +17,84 @@ template<unsigned int dim, typename CellS, unsigned int NNc_size, unsigned int i
CellS & cl;
size_t ele_id;
Point<dim+1, typename CellS::value_type> p;
size_t currentlevel;
std::vector<size_t> cellindices;
std::vector<size_t>::iterator cell_iter, cell_iter_end;
decltype(cl.getCellContents(0).first) ele_iter, ele_iter_end;
// TODO: theres got to be a more intelligent way of doing this.
inline void populateIndices(int current_dim, Point<dim+1, typename CellS::value_type>& p, typename CellS::value_type edgelength) {
if (current_dim == dim) {
if(cl.isInside(p))
cellindices.push_back(cl.findCellIndex(p));
//std::cout << p.toString() << std::endl;
}
else {
populateIndices(current_dim + 1, p, edgelength);
p.get(current_dim) -= edgelength;
populateIndices(current_dim + 1, p, edgelength);
p.get(current_dim) += edgelength+edgelength;
populateIndices(current_dim + 1, p, edgelength);
p.get(current_dim) -= edgelength;
}
}
inline void updateCellIterators() {
cellindices.resize(0);
cellindices.reserve(openfpm::math::pow(3,dim));
Point<dim+1, typename CellS::value_type> center(p);
center.get(dim) = cl.getRadiusForLevel(currentlevel);
typename CellS::value_type edgelength = cl.getEdgeLengthForLevel(currentlevel);
populateIndices(0,center,edgelength);
/*
std::cout << " -- On level " << currentlevel << ": ";
for (auto i : cellindices)
std::cout << i << ", ";
std::cout << std::endl;
//*/
cell_iter = cellindices.begin();
cell_iter_end = cellindices.end();
}
inline void updateEleIterators() {
//std::cout << "Request ele_iters for cell " << *cell_iter << std::endl;
auto iterpair = cl.getCellContents(*cell_iter);
ele_iter = iterpair.first;
ele_iter_end = iterpair.second;
//std::cout << "It contains " << std::distance(ele_iter, ele_iter_end) << " elements" << std::endl;
}
public:
/*! \brief
*
* Cell NN iterator
*
*/
AdaptiveCellNNIterator(CellS & cl)
:cl(cl), ele_id(0)
AdaptiveCellNNIterator(Point<dim+1, typename CellS::value_type> point, CellS & cl)
:cl(cl), p(point), currentlevel(0)
{
//std::cout << "NN of " << point.toString() << std::endl;
updateCellIterators();
updateEleIterators();
while(ele_iter == ele_iter_end && currentlevel <= cl.maxlevel() + 1) {
++cell_iter;
if(cell_iter == cell_iter_end) {
++currentlevel;
updateCellIterators();
}
updateEleIterators();
}
//std::cout << "done constructing with cell_iter at " << *cell_iter << std::endl;
}
/*! \brief
......@@ -35,7 +104,8 @@ public:
*/
bool isNext()
{
return ele_id < cl.size();
//std::cout << currentlevel << " < " << cl.maxlevel() << std::endl;
return currentlevel <= cl.maxlevel();
}
/*! \brief take the next element
......@@ -43,7 +113,18 @@ public:
*/
AdaptiveCellNNIterator & operator++()
{
if(ele_id < cl.size()) ele_id++;
// Condition: don't start in an empty cell, otherwise this `++` will trick the following `if`.
++ele_iter;
while(ele_iter == ele_iter_end && currentlevel <= cl.maxlevel() + 1) {
++cell_iter;
if(cell_iter == cell_iter_end) {
++currentlevel;
updateCellIterators();
}
updateEleIterators();
}
return *this;
}
......@@ -52,9 +133,10 @@ public:
* \return the actual element
*
*/
inline size_t& get()
inline std::pair<Point<dim+1,typename CellS::value_type>, size_t>& get()
{
return cl.get(ele_id);
//std::cout << "serving " << ele_iter->first.toString() << std::endl << std::flush;
return *ele_iter;
}
};
......
......@@ -28,6 +28,7 @@ Point<dim+1,double> giveRadius(Point<dim,double> const & p)
*/
template<unsigned int dim, typename T, typename CellS> void Test_ar_cell_s()
{
/*
//Space where is living the Cell list
SpaceBox<dim,T> box({0.0f,0.0f,0.0f},{1.0f,1.0f,1.0f});
......@@ -39,30 +40,32 @@ template<unsigned int dim, typename T, typename CellS> void Test_ar_cell_s()
// Iterate through each element
Point<dim,T> key({0.1,0.1,0.1});
cl1.add(giveRadius(key),13);
Point<dim,T> key({0.1,0.1,0.1});
cl1.add(giveRadius(key),13);
// particle two
key = {0.9,0.9,0.9};
cl1.add(giveRadius(key),42);
// particle two
key = {0.9,0.9,0.9};
cl1.add(giveRadius(key),42);
// check the cells are correctly filled
// Check we have 2 objects
auto NN = cl1.template getNNIterator<NO_CHECK>(4711);
size_t total = 0;
auto NN = cl1.template getNNIterator<NO_CHECK>(4711);
size_t total = 0;
while(NN.isNext())
{
// total
while(NN.isNext())
{
// total
total++;
total++;
++NN;
}
++NN;
}
BOOST_REQUIRE_EQUAL(total,2);
BOOST_REQUIRE_EQUAL(total,2);
//*/
}
BOOST_AUTO_TEST_SUITE( AdaptiveCellList_test )
......
......@@ -18,7 +18,16 @@
*
* has complexity O(1) in getting the cell id and the elements in a cell
* but with different constants
*
*
* So how does it work?
*
* I have a particle and want to know all its interactions with particles in its subdomain.
* Since I don't want to check against all other particles, I only check the proximity condition for those that are "nearby enough".
* We could do so by dividing the subdomain into many cells of equal size (maximal cutoff/radius of the particles as the length of each side).
* Once we accomplished that we only have to check for possible interactions with all the particles that are in neighbor cells of the cell our initially chosen particle was in.
* This is what the NNIterator does: given some start cell, list all elements that elements in this start cell could maybe interact with.
* There are multiple ways to define neighbor cells, I can check all neighbor cells (full) or only those that are "further ahaed" (sym).
* If we know that interactions are symmetric, this makes them unique (and saves us nearly half the computations!)
*/
// Fastest implementation of the Cell list
......
......@@ -7,6 +7,27 @@
#include <iostream>
#include <boost/mpl/int.hpp>
#include <typeinfo>
#include <unordered_set>
// http://stackoverflow.com/questions/20590656/error-for-hash-function-of-pair-of-ints
struct pairhash {
public:
template <typename T, typename U>
std::size_t operator()(const std::pair<T, U> &x) const
{
return (3 * std::hash<T>()(x.first)) ^ std::hash<U>()(x.second);
}
};
// Include tests
......@@ -118,6 +139,7 @@ std::vector<Point<dim, float>> getSamplePoints() {
std::vector<Point<2, float>> samplepoints2 = getSamplePoints<2>();
std::vector<Point<3, float>> samplepoints3 = getSamplePoints<3>();
std::unordered_multiset<std::pair<size_t, size_t>, pairhash> allInteractions;
template<int dim, class CellS>
void insertIntoCl(CellS &cl, std::vector<Point<dim, float>> &points) {
......@@ -129,8 +151,20 @@ void insertIntoCl(CellS &cl, std::vector<Point<dim, float>> &points) {
}
}
AdaptiveCellList<2, float> getSampleARCL() {
SpaceBox<2,float> box({0.0f,0.0f},{1.0f,1.0f});
Point<2,float> org({0.0,0.0});
AdaptiveCellList<2, float> arcl(box,org);
insertIntoCl<3, AdaptiveCellList<2, float>>(arcl, samplepoints3);
arcl.construct();
return arcl;
}
BOOST_AUTO_TEST_CASE( celllist_realloc)
{
return;
std::cout << "2D Cell List:" << std::endl;
//Space where is living the Cell list
......@@ -223,12 +257,13 @@ inline bool is_in_radius(const Point<3, float> &p1, const Point<3, float> &p2) {
BOOST_AUTO_TEST_CASE( get_all_interactions_classiccelllist)
{
return;
//return;
timestamp_t c0 = get_timestamp();
auto maxradiusiterator = std::max_element(samplepoints3.begin(), samplepoints3.end(), [](Point<3, float>& a, Point<3, float>& b){return a[2] < b[2];});
size_t gridsize = std::floor(1.0f / (*maxradiusiterator)[2]);
//gridsize = 7;
std::cout << "Please choose a gridsize < " << 1.0f / (*maxradiusiterator)[2] << ": " << gridsize << std::endl;
//std::cout << "Please choose a gridsize < " << 1.0f / (*maxradiusiterator)[2] << ": " << gridsize << std::endl;
const float epsilon = 0.0001f; // just a little bit bigger, so 1.0 is still inside.
//const float epsilon = 0.0f; // Iterating over padding cells crashes... due to missing boundary checks, I guess?
......@@ -241,11 +276,15 @@ BOOST_AUTO_TEST_CASE( get_all_interactions_classiccelllist)
insertIntoCl<2, CellList<2, float>>(cl, samplepoints2); // shares indices with ~3, so just use the 2d one for a simple celllist
timestamp_t c1 = get_timestamp();
std::cout << "Creation time (us): " << c1-c0 << std::endl;
size_t interactions_count = 0, interaction_candidates = 0;
size_t cell, i1, i2;
std::cout << "ready lets go" << std::endl;
timestamp_t t0 = get_timestamp();
for (int i=pad; i<div[0]+pad; i++) {
......@@ -261,9 +300,10 @@ BOOST_AUTO_TEST_CASE( get_all_interactions_classiccelllist)
i1 = iter_inner.get();
++interaction_candidates;
if(i1 != i2 && is_in_radius(samplepoints3[i1], samplepoints3[i2])) {
//std::cout << i1 << " and " << i2 << std::endl;
//std::cout << i1 << " and " << i2 << " with p2 = " << samplepoints3[i2].toString() << std::endl;
allInteractions.insert(std::make_pair(i1, i2));
++interactions_count;
}//else std::cout << i1 << " AND " << i2 << std::endl;
}
++iter_inner;
}
//std::cout << std::endl;
......@@ -273,64 +313,22 @@ BOOST_AUTO_TEST_CASE( get_all_interactions_classiccelllist)
}
timestamp_t t1 = get_timestamp();
std::cout << "that was easy (us): " << t1-t0 << std::endl;
std::cout << "found interactions: " << interactions_count
std::cout << "Interactions time (us): " << t1-t0 << std::endl;
std::cout << "Found interactions: " << interactions_count
<< " of " << interaction_candidates
<< " candidates (" << (static_cast<float>(interaction_candidates) / interactions_count)
<< "x)." << std::endl;
}
BOOST_AUTO_TEST_CASE( get_all_interactions_arcelllist)
BOOST_AUTO_TEST_CASE( find_inserted_items_in_arcl)
{
SpaceBox<2,float> box({0.0f,0.0f},{1.0f,1.0f});
Point<2,float> org({0.0,0.0});
AdaptiveCellList<2, float> arcl(box,org);
insertIntoCl<3, AdaptiveCellList<2, float>>(arcl, samplepoints3);
arcl.construct();
size_t interactions_count = 0, interaction_candidates = 0;
size_t i1, i2;
std::cout << "ready lets go (adaptive)" << std::endl;
timestamp_t t0 = get_timestamp();
///*
auto iter = arcl.getNNIterator<FAST>(4711); //full interactions
while (iter.isNext()) {
i2 = iter.get();
//std::cout << "Neighbor found: " << i2 << std::endl;
auto iter_inner = arcl.getNNIterator<FAST>(4711); //full interactions
while (iter_inner.isNext()) {
i1 = iter_inner.get();
++interaction_candidates;
if(i1 != i2 && is_in_radius(samplepoints3[i1], samplepoints3[i2])) {
++interactions_count;
//std::cout << i1 << " and " << i2 << std::endl;
}
++iter_inner;
}
//std::cout << std::endl;
++iter;
}
//*/
timestamp_t t1 = get_timestamp();
std::cout << "that was easy (us): " << t1-t0 << std::endl;
std::cout << "found interactions (should be 475560): " << interactions_count
<< " of " << interaction_candidates
<< " candidates (" << (static_cast<float>(interaction_candidates) / interactions_count)
<< "x)." << std::endl;
//auto printresult = arcl.printTree(0,0,0);
//std::cout << printresult.first << printresult.second << std::endl;
return;
auto arcl = getSampleARCL();
for(auto& p: samplepoints3) {
size_t cellindex = arcl.findCellIndex(p);
// check here, whether the calculated cell really contains our point! thats easy and will be a good base for debugging
// check here, whether the calculated cell really contains our point!
bool found = false;
auto iter = arcl.getCellContents(cellindex);
......@@ -341,39 +339,77 @@ BOOST_AUTO_TEST_CASE( get_all_interactions_arcelllist)
BOOST_REQUIRE(found);
}
/*
std::cout << "0: ";
for (auto& i : arcl.findChildrenIndices(0))
std::cout << i << " ";
std::cout << std::endl;
std::cout << "1: ";
for (auto& i : arcl.findChildrenIndices(1))
std::cout << i << " ";
std::cout << std::endl;
std::cout << "5: ";
for (auto& i : arcl.findChildrenIndices(5))
std::cout << i << " ";
std::cout << std::endl;
std::cout << "6: ";
for (auto& i : arcl.findChildrenIndices(6))
std::cout << i << " ";
std::cout << std::endl;
*/
//auto printresult = arcl.printTree(0,0,0);
//std::cout << printresult.first << printresult.second << std::endl;
}
BOOST_AUTO_TEST_CASE( findcellindex_and_findcellcenter_combine)
{
return;
auto arcl = getSampleARCL();
size_t max = -1; // = 18446744073709551615 on my system
size_t max = 164193736414550; // = theoretically we should choose -1, which is 18446744073709551615 on my system, but on high levels float fails, so we stop earlier.
size_t sqrt_max = sqrt(max);
//max = (1l << 32) + 1l;
for(size_t i = 0; i < 10000000; i++)
BOOST_REQUIRE_EQUAL(arcl.findCellIndex(arcl.findCellCenter(i)), i);
for(size_t i = 10000; i < max; i += 777777) {
std::cout << "Checked until " << 10000000 << std::endl;
for(size_t i = 100000000; i < sqrt_max+10000; i += 777777)
BOOST_REQUIRE_EQUAL(arcl.findCellIndex(arcl.findCellCenter(i)), i);
std::cout << "Checked until " << sqrt_max+100000 << std::endl;
for(size_t i = sqrt_max+10000; i < max; i += sqrt_max)
BOOST_CHECK_EQUAL(arcl.findCellIndex(arcl.findCellCenter(i)), i);
std::cout << "Checked until " << max << std::endl;
}
BOOST_AUTO_TEST_CASE( get_all_interactions_arcelllist)
{
timestamp_t c0 = get_timestamp();
auto arcl = getSampleARCL();
timestamp_t c1 = get_timestamp();
std::cout << "Creation time (us): " << c1-c0 << std::endl;
size_t interactions_count = 0, interaction_candidates = 0;
size_t i1, i2;
timestamp_t t0 = get_timestamp();
///*
for(std::pair<Point<3,float>, size_t>& p1 : arcl) {
i1 = p1.second;
auto iter_inner = arcl.getNNIterator<FAST>(p1.first); //full interactions
while (iter_inner.isNext()) {
auto p2 = iter_inner.get();
i2 = p2.second;
++interaction_candidates;
if(i1 != i2 && is_in_radius(p1.first, p2.first)) {
++interactions_count;
auto it = allInteractions.find(std::make_pair(i1, i2));
BOOST_REQUIRE(it != allInteractions.end());
allInteractions.erase(it);
//std::cout << i1 << " and " << i2 << " with p2 = " << p2.first.toString() << std::endl;
} //else std::cout << i1 << " | | " << i2 << " with p2 = " << p2.first.toString() << std::endl;