// ===================================================================================
// Copyright ScalFmm 2016 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.
// An extension to the license is given to allow static linking of scalfmm
// inside a proprietary application (no matter its license).
// See the main license file for more details.
//
// 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 <cstdio>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <limits>
#include "Utils/FParameters.hpp"
#include "Containers/FOctree.hpp"
#include "Components/FBasicCell.hpp"
#include "Components/FSimpleLeaf.hpp"
#include "Components/FBasicParticleContainer.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include "Utils/FMath.hpp"
#include "Files/FFmaGenericLoader.hpp"
#include "Utils/FParameterNames.hpp"
/**
* @brief Example of ScallFFM's use.
* @author B. Bramas, O. Coulaud, Q. Khan
*
* This example presents the basics of tree traversal with ScallFFM. Reading the
* source you can find out how to load a tree from an FMA formated file and how
* to move around the tree. Some statistics are produced to show how to get the
* information stored.
*
*/
int main(int argc, char ** argv)
{
typedef double FReal;
// Typedefs to shorten code
typedef FBasicCell CellClass;
typedef FBasicParticleContainer<FReal, 0, FReal > ContainerClass;
typedef FSimpleLeaf<FReal, ContainerClass > LeafClass;
typedef FOctree<FReal, CellClass, ContainerClass, LeafClass > OctreeClass;
namespace FPD = FParameterDefinitions;
// Constant strings for output
const std::string lhead("[STAT] ");
const std::string newline = "\n" + lhead;
// Parameter handling
// -------------------------------------------------------------------------
// Custom parameter
const FParameterNames LocalOptionHist = {
{"-histP", "--histogram-stat"}, // option enabling flags
"Build the histogram of the particle number per leaf." // option description
};
// ScalFFM generic options processing + program's arguments description
FHelpDescribeAndExit(argc, argv,
"Driver to obtain statistics on the octree.",
FPD::InputFile,
FPD::OctreeHeight,
FPD::OctreeSubHeight,
FPD::OutputFile,
LocalOptionHist);
// Octree, input and output parameters
const int TreeHeight = FParameters::getValue(argc,argv,FPD::OctreeHeight.options, 4);
const int SubTreeHeight = FParameters::getValue(argc,argv,FPD::OctreeSubHeight.options, 1);
const std::string inFileName = FParameters::getStr(argc, argv, FPD::InputFile.options, "../Data/test20k.fma");
const std::string outFileName =
FParameters::getStr(argc, argv, FPD::OutputFile.options, "output.txt");
std::cout << "Parameters " << std::endl
<< "\tOctree Depth : " << TreeHeight << std::endl
<< "\tSubOctree depth: " << SubTreeHeight << std::endl
<< "\tInput file : " << inFileName << std::endl
<< "\tOutput file : " << outFileName << std::endl
<< std::endl;
// Creating and Inserting particles in the tree
// -------------------------------------------------------------------------
// The loader reads FMA formated files. The file header is read
// automatically : we have basic information on the tree
// structure. Initially, the tree is empty.
FFmaGenericLoader<FReal> loader(inFileName.c_str());
OctreeClass tree(TreeHeight, SubTreeHeight,loader.getBoxWidth(),loader.getCenterOfBox());
std::cout << newline
<< "Particle file box." << newline
<< "\t\tCentre: " << loader.getCenterOfBox() << newline
<< "\t\tLength: " << loader.getBoxWidth()
<< std::endl << std::endl;
FReal physicalValue;
FPoint<FReal> particlePosition,
minPos(loader.getBoxWidth(),loader.getBoxWidth(),loader.getBoxWidth()),
maxPos(0., 0., 0.);
// Insertion of particles in the tree.
for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
// read next particle in file
loader.fillParticle(&particlePosition,&physicalValue);
// insert particle into tree
tree.insert(particlePosition);
// Box bounds stats
minPos.setX( FMath::Min(minPos.getX(), particlePosition.getX()) ) ;
minPos.setY( FMath::Min(minPos.getY(), particlePosition.getY()) ) ;
minPos.setZ( FMath::Min(minPos.getZ(), particlePosition.getZ()) ) ;
maxPos.setX( FMath::Max(maxPos.getX(), particlePosition.getX()) ) ;
maxPos.setY( FMath::Max(maxPos.getY(), particlePosition.getY()) ) ;
maxPos.setZ( FMath::Max(maxPos.getZ(), particlePosition.getZ()) ) ;
}
std::cout << lhead + "Particle bounding box " << newline
<< "\t\tMin corner: " << minPos << newline
<< "\t\tMax corner: " << maxPos << std::endl << std::endl;
// Statistics computation
// -------------------------------------------------------------------------
{ // Moving around the leaf level to get information on particles
long int allLeaves = (1 << (3 * (TreeHeight-1) )); // maximum number of leaves
FReal averageParticles = 0.0, varianceParticles = 0.0;
FSize nbLeafs = 0, nbPart = 0, nbTPart = 0; // number of particles, total number of particles
FSize minParticles = std::numeric_limits<FSize>::max(),
maxParticles = 0;
// Compute statistics on particles
// An octree iterator is used to get through the tree.
OctreeClass::Iterator octreeIterator(&tree);
// We move to the leftmost leaf.
octreeIterator.gotoBottomLeft();
do {
nbPart = octreeIterator.getCurrentListTargets()->getNbParticles() ;
minParticles = FMath::Min(minParticles,nbPart) ;
maxParticles = FMath::Max(maxParticles,nbPart) ;
nbTPart += nbPart;
varianceParticles += FReal(nbPart*nbPart) ;
++ nbLeafs;
} while(octreeIterator.moveRight()); // advancing through the level
averageParticles = FReal(nbTPart)/FReal(nbLeafs);
varianceParticles = varianceParticles/FReal(nbLeafs) - averageParticles*averageParticles;
std::cout.precision(4);
std::cout << lhead
<< "Leaf level : " << TreeHeight-1 << newline
<< "Max leaf count : " << allLeaves << newline
<< "Non empty leaves: " << nbLeafs
<< " (" << 100 * static_cast<FReal>(nbLeafs) / static_cast<FReal>(allLeaves)
<< "%)" << newline;
std::cout << "Particles in leaves:" << newline
<< "\t\tMin :" << std::setw(8) << minParticles
<< "\t\tMax :" << std::setw(8) << maxParticles << newline
<< "\t\tAverage :" << std::setw(8) << averageParticles
<< "\t\tVariance:" << std::setw(8) << varianceParticles
<< std::endl << std::endl;
// Histogram of particles per leaf
if( FParameters::existParameter(argc, argv, LocalOptionHist.options) ) {
std::vector<FSize> hist(maxParticles + 1, 0);
octreeIterator.gotoBottomLeft();
do {
nbPart = octreeIterator.getCurrentListSrc()->getNbParticles();
++hist[nbPart] ;
} while(octreeIterator.moveRight());
// write data
std::ofstream outfile(outFileName);
if( ! outfile.good() ) {
std::cerr << "Cannot open file " << outFileName << std::endl;
exit(EXIT_FAILURE);
}
outfile << "# Particle per leaf histogram. " << hist.size() << " chunk" << std::endl;
for(size_t i = 0 ; i < hist.size() ; ++i){
outfile << i << " " << hist[i] << std::endl;
}
}
// Statistics on particles done
FReal averageNeighbors = 0.0, varianceNeighbors = 0.0 ;
FSize nbBox, minBox = 100000, maxBox=0;
ContainerClass* neighborsP2P[27];
octreeIterator.gotoBottomLeft();
do {// P2P Neighbors
nbBox = tree.getLeafsNeighbors(neighborsP2P,
octreeIterator.getCurrentGlobalCoordinate(),
TreeHeight-1) ;
// need the current particles and neighbors particles
minBox = FMath::Min(minBox,nbBox) ;
maxBox = FMath::Max(maxBox,nbBox) ;
averageNeighbors += FReal(nbBox);
varianceNeighbors += FReal(nbBox*nbBox) ;
} while(octreeIterator.moveRight());
averageNeighbors /= FReal(nbLeafs) ;
varianceNeighbors = varianceNeighbors/FReal(nbLeafs)-averageNeighbors*averageNeighbors;
std::cout << lhead
<< "P2P Neighbors for each leaf " << newline
<< "\t\tMin :" << std::setw(8) << minBox
<< "\t\tMax :" << std::setw(8) << maxBox << newline
<< "\t\tAverage :" << std::setw(8) << averageNeighbors
<< "\t\tVariance:" << std::setw(8) << varianceNeighbors
<< std::endl << std::endl;
} // End of leaves statistics
{ // Internal cells statistics
long long int totalCells = 0;
long long int totalM2L = 0;
long long int totalM2ML2L = 0;
int nbCellsAtTop = 0;
int nbCellsAtBottom = 0;
// As for the leaf level, we use an iterator got run through the tree.
OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
// For each level, see end of loop to climb a level in the tree.
for(int idxLevel = TreeHeight - 1 ; idxLevel >= 1 ; --idxLevel){
int nbCellsAtLevel = 0;
int nbChildrenAtLevel = 0, adaptiveCell = 0, nbChildrenForMyCell;
int nbNeighborsAtLevel = 0;
int nbM2LNeighbors = 0, minM2L = 500, maxM2L = -1;
FReal averageM2LNeighbors = 0.0, varianceM2LNeighbors = 0.0;
const CellClass* neighborsM2L[343];
do {
++ nbCellsAtLevel;
// Count children
if( idxLevel != TreeHeight - 1 ) {
nbChildrenForMyCell = 0;
CellClass** children = octreeIterator.getCurrentChild();
for( int idxChild = 0 ; idxChild < 8 ; ++ idxChild ) {
// Non existing children are NULL
if( children[idxChild] )
++ nbChildrenForMyCell;
}
nbChildrenAtLevel += nbChildrenForMyCell;
if( nbChildrenForMyCell > 1 )
++ adaptiveCell;
}
// M2L Neighbors
nbM2LNeighbors =
tree.getInteractionNeighbors(neighborsM2L,
octreeIterator.getCurrentGlobalCoordinate(),
idxLevel);
nbNeighborsAtLevel += nbM2LNeighbors;
minM2L = FMath::Min(minM2L,nbM2LNeighbors) ;
maxM2L = FMath::Max(maxM2L,nbM2LNeighbors) ;
averageM2LNeighbors += FReal(nbM2LNeighbors) ;
varianceM2LNeighbors += FReal(nbM2LNeighbors*nbM2LNeighbors) ;
} while(octreeIterator.moveRight());
averageM2LNeighbors /= FReal(nbCellsAtLevel);
varianceM2LNeighbors = varianceM2LNeighbors / nbCellsAtLevel
- averageM2LNeighbors * averageM2LNeighbors;
totalCells += (long long int)(nbCellsAtLevel);
totalM2L += (long long int)(nbNeighborsAtLevel);
totalM2ML2L += (long long int)(nbChildrenAtLevel);
nbCellsAtTop = nbCellsAtLevel;
if( idxLevel == TreeHeight - 1 )
nbCellsAtBottom = nbCellsAtLevel;
std::cout << "[STAT] "
<< "Level " << idxLevel
<< newline
<< "\tCell count : " << nbCellsAtLevel
<< " (adaptive are cells with > 1 children)" << newline
<< "\t\tAdaptive:" << std::setw(8) << adaptiveCell
<< "\t\tNon adap:" << std::setw(8) << nbCellsAtLevel-adaptiveCell
<< newline
<< "\tM2M/L2L interact:" << std::setw(8) << nbChildrenAtLevel
<< "\t\tAverage :" << std::setw(8)
<< FReal(nbChildrenAtLevel)/FReal(nbCellsAtLevel)
<< newline
<< "\tM2L interactions:" << std::setw(8) << nbNeighborsAtLevel
<< newline;
std::cout << "\tM2L Neighbors / leaf: "
<< newline
<< "\t\tMin :" << std::setw(8) << minM2L
<< "\t\tMax :" << std::setw(8) << maxM2L
<< newline
<< "\t\tAverage :" << std::setw(8) << averageM2LNeighbors
<< "\t\tVariance:" << std::setw(8) << varianceM2LNeighbors
<< std::endl << std::endl;
// Go to next level
octreeIterator.moveUp();
// Move to leftmost cell on new level
octreeIterator.gotoLeft();
}
// Global statistics on the octree
std::cout
<< lhead + "Global stats" << newline
<< "\tCells = " << totalCells-nbCellsAtTop << newline
<< "\tM2L interactions / cell :" << std::setw(8)
<< totalM2L << newline
<< "\tAverage M2L inter / cell :" << std::setw(8)
<< FReal(totalM2L)/FReal(totalCells) << newline
<< "\tM2M/L2L interactions :" << std::setw(8)
<< totalM2ML2L << newline
<< "\tAverage M2M/L2L interact :" << std::setw(8)
<< FReal(totalM2ML2L-nbCellsAtTop) / FReal(totalCells-nbCellsAtBottom)
<< std::endl;
}
return 0;
}