guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of /home/lavaux/Dropbox/gitRoot/CosmoToolbox

Browse Source
This commit is contained in:
Guilhem Lavaux 2012-06-04 22:21:14 -04:00
commit c1f1cf5a38
10 changed files with 502 additions and 65 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
sample/CMakeLists.txt
View File

@ -22,6 +22,9 @@ target_link_libraries(testDelaunay ${tolink})
add_executable(testNewton testNewton.cpp)
target_link_libraries(testNewton ${tolink})
add_executable(testPool testPool.cpp)
target_link_libraries(testPool ${tolink})
if (HDF5_FOUND)
add_executable(testReadFlash testReadFlash.cpp)
target_link_libraries(testReadFlash ${tolink})
@ -29,4 +32,7 @@ endif (HDF5_FOUND)
add_executable(testEskow testEskow.cpp)
target_link_libraries(testEskow ${tolink})
target_link_libraries(testEskow ${tolink})
add_executable(testAlgo testAlgo.cpp)
target_link_libraries(testAlgo ${tolink})

20
sample/testAlgo.cpp Normal file
View File

@ -0,0 +1,20 @@
#include <iostream>
#include <stdio.h>
#include "algo.hpp"
using namespace CosmoTool;
using namespace std;
int main(int argc, char **argv)
{
cout << square(2) << endl;
cout << cube(2) << endl;
cout << square(2.1f) << endl;
cout << cube(2.1f) << endl;
cout << spower<2>(2.1f) << endl;
cout << spower<3>(2.1f) << endl;
cout << spower<4>(2.1f) << endl;
cout << spower<5>(2.1f) << endl;
return 0;
}

19
sample/testPool.cpp Normal file
View File

@ -0,0 +1,19 @@
#include "pool.hpp"
using namespace CosmoTool;
int main(int argc, char **argv)
{
MemoryPool<int> pool(1024);
int **j = new int *[3000];
for (int i = 0; i < 3000; i++)
{
j[i] = pool.alloc();
j[i][0] = i;
}
pool.free_all();
return 0;
}

42
sample/testkd2.cpp
View File

@ -5,14 +5,16 @@
#include <fstream>
#define __KD_TREE_NUMNODES
#include "mykdtree.hpp"
#include "kdtree_splitters.hpp"
#define NTRY 100
#define NTRY 10
#define ND 3
using namespace std;
using namespace CosmoTool;
typedef KDTree<ND,char> MyTree;
typedef KDTree<ND,char,ComputePrecision,KD_homogeneous_cell_splitter<ND, char> > MyTree;
//typedef KDTree<ND,char,ComputePrecision > MyTree;
typedef KDCell<ND,char> MyCell;
MyCell *findNearest(MyTree::coords& xc, MyCell *cells, uint32_t Ncells)
@ -38,7 +40,7 @@ MyCell *findNearest(MyTree::coords& xc, MyCell *cells, uint32_t Ncells)
int main()
{
uint32_t Ncells = 100000;
uint32_t Ncells = 10000000;
MyCell *cells = new MyCell[Ncells];
for (int i = 0; i < Ncells; i++)
@ -48,8 +50,16 @@ int main()
cells[i].coord[l] = drand48();
}
// Check timing
clock_t startTimer = clock();
MyTree tree(cells, Ncells);
clock_t endTimer = clock();
clock_t delta = endTimer-startTimer;
double myTime = delta*1.0/CLOCKS_PER_SEC * 1.0;
cout << "KDTree build = " << myTime << " s" << endl;
MyTree::coords *xc = new MyTree::coords[NTRY];
cout << "Generating seeds..." << endl;
@ -68,31 +78,39 @@ int main()
for (int k = 0; k < NTRY; k++) {
cout << "Seed = " << xc[k][0] << " " << xc[k][1] << " " << xc[k][2] << endl;
tree.getNearestNeighbours(xc[k], 12, ngb, distances);
int last = -1;
for (uint32_t i = 0; i < 12; i++)
{
if (ngb[i] == 0)
continue;
last = i;
double d2 = 0;
for (int l = 0; l < 3; l++)
d2 += ({double delta = xc[k][l] - ngb[i]->coord[l]; delta*delta;});
fngb << ngb[i]->coord[0] << " " << ngb[i]->coord[1] << " " << ngb[i]->coord[2] << " " << sqrt(d2) << endl;
}
fngb << endl << endl;
double farther_dist = distances[11];
fngb << endl << endl;
double farther_dist = distances[last];
for (uint32_t i = 0; i < Ncells; i++)
{
bool found = false;
// If the points is not in the list, it means it is farther than the farther point
for (int j =0; j < 12; j++)
// If the points is not in the list, it means it is farther than the farthest point
for (int j =0; j < 12; j++)
{
if (&cells[i] == ngb[j]) {
found = true;
break;
}
}
double dist_to_seed = 0;
for (int l = 0; l < 3; l++)
{ double delta = xc[k][l]-cells[i].coord[l];
dist_to_seed += delta*delta; }
double dist_to_seed = 0;
for (int l = 0; l < 3; l++)
{
double delta = xc[k][l]-cells[i].coord[l];
dist_to_seed += delta*delta;
}
if (!found)
{
if (dist_to_seed <= farther_dist)

59
src/algo.hpp Normal file
View File

@ -0,0 +1,59 @@
#ifndef __COSMOTOOL_ALGO_HPP
#define __COSMOTOOL_ALGO_HPP
#include "config.hpp"
namespace CosmoTool
{
template<typename T>
T cube(T a)
{
return a*a*a;
}
template<typename T>
T square(T a)
{
return a*a;
}
template<int N>
class SPowerBase
{
public:
template<typename T>
static T spower(T a)
{
if ((N%2)==0)
{
T b = SPowerBase<N/2>::spower(a);
return b*b;
}
T b = SPowerBase<(N-1)/2>::spower(a);
return a*b*b;
}
};
template<>
class SPowerBase<0>
{
public:
template<typename T>
static T spower(T a)
{
return T(1);
}
};
template<int N, typename T>
T spower(T a)
{
return SPowerBase<N>::spower(a);
}
};
#endif

2
src/dinterpolate.tcc
View File

@ -217,7 +217,7 @@ namespace CosmoTool {
kdc[i] = c[i];
// It may happen that we are unlucky and have to iterate to farther
// neighbors. It should happen, especially on the boundaries.
// neighbors. It is bound to happen, especially on the boundaries.
do
{
uint32_t i;

134
src/kdtree_splitters.hpp Normal file
View File

@ -0,0 +1,134 @@
#ifndef __KDTREE_SPLITTERS_HPP
#define __KDTREE_SPLITTERS_HPP
#include <algorithm>
namespace CosmoTool
{
template<int N, typename ValType, typename CType = ComputePrecision>
struct KD_homogeneous_cell_splitter
{
typedef typename KDDef<N,CType>::KDCoordinates coords;
typedef typename KDDef<N,CType>::CoordType ctype;
void check_splitting(KDCell<N,ValType,CType> **cells, uint32_t Ncells, int axis, uint32_t split_index, ctype midCoord)
{
ctype delta = std::numeric_limits<ctype>::max();
assert(split_index < Ncells);
assert(axis < N);
for (uint32_t i = 0; i < split_index; i++)
{
assert(cells[i]->coord[axis] <= midCoord);
delta = min(midCoord-cells[i]->coord[axis], delta);
}
for (uint32_t i = split_index+1; i < Ncells; i++)
{
assert(cells[i]->coord[axis] > midCoord);
delta = min(cells[i]->coord[axis]-midCoord, delta);
}
assert(delta >= 0);
assert (std::abs(cells[split_index]->coord[axis]-midCoord) <= delta);
}
void operator()(KDCell<N,ValType,CType> **cells, uint32_t Ncells, uint32_t& split_index, int axis, coords minBound, coords maxBound)
{
if (Ncells == 1)
{
split_index = 0;
return;
}
ctype midCoord = 0.5*(maxBound[axis]+minBound[axis]);
uint32_t below = 0, above = Ncells-1;
ctype delta_min = std::numeric_limits<ctype>::max();
uint32_t idx_min = std::numeric_limits<uint32_t>::max();
while (below < above)
{
ctype delta = cells[below]->coord[axis]-midCoord;
if (delta > 0)
{
if (delta < delta_min)
{
delta_min = delta;
idx_min = above;
}
std::swap(cells[below], cells[above--]);
}
else
{
if (-delta < delta_min)
{
delta_min = -delta;
idx_min = below;
}
below++;
}
}
// Last iteration
{
ctype delta = cells[below]->coord[axis]-midCoord;
if (delta > 0)
{
if (delta < delta_min)
{
delta_min = delta;
idx_min = above;
}
}
else
{
if (-delta < delta_min)
{
delta_min = -delta;
idx_min = above;
}
}
}
if (idx_min != above)
{
bool cond1 = cells[idx_min]->coord[axis] > midCoord;
bool cond2 = cells[above]->coord[axis] > midCoord;
if ((cond1 && cond2) || (!cond1 && !cond2))
{
split_index = above;
std::swap(cells[above], cells[idx_min]);
}
else if (cond2)
{
if (above >= 1)
{
split_index = above-1;
std::swap(cells[above-1], cells[idx_min]);
}
else
split_index = 0;
assert(split_index >= 0);
}
else
{
if (above+1 < Ncells)
{
split_index = above+1;
std::swap(cells[above+1], cells[idx_min]);
}
else
split_index = Ncells-1;
assert(split_index < Ncells);
}
}
else split_index = above;
// check_splitting(cells, Ncells, axis, split_index, midCoord);
}
};
};
#endif

11
src/mykdtree.hpp
View File

@ -59,6 +59,7 @@ namespace CosmoTool {
uint32_t currentRank;
uint32_t numCells;
};
template<int N, typename ValType, typename CType = ComputePrecision>
class RecursionMultipleInfo
@ -76,7 +77,13 @@ namespace CosmoTool {
}
};
template<int N, typename ValType, typename CType = ComputePrecision>
template<int N, typename ValType, typename CType = ComputePrecision>
struct KD_default_cell_splitter
{
void operator()(KDCell<N,ValType,CType> **cells, uint32_t Ncells, uint32_t& split_index, int axis, typename KDDef<N,CType>::KDCoordinates minBound, typename KDDef<N,CType>::KDCoordinates maxBound);
};
template<int N, typename ValType, typename CType = ComputePrecision, typename CellSplitter = KD_default_cell_splitter<N,ValType,CType> >
class KDTree
{
public:
@ -84,6 +91,8 @@ namespace CosmoTool {
typedef typename KDDef<N>::KDCoordinates coords;
typedef KDCell<N,ValType,CType> Cell;
typedef KDTreeNode<N,ValType,CType> Node;
CellSplitter splitter;
KDTree(Cell *cells, uint32_t Ncells);
~KDTree();

106
src/mykdtree.tcc
View File

@ -23,13 +23,13 @@ namespace CosmoTool {
int rank;
};
template<int N, typename ValType, typename CType>
KDTree<N,ValType,CType>::~KDTree()
template<int N, typename ValType, typename CType, typename CellSplitter>
KDTree<N,ValType,CType,CellSplitter>::~KDTree()
{
}
template<int N, typename ValType, typename CType>
KDTree<N,ValType,CType>::KDTree(Cell *cells, uint32_t Ncells)
template<int N, typename ValType, typename CType, typename CellSplitter>
KDTree<N,ValType,CType,CellSplitter>::KDTree(Cell *cells, uint32_t Ncells)
{
base_cell = cells;
@ -43,8 +43,8 @@ namespace CosmoTool {
optimize();
}
template<int N, typename ValType, typename CType>
void KDTree<N,ValType,CType>::optimize()
template<int N, typename ValType, typename CType, typename CellSplitter>
void KDTree<N,ValType,CType,CellSplitter>::optimize()
{
coords absoluteMin, absoluteMax;
@ -76,9 +76,9 @@ namespace CosmoTool {
std::cout << " done." << std::endl;
}
template<int N, typename ValType, typename CType>
uint32_t KDTree<N,ValType,CType>::getIntersection(const coords& x, CoordType r,
KDTree<N,ValType,CType>::Cell **cells,
template<int N, typename ValType, typename CType, typename CellSplitter>
uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r,
KDTree<N,ValType,CType,CellSplitter>::Cell **cells,
uint32_t numCells)
throw (NotEnoughCells)
{
@ -96,8 +96,8 @@ namespace CosmoTool {
return info.currentRank;
}
template<int N, typename ValType, typename CType>
uint32_t KDTree<N,ValType,CType>::getIntersection(const coords& x, CoordType r,
template<int N, typename ValType, typename CType, typename CellSplitter>
uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r,
Cell **cells,
CoordType *distances,
uint32_t numCells)
@ -117,8 +117,8 @@ namespace CosmoTool {
return info.currentRank;
}
template<int N, typename ValType, typename CType>
uint32_t KDTree<N,ValType,CType>::countCells(const coords& x, CoordType r)
template<int N, typename ValType, typename CType, typename CellSplitter>
uint32_t KDTree<N,ValType,CType,CellSplitter>::countCells(const coords& x, CoordType r)
{
RecursionInfoCells<N,ValType> info;
@ -134,9 +134,9 @@ namespace CosmoTool {
return info.currentRank;
}
template<int N, typename ValType, typename CType>
template<int N, typename ValType, typename CType, typename CellSplitter>
template<bool justCount>
void KDTree<N,ValType,CType>::recursiveIntersectionCells(RecursionInfoCells<N,ValType,CType>& info,
void KDTree<N,ValType,CType,CellSplitter>::recursiveIntersectionCells(RecursionInfoCells<N,ValType,CType>& info,
Node *node,
int level)
throw (NotEnoughCells)
@ -208,25 +208,31 @@ namespace CosmoTool {
}
template<int N, typename ValType, typename CType>
KDTreeNode<N,ValType,CType> *KDTree<N,ValType,CType>::buildTree(Cell **cell0,
uint32_t Ncells,
uint32_t depth,
coords minBound,
coords maxBound)
void KD_default_cell_splitter<N,ValType,CType>::operator()(KDCell<N,ValType,CType> **cells, uint32_t Ncells, uint32_t& split_index, int axis, typename KDDef<N,CType>::KDCoordinates minBound, typename KDDef<N,CType>::KDCoordinates maxBound)
{
CellCompare<N,ValType,CType> compare(axis);
std::sort(cells, cells+Ncells, compare);
split_index = Ncells/2;
}
template<int N, typename ValType, typename CType, typename CellSplitter>
KDTreeNode<N,ValType,CType> *KDTree<N,ValType,CType,CellSplitter>::buildTree(Cell **cell0,
uint32_t Ncells,
uint32_t depth,
coords minBound,
coords maxBound)
{
if (Ncells == 0)
return 0;
int axis = depth % N;
Node *node = &nodes[lastNode++];
uint32_t mid = Ncells/2;
uint32_t mid;
coords tmpBound;
// Isolate the environment
{
CellCompare<N,ValType,CType> compare(axis);
std::sort(cell0, cell0+Ncells, compare);
}
splitter(cell0, Ncells, mid, axis, minBound, maxBound);
node->value = *(cell0+mid);
memcpy(&node->minBound[0], &minBound[0], sizeof(coords));
@ -252,8 +258,8 @@ namespace CosmoTool {
return node;
}
template<int N, typename ValType, typename CType>
uint32_t KDTree<N,ValType,CType>::countActives() const
template<int N, typename ValType, typename CType, typename CellSplitter>
uint32_t KDTree<N,ValType,CType,CellSplitter>::countActives() const
{
uint32_t numActive = 0;
for (uint32_t i = 0; i < lastNode; i++)
@ -264,9 +270,9 @@ namespace CosmoTool {
return numActive;
}
template<int N, typename ValType, typename CType>
template<int N, typename ValType, typename CType, typename CellSplitter>
typename KDDef<N,CType>::CoordType
KDTree<N,ValType,CType>::computeDistance(const Cell *cell, const coords& x) const
KDTree<N,ValType,CType,CellSplitter>::computeDistance(const Cell *cell, const coords& x) const
{
CoordType d2 = 0;
@ -278,9 +284,9 @@ namespace CosmoTool {
return d2;
}
template<int N, typename ValType, typename CType>
template<int N, typename ValType, typename CType, typename CellSplitter>
void
KDTree<N,ValType,CType>::recursiveNearest(
KDTree<N,ValType,CType,CellSplitter>::recursiveNearest(
Node *node,
int level,
const coords& x,
@ -352,9 +358,9 @@ namespace CosmoTool {
}
}
template<int N, typename ValType, typename CType>
template<int N, typename ValType, typename CType, typename CellSplitter>
KDCell<N,ValType,CType> *
KDTree<N,ValType,CType>::getNearestNeighbour(const coords& x)
KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbour(const coords& x)
{
CoordType R2 = INFINITY;
Cell *best = 0;
@ -364,9 +370,9 @@ namespace CosmoTool {
return best;
}
template<int N, typename ValType, typename CType>
template<int N, typename ValType, typename CType, typename CellSplitter>
void
KDTree<N,ValType,CType>::recursiveMultipleNearest(RecursionMultipleInfo<N,ValType,CType>& info, Node *node,
KDTree<N,ValType,CType,CellSplitter>::recursiveMultipleNearest(RecursionMultipleInfo<N,ValType,CType>& info, Node *node,
int level)
{
CoordType d2 = 0;
@ -384,11 +390,11 @@ namespace CosmoTool {
// If not it is in 1.
go = node->children[1];
other = node->children[0];
if (go == 0)
{
go = other;
other = 0;
}
// if (go == 0)
// {
// go = other;
//other = 0;
//}
}
if (go != 0)
@ -401,8 +407,8 @@ namespace CosmoTool {
computeDistance(node->value, info.x);
info.queue.push(node->value, thisR2);
info.traversed++;
if (go == 0)
return;
// if (go == 0)
// return;
// Now we found the best. We check whether the hypersphere
// intersect the hyperplane of the other branch
@ -421,8 +427,8 @@ namespace CosmoTool {
}
}
template<int N, typename ValType, typename CType>
void KDTree<N,ValType,CType>::getNearestNeighbours(const coords& x, uint32_t N2,
template<int N, typename ValType, typename CType, typename CellSplitter>
void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint32_t N2,
Cell **cells)
{
RecursionMultipleInfo<N,ValType> info(x, cells, N2);
@ -435,8 +441,8 @@ namespace CosmoTool {
// std::cout << "Traversed = " << info.traversed << std::endl;
}
template<int N, typename ValType, typename CType>
void KDTree<N,ValType,CType>::getNearestNeighbours(const coords& x, uint32_t N2,
template<int N, typename ValType, typename CType, typename CellSplitter>
void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint32_t N2,
Cell **cells,
CoordType *distances)
{
@ -470,8 +476,8 @@ namespace CosmoTool {
long rootId;
};
template<int N, typename ValType, typename CType>
void KDTree<N,ValType,CType>::saveTree(std::ostream& o) const
template<int N, typename ValType, typename CType, typename CellSplitter>
void KDTree<N,ValType,CType,CellSplitter>::saveTree(std::ostream& o) const
{
KDTreeHeader h;
@ -505,8 +511,8 @@ namespace CosmoTool {
}
}
template<int N, typename ValType, typename CType>
KDTree<N,ValType,CType>::KDTree(std::istream& in, Cell *cells, uint32_t Ncells)
template<int N, typename ValType, typename CType, typename CellSplitter>
KDTree<N,ValType,CType,CellSplitter>::KDTree(std::istream& in, Cell *cells, uint32_t Ncells)
throw (InvalidOnDiskKDTree)
{
KDTreeHeader h;

166
src/pool.hpp Normal file
View File

@ -0,0 +1,166 @@
#ifndef __COSMO_POOL_HPP
#define __COSMO_POOL_HPP
#include <list>
#include "config.hpp"
namespace CosmoTool
{
template<typename T>
struct PoolNode
{
T *data;
uint32_t last_free, size;
PoolNode<T> *next;
};
template<typename T> class MemoryPool;
template<typename T>
class MemoryIterator
{
private:
friend class MemoryPool<T>;
PoolNode<T> *cur, *previous;
uint32_t in_node;
MemoryIterator(PoolNode<T> *h)
{
cur = h;
previous = h;
in_node = 0;
}
public:
MemoryIterator() { cur = 0; }
~MemoryIterator() {}
const MemoryIterator& operator=(const MemoryIterator& i)
{
cur = i.cur;
previous = i.previous;
in_node = i.in_node;
}
bool operator==(const MemoryIterator& i) const
{
return (cur == i.cur) && (in_node == i.in_node);
}
MemoryIterator& operator++()
{
if (cur == 0)
return *this;
in_node++;
if (in_node == cur->size)
{
in_node = 0;
previous = cur;
cur = cur->next;
}
return *this;
}
T& operator*()
{
return cur->data[in_node];
}
T& operator->()
{
return cur->data[in_node];
}
};
// This is bare simple memory pools
template<typename T>
class MemoryPool
{
private:
uint32_t m_allocSize;
PoolNode<T> *head, *current;
typedef MemoryIterator<T> iterator;
public:
MemoryPool(uint32_t allocSize)
: m_allocSize(allocSize), head(0), current(0) {}
~MemoryPool()
{
free_all();
}
void free_all()
{
PoolNode<T> *node = head;
while (node != 0)
{
PoolNode<T> *next = node->next;
delete[] node->data;
delete node;
node = next;
}
current = head = 0;
}
T *alloc()
{
T *ret = alloc_in_node();
return (ret == 0) ? alloc_new_in_node() : ret;
}
iterator begin()
{
return iterator(head);
}
iterator end()
{
return iterator(0);
}
protected:
T *alloc_in_node()
{
if (current == 0 || current->last_free == current->size)
return 0;
return &current->data[current->last_free++];
}
T *alloc_new_in_node()
{
PoolNode<T> *newNode = new PoolNode<T>;
if (newNode == 0)
return 0;
newNode->last_free = 1;
newNode->size = m_allocSize;
newNode->data = new T[m_allocSize];
if (newNode->data == 0)
{
delete newNode;
return 0;
}
newNode->next = 0;
if (current == 0)
current = head = newNode;
else
{
current->next = newNode;
current = newNode;
}
return &newNode->data[0];
}
};
};
#endif