// ===================================================================================
// 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;
}