statisticsOnOctree.cpp 13.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
// ===================================================================================
// 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 <cstdio>
#include <cstdlib>
#include <iostream>
#include <fstream>

#include "../../Src/Utils/FParameters.hpp"
#include "../../Src/Containers/FOctree.hpp"

#include "../../Src/Components/FBasicCell.hpp"
#include "../../Src/Components/FSimpleLeaf.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
26
#include "../../Src/Components/FBasicParticleContainer.hpp"
27 28 29


#include "../../Src/Utils/FMath.hpp"
30
#include "../../Src/Files/FFmaGenericLoader.hpp"
31 32


COULAUD Olivier's avatar
COULAUD Olivier committed
33
/// \file  statisticsOnOctree.cpp
34
//!
COULAUD Olivier's avatar
COULAUD Olivier committed
35 36
//! \brief Driver to generate statistics on the octree for the particle distribution given by parameter -infile
//!  \authors B. Bramas, O. Coulaud
37
//!
COULAUD Olivier's avatar
COULAUD Olivier committed
38 39 40 41 42 43 44 45 46 47 48
//!  This code gives you some statistics (Particles, Cells, P2P and M2L number of operators) on the octree
//!    Those statistics are shown level by level.
//!
//! The statistics are
//!
//!  <b>for particles:</b> min/max particles per leaf, the average number, the variance and the average number of P2P neighbors.
//!
//! <b>For each level in the octree</b>
//!
//!  \arg       The number of cells, of adaptive cells (cell with more one child)
//!  \arg          The average, min and max numbers of M2L operators and also its  variance.
49 50 51

//!
//!  <b> General arguments:</b>
COULAUD Olivier's avatar
COULAUD Olivier committed
52
//!     \param   -help(-h)      to see the parameters available in this driver
COULAUD Olivier's avatar
COULAUD Olivier committed
53 54
//!     \param   -depth          The depth of the octree
//!     \param   -subdepth     Specifies the size of the sub octree
COULAUD Olivier's avatar
COULAUD Olivier committed
55 56
//!
//!     \param   -infile name   Name of the particles file. The file have to be in our FMA format
57
//!     \param    -bin              if the input file in binary mode
COULAUD Olivier's avatar
COULAUD Olivier committed
58 59
//!     \param   -outfile name Generic name  for output file  (without extension)
//!
60
//!  <b> Statistics options:</b>
COULAUD Olivier's avatar
COULAUD Olivier committed
61
//!   \param -histP build a file to generate histogram of particles per leaf. The data are store in file given by -outfile arguments and  .txt extension
62 63

// Simply create particles and try the kernels
COULAUD Olivier's avatar
COULAUD Olivier committed
64 65 66 67
//
void usage() {
	std::cout << "Driver to obtain statistics on the octree" << std::endl;
	std::cout <<	 "Options  "<< std::endl
COULAUD Olivier's avatar
COULAUD Olivier committed
68 69 70 71
                   <<     "      -help              to see the parameters    " << std::endl
                   <<	     "      -depth            the depth of the octree   "<< std::endl
                   <<	     "      -subdepth      specifies the size of the sub octree   " << std::endl
                   <<     "      -infile name    specifies the name of the particle distribution" << std::endl
72
                   <<    "        -bin              if the input file in binary mode"<< std::endl
COULAUD Olivier's avatar
COULAUD Olivier committed
73
                   <<     "      -outfile name  specifies the file for the diagnostics" << std::endl
COULAUD Olivier's avatar
COULAUD Olivier committed
74
                   <<     "      -histP              build the histogram of the particle number per leaf"<<std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
75 76
}

77
int main(int argc, char ** argv){
COULAUD Olivier's avatar
COULAUD Olivier committed
78 79
	typedef FBasicParticleContainer<0>                                    ContainerClass;
	typedef FSimpleLeaf< ContainerClass >                              LeafClass;
80
	typedef FOctree< FBasicCell, ContainerClass , LeafClass >  OctreeClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
81 82 83 84 85 86 87
	if(FParameters::existParameter(argc, argv, "-h")||FParameters::existParameter(argc, argv, "-help")|| (argc < 3 )){
		usage() ;
		exit(-1);
	}
	//
	//   Octree parameters
	//
COULAUD Olivier's avatar
COULAUD Olivier committed
88 89
	const int TreeHeight        = FParameters::getValue(argc,argv,"-depth", 5);
	const int SubTreeHeight  = FParameters::getValue(argc,argv,"-subdepth", 3);
COULAUD Olivier's avatar
COULAUD Olivier committed
90 91 92
	//
	//  input and output  Files parameters
	//
93 94
	const char* const filename = FParameters::getStr(argc,argv,"-infile", "../Data/test20k.fma");
	const std::string genericFileName(FParameters::getStr(argc,argv,"-outfile",   "output"));
COULAUD Olivier's avatar
COULAUD Olivier committed
95 96 97 98 99 100
	//
	std::cout <<	 "Parameters  "<< std::endl
			<<     "      Octree Depth      "<< TreeHeight <<std::endl
			<<	  "      SubOctree depth " << SubTreeHeight <<std::endl
			<<std::endl;

COULAUD Olivier's avatar
COULAUD Olivier committed
101
	//
102
	std::cout << "Opening : " << filename << "\n";
103 104 105 106 107
	bool binaryMode = false;
	if(FParameters::existParameter(argc, argv, "-bin")){
		binaryMode = true;
	}
	FFmaGenericLoader loader(filename,binaryMode);
108 109 110 111 112
	if(!loader.isOpen()){
		std::cout << "Loader Error, " << filename << " is missing\n";
		return 1;
	}
	// -----------------------------------------------------
COULAUD Olivier's avatar
COULAUD Olivier committed
113
	OctreeClass tree(TreeHeight, SubTreeHeight,loader.getBoxWidth(),loader.getCenterOfBox());
114 115 116 117 118
	//
	// -----------------------------------------------------
	//     Creating and Inserting particles in the tree
	// -----------------------------------------------------
	//
COULAUD Olivier's avatar
COULAUD Olivier committed
119 120 121 122
	std::cout << "Tree box is cubic "<<std::endl
			<< "         Centre:   "<< loader.getCenterOfBox() <<std::endl
			<< "         Length:  "<< loader.getBoxWidth()       <<std::endl <<std::endl;
	//
123
	std::cout << "Creating and Inserting " << loader.getNumberOfParticles() << " particles ..." << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
124
	std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
125 126
	FPoint particlePosition, minPos, maxPos;
	FReal physicalValue;
127 128
	for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
		loader.fillParticle(&particlePosition,&physicalValue);
COULAUD Olivier's avatar
COULAUD Olivier committed
129 130 131 132 133 134 135 136
		//
		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.setX(FMath::Max(maxPos.getY(),particlePosition.getY())) ;
		maxPos.setX(FMath::Max(maxPos.getZ(),particlePosition.getZ())) ;
		//
137 138
		tree.insert(particlePosition );
	}
COULAUD Olivier's avatar
COULAUD Olivier committed
139 140 141
	std::cout << "Data are inside the box delimited by "<<std::endl
			<< "         Min corner:  "<< minPos<<std::endl
			<< "         Max corner:  "<< maxPos<<std::endl <<std::endl;
142 143 144 145 146 147 148
	//
	// -----------------------------------------------------
	//     Start statistics
	// -----------------------------------------------------
	//
	{ // get stats
		{    // get stats on the leaf level (Particles)
COULAUD Olivier's avatar
COULAUD Olivier committed
149 150
			long int allLeaves =  (1 << (3* (TreeHeight-1) )) ;
			std::cout << std::endl<< "[STAT] Leaf level "  << " is  " << TreeHeight -1<< std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
151
			std::cout << "[STAT] potentials leafs number is " << allLeaves<< std::endl;
152 153 154 155

			FReal averageParticles = 0.0, varianceParticles = 0.0 ;
			int nbLeafs = 0,minParticles = 1000000.0, maxParticles = 0.0 ;
			OctreeClass::Iterator octreeIterator(&tree);
COULAUD Olivier's avatar
COULAUD Olivier committed
156 157
			//
			// Start to compute statistics on particles
158 159 160 161 162 163 164 165
			//
			int  nbPart ,nbTPart=0;
			octreeIterator.gotoBottomLeft();
			do{
				nbPart                   = octreeIterator.getCurrentListTargets()->getNbParticles() ;
				minParticles          =  FMath::Min(minParticles,nbPart) ;
				maxParticles         =  FMath::Max(maxParticles,nbPart) ;
				nbTPart              += nbPart;
COULAUD Olivier's avatar
COULAUD Olivier committed
166
				varianceParticles += FReal(nbPart*nbPart) ;
167 168 169
				++nbLeafs;
			} while(octreeIterator.moveRight());
			averageParticles   = nbTPart/FReal(nbLeafs);
COULAUD Olivier's avatar
COULAUD Olivier committed
170
			varianceParticles  = varianceParticles/FReal(nbLeafs) - averageParticles*averageParticles;
171
			//
COULAUD Olivier's avatar
COULAUD Olivier committed
172 173 174 175 176 177 178
			std::cout.precision(4);
			std::cout << "[STAT] Non empty leafs: " << nbLeafs << " % of non empty leaves: "<<100*static_cast<FReal>(nbLeafs)/static_cast<FReal>(allLeaves)<<" %" << std::endl;
			std::cout << "[STAT]  Particles on leafs:"  << std::endl
					<<   "[STAT]           Min:         "<< minParticles << std::endl
					<<   "[STAT]           Max:        "<< maxParticles << std::endl
					<<   "[STAT]           Average:  "<< averageParticles << std::endl
					<<   "[STAT]           Variance: " << varianceParticles << std::endl;
179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194
			//
			//  Histogram of particles per leaf
			//

			if(FParameters::existParameter(argc, argv, "-histP")){
				int size = maxParticles+1;
				int * hist = new int [size] ;
				memset(hist,0,(size)*sizeof(int));
				octreeIterator.gotoBottomLeft();
				do{
					nbPart  = octreeIterator.getCurrentListTargets()->getNbParticles() ;
					++hist[nbPart] ;
				} while(octreeIterator.moveRight());
				//
				// write data
				//
195 196
				std::string name(genericFileName + ".txt");
				std::ofstream outfile( name.c_str(), std::ofstream::out);
197 198 199 200 201 202 203 204 205 206
				if(!outfile) {
					std::cout << "Cannot open file "<< std::endl;
					exit(-1)	;
				}	//
				outfile << "# Particle histogram.   "<< size << " chunk" <<std::endl;
				for(int i=0 ; i < size ; ++i){
					outfile << i << "  " << hist[i] <<std::endl;
				}
				delete [] hist ;
			}
COULAUD Olivier's avatar
COULAUD Olivier committed
207 208
			FReal averageNeighbors = 0.0, varianceNeighbors =0.0 ;
			int nbBox,minBox=30,maxBox=0;
209
			octreeIterator.gotoBottomLeft();
COULAUD Olivier's avatar
COULAUD Olivier committed
210
			ContainerClass*  neighborsP2P[27];
211
			do{
COULAUD Olivier's avatar
COULAUD Olivier committed
212 213 214
				//
				//  P2P Neighbors
				//
COULAUD Olivier's avatar
COULAUD Olivier committed
215
				nbBox = tree.getLeafsNeighbors(neighborsP2P, octreeIterator.getCurrentGlobalCoordinate(),TreeHeight-1) ;
216
				// need the current particles and neighbors particles
COULAUD Olivier's avatar
COULAUD Olivier committed
217 218 219 220
				minBox                      =  FMath::Min(minBox,nbBox) ;
				maxBox                     =  FMath::Max(maxBox,nbBox) ;
				averageNeighbors   += FReal(nbBox);
				varianceNeighbors  += FReal(nbBox*nbBox) ;
221
			} while(octreeIterator.moveRight());
COULAUD Olivier's avatar
COULAUD Olivier committed
222 223 224 225 226 227 228 229 230
			//
			averageNeighbors/=FReal(nbLeafs) ;
			varianceNeighbors = varianceNeighbors/nbLeafs-averageNeighbors*averageNeighbors;
			//
			std::cout << "[STAT]  P2P Neighbors for each leaf " << std::endl
					<< "[STAT]           Min:       " <<  minBox << std::endl
					<< "[STAT]           Max:      " <<  maxBox << std::endl
					<< "[STAT]           Average: " <<  averageNeighbors<< std::endl
					<< "[STAT]           Variance: " <<  varianceNeighbors << std::endl<< std::endl;
231
		}
COULAUD Olivier's avatar
COULAUD Olivier committed
232
		//
233
		//    ---------------  END LEAVES ---------------
COULAUD Olivier's avatar
COULAUD Olivier committed
234
		//
235 236 237 238 239 240 241 242 243 244 245
		{
			long long int totalCells = 0;
			long long int totalM2L = 0;
			long long int totalM2ML2L = 0;

			int nbCellsAtTop = 0;
			int nbCellsAtBottom = 0;

			OctreeClass::Iterator octreeIterator(&tree);
			octreeIterator.gotoBottomLeft();

COULAUD Olivier's avatar
COULAUD Olivier committed
246
			for(int idxLevel = TreeHeight - 1 ; idxLevel > 1 ; --idxLevel){
247 248

				int nbCellsAtLevel = 0;
COULAUD Olivier's avatar
COULAUD Olivier committed
249
				int nbChildAtLevel = 0, adaptiveCell=0 ,nbChildForMyCell;
250
				int nbNeighborsAtLevel = 0;
COULAUD Olivier's avatar
COULAUD Olivier committed
251 252 253 254 255
				//
				int nbM2LNeighbors, minM2L=500,maxM2L=-1;
				FReal averageM2LNeighbors=0.0, varianceM2LNeighbors=0.0	;
				//
				const FBasicCell* neighborsM2L[343];
256 257
				do{
					++nbCellsAtLevel;
COULAUD Olivier's avatar
COULAUD Olivier committed
258
					// Check number of
COULAUD Olivier's avatar
COULAUD Olivier committed
259
					if( idxLevel != TreeHeight - 1 ){
COULAUD Olivier's avatar
COULAUD Olivier committed
260
						nbChildForMyCell=0 ;
261 262
						FBasicCell** child = octreeIterator.getCurrentChild();
						for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
COULAUD Olivier's avatar
COULAUD Olivier committed
263
							if(child[idxChild]) ++nbChildForMyCell;
264
						}
COULAUD Olivier's avatar
COULAUD Olivier committed
265 266
						nbChildAtLevel += nbChildForMyCell ;
						if(nbChildForMyCell>1) ++adaptiveCell ;
267 268 269
					}
					const FBasicCell* neighbors[343];
					nbNeighborsAtLevel += tree.getInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(),idxLevel);
COULAUD Olivier's avatar
COULAUD Olivier committed
270 271 272 273 274 275 276 277
					//
					//  M2L Neighbors
					//
					nbM2LNeighbors = tree.getInteractionNeighbors(neighborsM2L, octreeIterator.getCurrentGlobalCoordinate(),idxLevel);
					minM2L                          =  FMath::Min(minM2L,nbM2LNeighbors) ;
					maxM2L                         =  FMath::Max(maxM2L,nbM2LNeighbors) ;
					averageM2LNeighbors  += FReal(nbM2LNeighbors) ;
					varianceM2LNeighbors += FReal(nbM2LNeighbors*nbM2LNeighbors) ;
278
				} while(octreeIterator.moveRight());
COULAUD Olivier's avatar
COULAUD Olivier committed
279 280 281 282 283 284 285 286 287 288 289 290 291 292 293
				//
				averageM2LNeighbors/=FReal(nbCellsAtLevel) ;
				varianceM2LNeighbors = varianceM2LNeighbors/nbCellsAtLevel-averageM2LNeighbors*averageM2LNeighbors;

				std::cout << "[STAT] Level = " << idxLevel << std::endl
				               << "[STAT]     >> Nb Cells =                                 \t " << nbCellsAtLevel << std::endl
				               << "[STAT]     >> Nb Adaptive Cells =                   \t" << adaptiveCell    << "  Non Adaptive (1 son): " <<100*FReal(nbCellsAtLevel-adaptiveCell)/nbCellsAtLevel<<std::endl
				               << "[STAT]     >> Nb M2M/L2L interactions =        \t" << nbChildAtLevel << std::endl
				               << "[STAT]     >> Average M2M/L2L interactions = \t" << FReal(nbChildAtLevel)/FReal(nbCellsAtLevel) << std::endl
				               << "[STAT]     >> Nb M2L interactions =                 \t" << nbNeighborsAtLevel << std::endl;
				std::cout << "[STAT]     >> M2L Neighbors for each leaf " << std::endl
						       << "[STAT]             >>  Min:        " <<  minM2L << std::endl
					          <<  "[STAT]             >> Max:       " <<  maxM2L << std::endl
						      <<  "[STAT]             >> Average: " <<  averageM2LNeighbors<< std::endl
						      <<  "[STAT]             >> Variance: " <<  varianceM2LNeighbors << std::endl<< std::endl;
294 295 296 297 298

				totalCells += (long long int)(nbCellsAtLevel);
				totalM2L += (long long int)(nbNeighborsAtLevel);
				totalM2ML2L += (long long int)(nbChildAtLevel);
				nbCellsAtTop = nbCellsAtLevel;
COULAUD Olivier's avatar
COULAUD Olivier committed
299
				if( idxLevel == TreeHeight - 1 ) nbCellsAtBottom = nbCellsAtLevel;
COULAUD Olivier's avatar
COULAUD Olivier committed
300 301 302 303 304 305
				std::cout << std::endl;
				//
				//  Go to next level
				octreeIterator.moveUp();
				octreeIterator.gotoLeft();
               //
306
			}
COULAUD Olivier's avatar
COULAUD Olivier committed
307 308 309
			//
			// Global statistics on the octree
			//
310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327
			std::cout << "[STAT] For all the tree\n";
			std::cout << "[STAT] >> Total Nb Cells = " << totalCells-nbCellsAtTop << "\n";
			std::cout << "[STAT] >> Total Nb M2M/L2L interactions = " << totalM2ML2L << "\n";
			std::cout << "[STAT] >> Total Average M2M/L2L interactions = " << FReal(totalM2ML2L-nbCellsAtTop)/FReal(totalCells-nbCellsAtBottom) << "\n";
			std::cout << "[STAT] >> Total Nb M2L interactions per cell = " << totalM2L << "\n";
			std::cout << "[STAT] >> Total Average M2L interactions per cell = " << FReal(totalM2L)/FReal(totalCells) << "\n";

		}
	}

	// -----------------------------------------------------


	return 0;
}