Commit 828203e7 authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille
Browse files

FTreeBuilder and two tests added. Efficient Insert is available

parent 576ead56
......@@ -65,6 +65,12 @@ protected:
void addParticleValue(const int /*insertPosition*/){
}
/** Ending call for pushing array of attributes */
template<int index>
void addParticleValueS(const int /*insertPosition*/,const int /*nbParticles*/){
}
/** Filling call for each attributes values */
template<int index, typename... Args>
void addParticleValue(const int insertPosition, const AttributeClass value, Args... args){
......@@ -76,6 +82,23 @@ protected:
addParticleValue<index+1>( insertPosition, args...);
}
/** Filling call for each attributes values
* add multiples attributes from arrays
*/
template<int index, typename... Args>
void addParticleValueS(const int insertPosition, const int nbParts, const AttributeClass* value, Args... args){
// Compile test to ensure indexing
static_assert(index < NbAttributesPerParticle, "Index to get attributes is out of scope.");
for(int idxPart = 0; idxPart<nbParts ; ++idxPart){
// insert the value
attributes[index][insertPosition+idxPart] = value[idxPart];
// Continue for reamining values
}
addParticleValueS<index+1>( insertPosition, nbParts, args...);
}
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
......@@ -205,20 +228,23 @@ public:
*/
template<typename... Args>
void pushArray(const FPoint * inParticlePosition, FSize numberOfParts, Args... args){
const int positionToInsert = nbParticles;
//Tests if enough space
increaseSizeIfNeeded(numberOfParts);
for(int idxPart = 0; idxPart<numberOfParts ; ++idxPart){
// insert particle data
positions[0][nbParticles] = inParticlePosition[idxPart].getX();
positions[1][nbParticles] = inParticlePosition[idxPart].getY();
positions[2][nbParticles] = inParticlePosition[idxPart].getZ();
// insert attribute data
addParticleValue<0>( nbParticles, args[idxPart]...);
nbParticles += 1;
positions[0][positionToInsert + idxPart] = inParticlePosition[idxPart].getX();
positions[1][positionToInsert + idxPart] = inParticlePosition[idxPart].getY();
positions[2][positionToInsert + idxPart] = inParticlePosition[idxPart].getZ();
}
// insert attribute data
addParticleValueS<0>( nbParticles, numberOfParts ,args...);
nbParticles += numberOfParts;
}
/**
* Push called bu FSimpleLeaf
* Should have a particle position fallowed by attributes
......
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FTREEBUILDER_H
#define FTREEBUILDER_H
#include <omp.h>
#include "../Utils/FLog.hpp"
#include "../Utils/FQuickSort.hpp"
#include "../Utils/FTic.hpp"
#include "../Utils/FAssert.hpp"
#include "../Containers/FOctree.hpp"
#include "ScalFmmConfig.h"
/**
* @author Cyrille Piacibello
* @class FTreeBuilder
* @brief
* Please read the license
*
* This class provides a way to insert efficiently large amount of
* particles inside a tree.
*
* This is a static class. It's useless to instance it. This class use
* the Threaded QuickSort or the output of FMpiTreeBuilder in order to
* sort the parts and insert them.
*
*/
template<class ParticleClass, class OctreeClass, class LeafClass>
class FTreeBuilder{
private:
/**
* This method has been taken from the octree class,
* it computes a tree coordinate (x or y or z) from real cartesian position
*/
static int GetTreeCoordinate(const FReal inRelativePosition, const FReal boxWidthAtLeafLevel,
const FReal boxWidth, const int height) {
FAssertLF( (inRelativePosition >= 0 && inRelativePosition <= boxWidth), "inRelativePosition : ",inRelativePosition );
if(inRelativePosition == boxWidth){
return FMath::pow2(height-1)-1;
}
const FReal indexFReal = inRelativePosition / boxWidthAtLeafLevel;
return static_cast<int>(indexFReal);
}
/**
* A particle may not have a MortonIndex Method (set/get morton
* index) But in this algorithm they are sorted based on their
* morton indexes. So an IndexedParticle is storing a real
* particle + its index.
*
*/
struct IndexedParticle{
public:
MortonIndex index;
ParticleClass particle;
operator MortonIndex() const {
return this->index;
}
bool operator<=(const IndexedParticle& rhs){
return this->index <= rhs.index;
}
};
public:
/**
*
* Fill the tree with the parts from the array.
* The Particles must provide a getPosition and a getPhysicalValue(), method.
*/
template<class ContainerClass>
static void BuildTreeFromArray(ContainerClass* arrayToBeInserted, const FSize numberOfParticle, OctreeClass * tree,
bool isAlreadySorted=false){
//If the parts are already sorted, no need to sort again
FLOG(FTic enumTimer);
FLOG(FTic leavesPtr);
FLOG(FTic leavesOffset);
FLOG(FTic insertTimer);
FLOG(FTic copyTimer);
FLOG(FTic sortTimer);
//General values needed
int NbLevels = tree->getHeight();
FPoint centerOfBox = tree->getBoxCenter();
FReal boxWidth = tree->getBoxWidth();
FReal boxWidthAtLeafLevel = boxWidth/FReal(1 << (NbLevels - 1));
FPoint boxCorner = centerOfBox - boxWidth/2;
//We need to sort it in order to insert efficiently
//First, copy datas into an array that will be sorted and
//set Morton index for each particle
//Temporary FTreeCoordinate
FTreeCoordinate host;
IndexedParticle * toBeSorted = new IndexedParticle[numberOfParticle];
//Things to keep
FPoint* posToBeIn = new FPoint[numberOfParticle];
FReal * phyToBeIn = new FReal[numberOfParticle];
FLOG(copyTimer.tic());
for(int idxParts=0; idxParts<numberOfParticle ; ++idxParts ){
toBeSorted[idxParts].particle = arrayToBeInserted->data()[idxParts];
host.setX( GetTreeCoordinate(arrayToBeInserted->data()[idxParts].getPosition().getX() - boxCorner.getX(),
boxWidthAtLeafLevel, boxWidth, NbLevels));
host.setY( GetTreeCoordinate(arrayToBeInserted->data()[idxParts].getPosition().getY() - boxCorner.getY(),
boxWidthAtLeafLevel, boxWidth, NbLevels ));
host.setZ( GetTreeCoordinate(arrayToBeInserted->data()[idxParts].getPosition().getZ() - boxCorner.getZ(),
boxWidthAtLeafLevel, boxWidth, NbLevels ));
toBeSorted[idxParts].index = host.getMortonIndex(NbLevels-1);
}
FLOG(copyTimer.tac());
FLOG(FLog::Controller<<"Time needed for copying "<< numberOfParticle<<" particles : "<<copyTimer.elapsed() << " secondes !\n");
if(!isAlreadySorted){
//Sort dat array
FLOG(sortTimer.tic());
FQuickSort<IndexedParticle,MortonIndex>::QsOmp(toBeSorted,numberOfParticle);
FLOG(sortTimer.tac());
FLOG(FLog::Controller << "Time needed for sorting the particles : "<< sortTimer.elapsed() << " secondes !\n");
}
//Enumerate the different leaves AND copy the positions
unsigned int numberOfLeaves = 1;
FLOG(enumTimer.tic());
//First Values are copied
posToBeIn[0] = toBeSorted[0].particle.getPosition();
phyToBeIn[0] = toBeSorted[0].particle.getPhysicalValue();
for(int idxParts = 1 ; idxParts < numberOfParticle ; ++idxParts){
//Copy
posToBeIn[idxParts] = toBeSorted[idxParts].particle.getPosition();
phyToBeIn[idxParts] = toBeSorted[idxParts].particle.getPhysicalValue();
if(toBeSorted[idxParts].index != toBeSorted[idxParts-1].index){
numberOfLeaves++;
}
}
FLOG(enumTimer.tac());
FLOG(FLog::Controller << "Time needed for enumerate the leaves : "<< enumTimer.elapsed() << " secondes !\n");
FLOG(FLog::Controller << "Found " << numberOfLeaves << " leaves differents. \n");
//Store the size of each leaves
int * arrayOfSizeNbLeaves = new int[numberOfLeaves];
memset(arrayOfSizeNbLeaves,0,sizeof(int)*(numberOfLeaves));
//Init
int indexInLeafArray = -1;
arrayOfSizeNbLeaves[0] = 1;
FLOG(leavesOffset.tic());
MortonIndex currIndex = -1;
for(int idxParts = 0 ; idxParts < numberOfParticle ; ++idxParts){
if(toBeSorted[idxParts].index == currIndex){
arrayOfSizeNbLeaves[indexInLeafArray]++;
}
else{
FAssertLF(FSize(indexInLeafArray)<numberOfLeaves,"Problem there : ",indexInLeafArray);
indexInLeafArray++;
//There is alway at least 1 part in each leaf
arrayOfSizeNbLeaves[indexInLeafArray] = 1;
currIndex = toBeSorted[idxParts].index;
}
}
//Debug
int acc = 0;
for(unsigned int i=0 ; i<numberOfLeaves ; ++i){
//printf("acc : %d arrayofsize[%d] = %d \n",acc,i,arrayOfSizeNbLeaves[i]);
acc += arrayOfSizeNbLeaves[i];
}
printf("Tot : %d/%lld\n",acc,numberOfParticle);
FLOG(leavesOffset.tac());
FLOG(FLog::Controller << "Time needed for setting the offset of each leaves : "<< leavesOffset.elapsed() << " secondes !\n");
//Then, we create the leaves inside the tree
//Idx of the first part in this leaf in the array of part
int idxOfFirstPartInLeaf = 0;
//struct to store leaves and idx
struct LeafToFill{
LeafClass * leaf;
FSize idxOfLeafInPartArray;
};
FLOG(leavesPtr.tic());
LeafToFill * leavesToFill = new LeafToFill[numberOfLeaves];
memset(leavesToFill,0,sizeof(struct LeafToFill)*numberOfLeaves);
for(FSize idxLeaf = 0; idxLeaf < numberOfLeaves ; ++idxLeaf){
leavesToFill[idxLeaf].leaf = tree->createLeaf(toBeSorted[idxOfFirstPartInLeaf].index);
leavesToFill[idxLeaf].idxOfLeafInPartArray = idxOfFirstPartInLeaf;
idxOfFirstPartInLeaf += arrayOfSizeNbLeaves[idxLeaf];
}
FLOG(leavesPtr.tac());
FLOG(FLog::Controller << "Time needed for creating empty leaves : "<< leavesPtr.elapsed() << " secondes !\n");
FLOG(insertTimer.tic());
//Copy each parts into corresponding Leaf
#pragma omp parallel for schedule(auto)
for(FSize idxLeaf=0 ; idxLeaf<numberOfLeaves ; ++idxLeaf ){
//Task consists in copy the parts inside the leaf
leavesToFill[idxLeaf].leaf->pushArray(&posToBeIn[leavesToFill[idxLeaf].idxOfLeafInPartArray],
arrayOfSizeNbLeaves[idxLeaf],
&phyToBeIn[leavesToFill[idxLeaf].idxOfLeafInPartArray]);
}
FLOG(insertTimer.tac());
FLOG(FLog::Controller << "Time needed for inserting the parts into the leaves : "<< insertTimer.elapsed() << " secondes !\n");
//Clean the mess
delete [] leavesToFill;
delete [] arrayOfSizeNbLeaves;
delete [] posToBeIn;
delete [] phyToBeIn;
delete [] toBeSorted;
}
};
#endif //FTREEBUILDER_H
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#include <iostream>
#include <cstdlib>
#include <string>
#include <stdexcept>
#include <algorithm>
#include <vector>
#include <cstdio>
#include "ScalFmmConfig.h"
#include "Files/FFmaGenericLoader.hpp"
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include "../../Src/Kernels/P2P/FP2PParticleContainer.hpp"
#include "../../Src/Kernels/Rotation/FRotationKernel.hpp"
#include "../../Src/Kernels/Rotation/FRotationCell.hpp"
#include "../Src/BalanceTree/FLeafBalance.hpp"
#include "../Src/Containers/FTreeCoordinate.hpp"
#include "Utils/FTic.hpp"
#include "Utils/FQuickSort.hpp"
#include "Utils/FParameters.hpp"
#include "../../Src/Utils/FParameterNames.hpp"
#include "Containers/FOctree.hpp"
#ifdef _OPENMP
#include "Core/FFmmAlgorithmThread.hpp"
#else
#include "Core/FFmmAlgorithm.hpp"
#endif
#include "Utils/FTemplate.hpp"
/** This method has been tacken from the octree
* it computes a tree coordinate (x or y or z) from real position
*/
static int getTreeCoordinate(const FReal inRelativePosition, const FReal boxWidthAtLeafLevel) {
const FReal indexFReal = inRelativePosition / boxWidthAtLeafLevel;
const int index = int(FMath::dfloor(indexFReal));
if( index && FMath::LookEqual(inRelativePosition, boxWidthAtLeafLevel * FReal(index) ) ){
return index - 1;
}
return index;
}
/**
* This program build a tree and insert the parts inside.
* Time needed for the insert is outputed
*
*/
int main(int argc, char** argv){
struct TestParticle{
MortonIndex index;
FSize indexInFile;
FPoint position;
FReal physicalValue;
const FPoint& getPosition()const{
return position;
}
TestParticle& operator=(const TestParticle& other){
index=other.index;
indexInFile=other.indexInFile;
position=other.position;
physicalValue=other.physicalValue;
return *this;
}
bool operator<=(const TestParticle& rhs)const{
if(rhs.index < this->index){return false;}
else{
if(rhs.index > this->index){return true;}
else{
if(rhs.indexInFile == this->indexInFile){
return true;
}
else {
return rhs.indexInFile> this->indexInFile ;
}
}
}
}
};
static const int P = 9;
typedef FRotationCell<P> CellClass;
typedef FP2PParticleContainer<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FOctree< CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FRotationKernel< CellClass, ContainerClass , P> KernelClass;
typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClassThread;
const int NbLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options, 5);
const int SizeSubLevels = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options, 3);
const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/test20k.fma");
const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, omp_get_max_threads());
omp_set_num_threads(NbThreads);
std::cout << "Using " << omp_get_max_threads() <<" threads" << std::endl;
std::cout << "Opening : " << filename << "\n";
FFmaGenericLoader loaderRef(filename);
if(!loaderRef.isOpen()){
std::cout << "LoaderRef Error, " << filename << " is missing\n";
return 1;
}
FTic regInsert;
// -----------------------------------------------------
{
OctreeClass treeRef(NbLevels, SizeSubLevels, loaderRef.getBoxWidth(), loaderRef.getCenterOfBox());
std::cout << "Creating & Inserting " << loaderRef.getNumberOfParticles() << " particles ..." << std::endl;
std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
regInsert.tic();
for(int idxPart = 0 ; idxPart < loaderRef.getNumberOfParticles() ; ++idxPart){
FPoint particlePosition;
FReal physicalValue;
loaderRef.fillParticle(&particlePosition,&physicalValue);
treeRef.insert(particlePosition, physicalValue );
}
regInsert.tac();
std::cout << "Time needed for regular insert : " << regInsert.elapsed() << " secondes" << std::endl;
}
//Second solution, parts must be sorted for that
//Timers :
FTic sortTimer;
FTic insertTimer;
FTic fillTimer;
FTic enumTimer;
FTic counter;
FTic leavesOffset;
FTic leavesPtr;
FFmaGenericLoader loader(filename);
if(!loader.isOpen()){
std::cout << "Loader Error, " << filename << " is missing\n";
return 1;
}
//Get the needed informations
FReal boxWidth = loader.getBoxWidth();
FReal boxWidthAtLeafLevel = boxWidth/FReal(1 << (NbLevels - 1));
FPoint centerOfBox = loader.getCenterOfBox();
FPoint boxCorner = centerOfBox - boxWidth/2;
FSize nbOfParticles = loader.getNumberOfParticles();
//Temporary TreeCoordinate
FTreeCoordinate host;
TestParticle * arrayOfParts = new TestParticle[nbOfParticles];
memset(arrayOfParts,0,sizeof(TestParticle)*nbOfParticles);
fillTimer.tic();
for(int idxPart = 0 ; idxPart < nbOfParticles ; ++idxPart){
loader.fillParticle(&arrayOfParts[idxPart].position,&arrayOfParts[idxPart].physicalValue);
//Build temporary TreeCoordinate
host.setX( getTreeCoordinate( arrayOfParts[idxPart].getPosition().getX() - boxCorner.getX(), boxWidthAtLeafLevel ));
host.setY( getTreeCoordinate( arrayOfParts[idxPart].getPosition().getY() - boxCorner.getY(), boxWidthAtLeafLevel ));
host.setZ( getTreeCoordinate( arrayOfParts[idxPart].getPosition().getZ() - boxCorner.getZ(), boxWidthAtLeafLevel ));
//Set Morton index from Tree Coordinate
arrayOfParts[idxPart].index = host.getMortonIndex(NbLevels - 1);
arrayOfParts[idxPart].indexInFile = idxPart;
}
fillTimer.tac();
std::cout << "Time needed for filling the array : "<< fillTimer.elapsed() << " secondes !" << std::endl;
sortTimer.tic();
//std::sort(arrayOfParts,&arrayOfParts[nbOfParticles-1]);
FQuickSort<TestParticle,MortonIndex>::QsOmp(arrayOfParts,nbOfParticles);
sortTimer.tac();
std::cout << "Time needed for sorting the array : "<< sortTimer.elapsed() << " secondes !" << std::endl;
//Enumerate the different leaves
unsigned int numberOfLeaves = 0;
enumTimer.tic();
for(int idxParts = 1 ; idxParts < nbOfParticles ; ++idxParts){
if(arrayOfParts[idxParts].index != arrayOfParts[idxParts-1].index){
numberOfLeaves++;
}
}
enumTimer.tac();
std::cout << "Time needed for enumerate the leaves : "<< enumTimer.elapsed() << " secondes !" << std::endl;
std::cout << "Found " << numberOfLeaves << " leaves differents." << std::endl;
//Store the size of each subOctree
int * arrayOfSizeNbLeaves = new int[numberOfLeaves];
memset(arrayOfSizeNbLeaves,0,sizeof(int)*(numberOfLeaves));
//Init
int indexInLeafArray = 0;
arrayOfSizeNbLeaves[0] = 1;
leavesOffset.tic();
for(int idxParts = 1 ; idxParts < nbOfParticles ; ++idxParts){
if(arrayOfParts[idxParts].index == arrayOfParts[idxParts-1].index){
arrayOfSizeNbLeaves[indexInLeafArray]++;
}
else{
indexInLeafArray++;
}
}
leavesOffset.tac();
std::cout << "Time needed for setting the offset of each leaves : "<< leavesOffset.elapsed() << " secondes !" << std::endl;
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
//Then, we create the leaves inside the tree
//Idx of the first part in this leaf in the array of part
int idxOfFirstPartInLeaf = 0;
struct LeafToFill{
LeafClass * leaf;
FSize idxOfLeafInPartArray;
};
leavesPtr.tic();
LeafToFill * leavesToFill = new LeafToFill[numberOfLeaves];
memset(leavesToFill,0,sizeof(struct LeafToFill)*numberOfLeaves);
insertTimer.tic();
for(FSize idxLeaf = 0; idxLeaf < numberOfLeaves ; ++idxLeaf){
leavesToFill[idxLeaf].leaf = tree.createLeaf(arrayOfParts[idxOfFirstPartInLeaf].index);
idxOfFirstPartInLeaf += arrayOfSizeNbLeaves[idxLeaf];
leavesToFill[idxLeaf].idxOfLeafInPartArray += idxOfFirstPartInLeaf;
}
leavesPtr.tac();
std::cout << "Time needed for creating empty leaves : "<< leavesPtr.elapsed() << " secondes !" << std::endl;
#pragma omp parallel for schedule(auto)
for(FSize idxLeaf=0 ; idxLeaf<numberOfLeaves ; ++idxLeaf ){
//Task consists in copy the parts inside the leaf
leavesToFill[idxLeaf].leaf->pushArray(&arrayOfParts[leavesToFill[idxLeaf].idxOfLeafInPartArray].position,
arrayOfSizeNbLeaves[idxLeaf],
&arrayOfParts[leavesToFill[idxLeaf].idxOfLeafInPartArray].physicalValue);
}