statisticsOnOctree.cpp 13.9 KB
Newer Older
1
// See LICENCE file at project root
2 3 4 5
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <fstream>
6
#include <limits>
7

8 9
#include "Utils/FParameters.hpp"
#include "Containers/FOctree.hpp"
10

11 12 13
#include "Components/FBasicCell.hpp"
#include "Components/FSimpleLeaf.hpp"
#include "Components/FBasicParticleContainer.hpp"
14
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
15

16 17
#include "Utils/FMath.hpp"
#include "Files/FFmaGenericLoader.hpp"
18

19
#include "Utils/FParameterNames.hpp"
20

21 22 23 24 25 26 27 28 29 30 31
/**
 * @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.
 * 
 */

COULAUD Olivier's avatar
COULAUD Olivier committed
32

33 34
int main(int argc, char ** argv) 
{
BRAMAS Berenger's avatar
BRAMAS Berenger committed
35 36 37 38 39 40 41 42 43 44 45 46 47
    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;

48 49 50
    //     Parameter handling
    // -------------------------------------------------------------------------
    // Custom parameter
COULAUD Olivier's avatar
COULAUD Olivier committed
51
    const FParameterNames LocalOptionHist = {
52 53
        {"-histP", "--histogram-stat"},                        // option enabling flags
        "Build the histogram of the particle number per leaf." // option description
54
    };
55 56

    // ScalFFM generic options processing + program's arguments description
57 58
    FHelpDescribeAndExit(argc, argv,
                         "Driver to obtain statistics on the octree.",
59 60 61 62 63 64 65
                         FPD::InputFile,
                         FPD::OctreeHeight,
                         FPD::OctreeSubHeight,
                         FPD::OutputFile, 
                         LocalOptionHist);

    // Octree, input and output parameters
BRAMAS Berenger's avatar
BRAMAS Berenger committed
66 67 68
    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");
69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
    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.
85
    FFmaGenericLoader<FReal> loader(inFileName.c_str());
86 87 88 89 90 91 92 93 94
    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;
95
    FPoint<FReal> particlePosition,
96 97 98 99
        minPos(loader.getBoxWidth(),loader.getBoxWidth(),loader.getBoxWidth()),
        maxPos(0., 0., 0.);

    // Insertion of particles in the tree.
100
    for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
        // 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;
125 126
        FSize nbLeafs = 0, nbPart = 0, nbTPart = 0; // number of particles, total number of particles
        FSize minParticles = std::numeric_limits<FSize>::max(),
127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
            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

143
        averageParticles  = FReal(nbTPart)/FReal(nbLeafs);
144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
        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) ) {
164
            std::vector<FSize> hist(maxParticles + 1, 0);
165 166 167 168 169 170 171 172 173 174 175 176 177 178

            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;
179
            for(size_t i = 0 ; i < hist.size() ; ++i){
180 181 182 183 184 185
                outfile << i << "  " << hist[i] << std::endl;
            }
        }
        // Statistics on particles done

        FReal averageNeighbors = 0.0, varianceNeighbors = 0.0 ;
186
        FSize nbBox, minBox = 100000, maxBox=0;
187 188 189 190 191 192 193 194 195 196 197 198 199 200 201
        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) ;
202
        varianceNeighbors  = varianceNeighbors/FReal(nbLeafs)-averageNeighbors*averageNeighbors;
203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322


        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;
323 324 325 326
}