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

ARCL SYM scheme is also working!

parent c54aef88
......@@ -111,6 +111,9 @@ public:
inline size_t maxlevel() {return max_level;}
// 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();}
inline size_t size() {return all_eles.size();}
inline T getEdgeLengthForLevel(size_t level) {
......@@ -118,7 +121,11 @@ public:
}
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
return getEdgeLengthForLevel(level+1);
}
inline unsigned int getLevelOfRadius(T radius) {
return std::ceil(std::log2(D_m / radius)) - 1;
}
/*! \warning Checks for < low bound and >= (!) high bound
......@@ -131,10 +138,6 @@ public:
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
......@@ -203,7 +206,7 @@ public:
T minradius = all_eles.last().first.get(dim);
//std::cout << "Min. radius: " << minradius << std::endl;
max_level = std::ceil(std::log2(D_m / minradius)) - 1;
max_level = getLevelOfRadius(minradius);
// 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):
......@@ -277,7 +280,7 @@ public:
throw std::invalid_argument("Point with cutoff radius > D_m searched!");
// Alternatively we could just limit it to D_m, but... that might yield wrong results.
unsigned int target_level = std::ceil(std::log2(D_m / radius)) - 1;
unsigned int target_level = getLevelOfRadius(radius);
//std::cout << p.toString() << " is on level " << target_level << std::endl;
if(target_level == 0)
......@@ -391,21 +394,21 @@ public:
* \param cell cell id
*
*/
template<unsigned int impl> inline AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,FULL,impl> getNNIterator(Point<dim+1,T> p)
template<unsigned int impl> inline AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,FULL,impl> getNNIterator(std::pair<Point<dim+1,T>, size_t> p)
{
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(Point<dim+1,T> p)
template<unsigned int impl> inline AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,SYM,impl> getNNIteratorSym(std::pair<Point<dim+1,T>, size_t> p)
{
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(Point<dim+1,T> p)
template<unsigned int impl> inline AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,CRS,impl> getNNIteratorCross(std::pair<Point<dim+1,T>, size_t> p)
{
AdaptiveCellNNIterator<dim,AdaptiveCellList<dim,T,BALANCED,ele_container>,CRS,impl> cln(p, *this);
......
......@@ -12,7 +12,12 @@
#include <vector>
#include <iostream>
template<unsigned int dim, typename CellS, unsigned int NNc_size, unsigned int impl> class AdaptiveCellNNIterator
// TODO: copied from CellNNIterator (including the fixed +1 at SYM!), doesn't look like a very good idea to me
#define FULL openfpm::math::pow(3,dim)
#define SYM openfpm::math::pow(3,dim)/2+1
#define CRS openfpm::math::pow(2,dim)
template<unsigned int dim, typename CellS, unsigned int NNc_size_original, unsigned int impl> class AdaptiveCellNNIterator
{
CellS & cl;
......@@ -21,42 +26,55 @@ template<unsigned int dim, typename CellS, unsigned int NNc_size, unsigned int i
size_t currentlevel;
// With FULL it always stays full, but with SYM we have to change it to FULL after the first level!
unsigned int NNc_size;
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) {
// Also: we pass NNc as a parameter, to use FULL inside SYM. See code below.
inline void populateIndices(int current_dim, Point<dim+1, typename CellS::value_type>& p, typename CellS::value_type edgelength, unsigned int NNc) {
if (current_dim < 0) {
if(cl.isInside(p))
cellindices.push_back(cl.findCellIndex(p));
//std::cout << p.toString() << std::endl;
}
else {
populateIndices(current_dim + 1, p, edgelength);
else if(NNc == SYM) {
// This means that the first cell we start in is the cell of the origin particle itself!
populateIndices(current_dim - 1, p, edgelength, SYM);
p.get(current_dim) += edgelength;
populateIndices(current_dim - 1, p, edgelength, FULL);
p.get(current_dim) -= edgelength;
populateIndices(current_dim + 1, p, edgelength);
}
else { // Assumption: FULL
populateIndices(current_dim - 1, p, edgelength, NNc);
p.get(current_dim) -= edgelength;
populateIndices(current_dim - 1, p, edgelength, NNc);
p.get(current_dim) += edgelength+edgelength;
populateIndices(current_dim + 1, p, edgelength);
populateIndices(current_dim - 1, p, edgelength, NNc);
p.get(current_dim) -= edgelength;
}
}
inline void updateCellIterators() {
cellindices.resize(0);
cellindices.reserve(openfpm::math::pow(3,dim));
cellindices.reserve(NNc_size);
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);
populateIndices(dim - 1, center, edgelength, NNc_size);
/*
std::cout << " -- On level " << currentlevel << ": ";
for (auto i : cellindices)
std::cout << i << ", ";
std::cout << std::endl;
if(NNc_size == SYM) {
std::cout << " -- On level " << currentlevel << ": ";
for (auto i : cellindices)
std::cout << i << ", ";
std::cout << std::endl;
}
//*/
cell_iter = cellindices.begin();
......@@ -64,7 +82,7 @@ template<unsigned int dim, typename CellS, unsigned int NNc_size, unsigned int i
}
inline void updateEleIterators() {
//std::cout << "Request ele_iters for cell " << *cell_iter << std::endl;
//if(NNc_size == SYM) 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;
......@@ -77,20 +95,42 @@ public:
* Cell NN iterator
*
*/
AdaptiveCellNNIterator(Point<dim+1, typename CellS::value_type> point, CellS & cl)
:cl(cl), p(point), currentlevel(0)
AdaptiveCellNNIterator(std::pair<Point<dim+1, typename CellS::value_type>, size_t> pointpair, CellS & cl)
:cl(cl), p(pointpair.first), currentlevel(0), NNc_size(NNc_size_original)
{
//std::cout << "NN of " << point.toString() << std::endl;
if(NNc_size != FULL)
currentlevel = cl.getLevelOfRadius(p.get(dim));
//std::cout << "NN of " << pointpair.second << " at " << p.toString() << " starting on lv " << currentlevel << std::endl;
updateCellIterators();
updateEleIterators();
while(ele_iter == ele_iter_end && currentlevel <= cl.maxlevel() + 1) {
if(NNc_size == SYM) {
// We know this cell contains an element - the original particle.
assert(ele_iter != ele_iter_end);
//for(auto it = ele_iter; it != ele_iter_end; ++it) std::cout << it->second << " at " << it->first.toString() << std::endl;
// So progress to it to ensure uniqueness for SYM.
//std::cout << "> ";
while (ele_iter->second != pointpair.second) {
++ele_iter;
//std::cout << "+";
}
++ele_iter;
//std::cout << " <, now at " << ele_iter->second << std::endl;
}
// Even in SYM the original particle might be the only one existing in the cell!
while(ele_iter == ele_iter_end && currentlevel <= cl.maxlevel()) {
++cell_iter;
if(cell_iter == cell_iter_end) {
++currentlevel;
if(NNc_size_original == SYM)
NNc_size = FULL;
updateCellIterators();
}
updateEleIterators();
}
......@@ -104,7 +144,7 @@ public:
*/
bool isNext()
{
//std::cout << currentlevel << " < " << cl.maxlevel() << std::endl;
//if(NNc_size == SYM) std::cout << currentlevel << " <= " << cl.maxlevel() << std::endl;
return currentlevel <= cl.maxlevel();
}
......@@ -113,15 +153,19 @@ public:
*/
AdaptiveCellNNIterator & operator++()
{
// Condition: don't start in an empty cell, otherwise this `++` will trick the following `if`.
// Condition: don't start in an empty cell, otherwise this `++` will trick the following `while`.
++ele_iter;
while(ele_iter == ele_iter_end && currentlevel <= cl.maxlevel() + 1) {
while(ele_iter == ele_iter_end && currentlevel <= cl.maxlevel()) {
++cell_iter;
if(cell_iter == cell_iter_end) {
++currentlevel;
if(NNc_size_original == SYM)
NNc_size = FULL;
updateCellIterators();
}
updateEleIterators();
}
......
......@@ -9,12 +9,12 @@
#define ADAPTIVECELLLIST_TEST_HPP_
#include "AdaptiveCellList.hpp"
#include "../CellList/CellListFast.hpp"
#include "Grid/grid_sm.hpp"
#include <typeinfo>
#include <unordered_set>
// http://stackoverflow.com/questions/20590656/error-for-hash-function-of-pair-of-ints
struct pairhash {
public:
......@@ -24,6 +24,17 @@ public:
return (3 * std::hash<T>()(x.first)) ^ std::hash<U>()(x.second);
}
};
struct pointhash {
public:
template <unsigned int dim, typename T>
std::size_t operator()(const Point<dim, T> &p) const
{
size_t res = 0;
for(unsigned int d = 0; d < dim; ++d)
res ^= std::hash<T>()(p.get(d));
return res;
}
};
#include <sys/time.h>
typedef unsigned long long timestamp_t;
......@@ -48,8 +59,11 @@ std::vector<Point<dim, float>> getSamplePoints() {
std::string rline, xline, yline;
int i=0;
for(; std::getline(rsource, rline); )
{
++i;
std::getline(xsource, xline);
std::getline(ysource, yline);
......@@ -63,10 +77,16 @@ std::vector<Point<dim, float>> getSamplePoints() {
xin >> x;
yin >> y;
// TODO: we actually lower these values, otherwise things break. i'm not sure that this should be happening.
if(x == 1.0f)
x -= std::numeric_limits<float>::epsilon();
if(y == 1.0f)
y -= std::numeric_limits<float>::epsilon();
if (dim == 2)
result.push_back({x, y});
else
result.push_back({y,x, r}); // TODO: change this back ;)
result.push_back({x,y, r});
float size = 0.05f;
int cr,cg,cb, index, nx,ny;
......@@ -138,8 +158,220 @@ AdaptiveCellList<2, float> getSampleARCL() {
return arcl;
}
inline bool is_in_radius(const Point<3, float> &p1, const Point<3, float> &p2) {
Point<3, float> diff = p1 - p2;
//std::cout << p1.toString() << " - " << p2.toString() << " = " << diff.toString() << std::endl << std::flush;
//std::cout << p2.get(2) << std::endl << std::flush;
//std::cout << diff[0]*diff[0] + diff[1]*diff[1] <<std::flush<< " <= " << (p1[2] < p2[2] ? p1[2] : p2[2]) << std::endl << std::flush;
float minradius = (p1[2] < p2[2] ? p1[2] : p2[2]);
return diff[0]*diff[0] + diff[1]*diff[1] <= minradius*minradius;
}
void interactionsOfClassicCL(unsigned int NNc_size)
{
std::cout << "## Classic CL:" << std::endl;
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]);
//std::cout << "Please choose a gridsize < " << 1.0f / (*maxradiusiterator)[2] << ": " << gridsize << std::endl;
const unsigned int dim = 2;
SpaceBox<2,float> box({0.0f,0.0f},{1.0f,1.0f});
size_t div[2] = {gridsize,gridsize};
grid_sm<2,void> g_info(div);
Point<2,float> org({0.0,0.0});
const float pad = 1.0f;
CellList<2, float> cl(box,div,org,pad, 30000);
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;
timestamp_t t0 = get_timestamp();
for (int i=pad; i<div[0]+pad; i++) {
for (int j=pad; j<div[1]+pad; j++) {
cell = cl.getGrid().LinId(grid_key_dx<2>(i, j));
//std::cout << "getting iterator for cell " << cell << " with elems: " << nels << std::endl;
switch(NNc_size) {
// differences in both blocks:
// * first line (I tried to use first templates, then function pointers, but nothing seemed to work)
// * SYM sorts the interaction pair
case FULL:
{
CellNNIterator<dim,CellList<dim,float,FAST>,FULL,FAST> iter = cl.getNNIterator<FAST>(cell);
while (iter.isNext()) {
i2 = iter.get();
//std::cout << "Neighbor found: " << i2 << std::endl;
CellIterator<CellList<2, float>> iter_inner(cell, cl);
while (iter_inner.isNext()) {
i1 = iter_inner.get();
++interaction_candidates;
if(i1 != i2 && is_in_radius(samplepoints3[i1], samplepoints3[i2])) {
//std::cout << i1 << " and " << i2 << " with p2 = " << samplepoints3[i2].toString() << std::endl;
allInteractions.insert(std::make_pair(i1, i2));
++interactions_count;
}
++iter_inner;
}
//std::cout << std::endl;
++iter;
}
}
break;
case SYM:
{
auto iter = cl.getNNIteratorSym<FAST>(cell);
while (iter.isNext()) {
i2 = iter.get();
//std::cout << "Neighbor found: " << i2 << std::endl;
CellIterator<CellList<2, float>> iter_inner(cell, cl);
while (iter_inner.isNext()) {
i1 = iter_inner.get();
++interaction_candidates;
/*
* We have a bit of a problem here.
* The iterator will at one point in time consider all particles in the current cell to be neighbors.
* Since we want symmetric interactions only once, we would like to filter them.
* The thing is: if we use it like here, the iterator doesn't know which element we are looking at (since we consider all, obviously)
* But we also cannot implement some filtering here, since we don't know where the elements come from.
*
* The way it is implemented now uses the auxiliary interaction storage vector, but this cannot be a proper solution.
*/
if(i1 != i2 && is_in_radius(samplepoints3[i1], samplepoints3[i2])) {
std::pair<size_t, size_t> interaction = std::make_pair(i1 < i2 ? i1 : i2, i1 < i2 ? i2 : i1);
auto isitinthere = allInteractions.find(interaction);
if(isitinthere == allInteractions.end()) {
//std::cout << i1 << " and " << i2 << " with p2 = " << samplepoints3[i2].toString() << std::endl;
allInteractions.insert(interaction);
++interactions_count;
}
}
++iter_inner;
}
//std::cout << std::endl;
++iter;
}
}
break;
}
}
}
timestamp_t t1 = get_timestamp();
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;
}
void interactionsOfARCL(unsigned int NNc_size)
{
std::cout << "## AR-CL:" << std::endl;
timestamp_t c0 = get_timestamp();
const unsigned int dim = 2; // needed for the FULL/SYM macros
auto arcl = getSampleARCL();
timestamp_t c1 = get_timestamp();
std::cout << "Creation time (us): " << c1-c0 << std::endl;
BOOST_REQUIRE_EQUAL(arcl.size(), samplepoints3.size());
size_t interactions_count = 0, interaction_candidates = 0;
size_t i1, i2;
const auto allInteractions_all(allInteractions);
timestamp_t t0 = get_timestamp();
std::unordered_set<Point<3,float>, pointhash> unipoints;
for(std::pair<Point<3,float>, size_t>& p1 : arcl) {
BOOST_REQUIRE(std::find(samplepoints3.begin(), samplepoints3.end(), p1.first) != samplepoints3.end());
unipoints.insert(p1.first);
}
//BOOST_REQUIRE_EQUAL(unipoints.size(), arcl.size());
// This is actually not the case, and it's causing problems.
arcl.construct(); // again
for(std::pair<Point<3,float>, size_t>& p1 : arcl) {
i1 = p1.second;
//if(i1 != 547) continue;
switch(NNc_size) {
// differences in both blocks:
// * first line (I tried to use first templates, then function pointers, but nothing seemed to work)
// * SYM sorts the interaction pair
case FULL:
{
auto iter = arcl.getNNIterator<FAST>(p1); //FULL interactions
while (iter.isNext()) {
auto p2 = iter.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;
++iter;
}
//std::cout << std::endl;
}
break;
case SYM:
{
//std::cout << "NN of pt " << p1.second << " at " << p1.first.toString() << std::endl << std::endl;
auto iter = arcl.getNNIteratorSym<FAST>(p1); //SYM interactions
while (iter.isNext()) {
auto p2 = iter.get();
//std::cout << "got for crunching: pt " << p2.second << " at " << p2.first.toString() << std::endl;
i2 = p2.second;
++interaction_candidates;
if(i1 != i2 && is_in_radius(p1.first, p2.first)) {
++interactions_count;
std::pair<size_t, size_t> interaction = std::make_pair(i1 < i2 ? i1 : i2, i1 < i2 ? i2 : i1);
auto it = allInteractions.find(interaction);
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;
//std::cout << "~~~< just crunched " << p2.first.toString() << std::endl;
++iter;
//std::cout << "~~~> finished counting up." << std::endl;
}
//std::cout << "\n/\\\n\\/\n" << std::endl;
}
break;
}
}
BOOST_REQUIRE_EQUAL(allInteractions.size(), 0);
timestamp_t t1 = get_timestamp();
std::cout << "Interaction 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_SUITE( ARCL_and_ARCLvsCL )
/*
// It doesn't really belong here, but I didn't want to just delete it.
BOOST_AUTO_TEST_CASE( celllist_realloc)
{
std::cout << "2D Cell List:" << std::endl;
......@@ -190,23 +422,21 @@ BOOST_AUTO_TEST_CASE( celllist_realloc)
int sum4 = 0;
for (unsigned int i=0; i<div[0]+2; i++) for (unsigned int j=0; j<div[1]+2; j++) sum4 += cl4.getNelements(cl4.getGrid().LinId(grid_key_dx<2>(i, j)));
/*
for (int i=0; i<div[0]+2; i++) {
for (int j=0; j<div[1]+2; j++) {
::printf("%6d ",cl0.getNelements(cl0.getGrid().LinId(grid_key_dx<2>(i, j))));
}
std::cout << std::endl;
}
std::cout << std::endl;
for (int i=0; i<div[0]+2; i++) {
for (int j=0; j<div[1]+2; j++) {
::printf("%6d ",cl1.getNelements(cl1.getGrid().LinId(grid_key_dx<2>(i, j))));
}
std::cout << std::endl;
}
std::cout << std::endl;
*/
// for (int i=0; i<div[0]+2; i++) {
// for (int j=0; j<div[1]+2; j++) {
// ::printf("%6d ",cl0.getNelements(cl0.getGrid().LinId(grid_key_dx<2>(i, j))));
// }
// std::cout << std::endl;
// }
// std::cout << std::endl;
//
// for (int i=0; i<div[0]+2; i++) {
// for (int j=0; j<div[1]+2; j++) {
// ::printf("%6d ",cl1.getNelements(cl1.getGrid().LinId(grid_key_dx<2>(i, j))));
// }
// std::cout << std::endl;
// }
// std::cout << std::endl;
size_t s = samplepoints2.size();
......@@ -223,78 +453,7 @@ BOOST_AUTO_TEST_CASE( celllist_realloc)
std::cout << "Without realloc, overestimate (us): " << t4-t3 << std::endl;
}
inline bool is_in_radius(const Point<3, float> &p1, const Point<3, float> &p2) {
Point<3, float> diff = p1 - p2;
//std::cout << p1.toString() << " - " << p2.toString() << " = " << diff.toString() << std::endl << std::flush;
//std::cout << p2.get(2) << std::endl << std::flush;
//std::cout << diff[0]*diff[0] + diff[1]*diff[1] <<std::flush<< " <= " << (p1[2] < p2[2] ? p1[2] : p2[2]) << std::endl << std::flush;
float minradius = (p1[2] < p2[2] ? p1[2] : p2[2]);
return diff[0]*diff[0] + diff[1]*diff[1] <= minradius*minradius;
}
BOOST_AUTO_TEST_CASE( get_all_interactions_classiccelllist)
{
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]);
//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.