#include <config.h>

#include <cmath>
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>

/*
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
#include <dune/grid/uggrid/ugwrapper.hh>
#pragma GCC diagnostic pop
*/

// dune-common includes
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/timer.hh>

// dune-istl includes
#include <dune/istl/matrix.hh>
#include <dune/istl/bcrsmatrix.hh>

// dune-fufem includes
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>

// dune-grid includes
#include <dune/grid/uggrid.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>

// dune-solver includes
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/norms/twonorm.hh>

// dune-faultnetworks includes
#include <dune/faultnetworks/utils/debugutils.hh>
#include <dune/faultnetworks/utils/matrixreader.hh>
#include <dune/faultnetworks/utils/vectorreader.hh>
#include <dune/faultnetworks/utils/levelinterfacenetworkwriter.hh>

#include <dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh>
#include <dune/faultnetworks/solvers/osccgsolver.hh>

#include <dune/faultnetworks/assemblers/intersectionjumpmassassembler.hh>
#include <dune/faultnetworks/assemblers/osclocalassembler.hh>

#include <dune/faultnetworks/oscrhs.hh>
#include <dune/faultnetworks/oscdata.hh>
#include <dune/faultnetworks/oscdatahierarchy.hh>

#include <dune/faultnetworks/faultp1nodalbasis.hh>
#include <dune/faultnetworks/interfacenetwork.hh>
#include <dune/faultnetworks/levelinterfacenetwork.hh>
#include <dune/faultnetworks/levelinterfacenetworkproblem.hh>
#include <dune/faultnetworks/dgmgtransfer.hh>

#include <dune/faultnetworks/faultfactories/cantorconvergencefaultfactory.hh>

int main(int argc, char** argv) { try
{
    MPIHelper::instance(argc, argv);

    const int dim = 2;
    //typedef double ctype;


    typedef typename Dune::FieldVector<double, dim> WorldVectorType;
    typedef typename Dune::BlockVector<Dune::FieldVector<double,1>> VectorType;
    typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>> MatrixType;


    #if HAVE_UG
    typedef typename Dune::UGGrid<dim> GridType;
    #else
    #error No UG!
    #endif


    typedef typename GridType::LevelGridView GV;
    typedef FaultP1NodalBasis<GV, double> DGBasis;
    typedef typename DGBasis::LocalFiniteElement DGFE;

    typedef OscLocalAssembler<GridType, DGFE, DGFE> LocalOscAssembler;
    typedef IntersectionJumpMassAssembler<GV, DGFE, DGFE> LocalInterfaceAssembler;

  // parse parameter file
    ParameterTree parameterSet;
    if (argc==2)
        ParameterTreeParser::readINITree(argv[1], parameterSet);
    else
        ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/faultnetworks/src/cantorfaultnetworks/cantorconvergence.parset", parameterSet);

    const std::string path = parameterSet.get<std::string>("path");
    const std::string resultPath = parameterSet.get<std::string>("resultPath");

    int problemCount = 0;
    while (parameterSet.hasSub("problem" + std::to_string(problemCount))) {

    ParameterTree& problemParameters = parameterSet.sub("problem" + std::to_string(problemCount));

    const std::string oscDataFile = problemParameters.get<std::string>("oscDataFile");

    const int minResolution = problemParameters.get<int>("minResolution");
    const int maxResolution = problemParameters.get<int>("maxResolution");
    const int exactResolution = problemParameters.get<int>("exactResolution");
    const int minCantorResolution = problemParameters.get<int>("minCantorResolution");

    const double c = problemParameters.get<double>("penaltyFactor");

    const int maxGridN = std::pow(2, maxResolution);


    std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(minResolution) + "_" + std::to_string(maxResolution) + "_" + std::to_string(exactResolution) + "_" + std::to_string(minCantorResolution) +  "_" + std::to_string((int) c) + ".log");
    std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf
    std::cout.rdbuf(out.rdbuf()); //redirect std::cout to out.txt


    std::cout << std::endl;
    std::cout << "Parameter file read successfully!" << std::endl;

    typedef Dune::Matrix<Dune::FieldMatrix<double,1,1>, std::allocator<Dune::FieldMatrix<double,1,1>>> OscMatrixType;
    OscMatrixType matrix;

    // interface integral jump penalty, B:
    OscMatrixType B(2,2);
    B[0][0] = 1;
    B[1][1] = 1;
    B[0][1] = -1;
    B[1][0] = -1;

    MatrixReader<OscMatrixType> matrixReader(path + oscDataFile);
    matrixReader.read(matrix);

    const int oscGridN = matrix.N();
    const int oscGridLevelIdx = std::round(std::log2(oscGridN)) - minResolution;

    if (oscGridN>maxGridN)
        DUNE_THROW(Dune::Exception, "Provided oscData too large!");

    std::cout << "OscData file read successfully!" << std::endl << std::endl;

    const int maxLevelIdx = maxResolution - minResolution;
    const int exactLevelIdx = exactResolution - minResolution;
    std::vector<size_t> maxCantorLevels(exactLevelIdx+1);
    for (int i=0; i<=exactLevelIdx; i++) {
        maxCantorLevels[i] = i+minCantorResolution;
    }

    // build multilevel cantor fault network
    CantorConvergenceFaultFactory<GridType> faultFactory(minResolution, exactLevelIdx, minCantorResolution, maxCantorLevels);
    const std::vector<std::shared_ptr<InterfaceNetwork<GridType>>>& interfaceNetworks = faultFactory.interfaceNetworks();

    /*
    for (size_t i=0; i<interfaceNetworks[1]->size(); i++) {
        LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(i) + ".tikz");
        networkWriter.write(interfaceNetworks[1]->levelInterfaceNetwork(i), true);
    }*/

    std::cout << "Interface networks built successfully!" << std::endl << std::endl;

    const GridType& grid = faultFactory.grid();

    OscData<GridType, OscMatrixType> oscData(grid, oscGridLevelIdx);
    oscData.set(matrix);
    OscDataHierarchy<GridType, OscMatrixType> oscDataHierarchy(oscData);

    Timer totalTimer;
    totalTimer.reset();
    totalTimer.start();


    // ----------------------------------------------------------------
    // ---              compute direct solutions                    ---
    // ----------------------------------------------------------------
    std::cout << std::endl;
    std::cout << "Computing direct solutions!" << std::endl;

    std::vector<VectorType> solutions;
    solutions.resize(exactLevelIdx+1);

    std::vector<double> approxErrors;
    approxErrors.resize(maxLevelIdx+1);

    double discError = NAN;

    OscRhs<WorldVectorType, Dune::FieldVector<double,1>> f;
    L2FunctionalAssembler<GridType, DGFE> l2FunctionalAssembler(f);

    const LevelInterfaceNetwork<GV>& exactLevelInterfaceNetwork = interfaceNetworks.back()->levelInterfaceNetwork(exactLevelIdx);
    LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> exactGlobalAssembler(exactLevelInterfaceNetwork);
    LocalOscAssembler exactLocalAssembler(oscDataHierarchy[exactLevelIdx]->data(), oscDataHierarchy[exactLevelIdx]->mapper());

    std::vector<std::shared_ptr<LocalInterfaceAssembler>> exactLocalInterfaceAssemblers(exactLevelInterfaceNetwork.size());
    for (size_t j=0; j<exactLocalInterfaceAssemblers.size(); j++) {
        //const double pen = penaltyFactor;
        const int k = minCantorResolution + exactLevelInterfaceNetwork.getIntersectionsLevels().at(j);
        const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
        exactLocalInterfaceAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
    }

    exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler);
    exactGlobalAssembler.solve();

    solutions[exactLevelIdx] = exactGlobalAssembler.getSol();
    approxErrors[maxLevelIdx] = 0;

    const DGBasis& exactBasis = exactGlobalAssembler.basis();

    EnergyNorm<MatrixType, VectorType> energyNorm(exactGlobalAssembler.matrix());

    for (int i=maxLevelIdx; i>=0; i--) {
        const LevelInterfaceNetwork<GV>& levelInterfaceNetwork = interfaceNetworks[i]->levelInterfaceNetwork(maxLevelIdx);
        LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> globalAssembler(levelInterfaceNetwork);
        LocalOscAssembler localAssembler(oscDataHierarchy[maxLevelIdx]->data(), oscDataHierarchy[maxLevelIdx]->mapper());

        std::vector<std::shared_ptr<LocalInterfaceAssembler>> localInterfaceAssemblers(levelInterfaceNetwork.size());
        for (size_t j=0; j<localInterfaceAssemblers.size(); j++) {
            //const double pen = penaltyFactor;
            const int k = minCantorResolution + levelInterfaceNetwork.getIntersectionsLevels().at(j);
            const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1);
            localInterfaceAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen);
        }

        globalAssembler.assemble(localAssembler, localInterfaceAssemblers, l2FunctionalAssembler);
        globalAssembler.solve();

        const VectorType& sol = globalAssembler.getSol();
        const DGBasis& basis = globalAssembler.basis();

        DGMGTransfer<DGBasis> levelTransfer(basis, exactBasis);
        levelTransfer.prolong(sol, solutions[i]);

        if (i==maxLevelIdx) {
            VectorType discErrorV(solutions[exactLevelIdx]);
            discErrorV -= solutions[maxLevelIdx];
            //approxErrors[i] = approxError.two_norm();
            discError= energyNorm.operator ()(discErrorV);
        } else {
            VectorType approxError(solutions[maxLevelIdx]);
            approxError -= solutions[i];
            //approxErrors[i] = approxError.two_norm();
            approxErrors[i] = energyNorm.operator ()(approxError);
        }
    }

    // print results
    std::cout << std::endl << std::endl;
    std::cout << "cantorLevel  h     approxError"<< std::endl;
    std::cout << "------------------------------" << std::endl;
    for (size_t i=0; i<approxErrors.size(); i++) {
        std::cout << (minCantorResolution+i) << "            2^-" << maxResolution << "  " << approxErrors[i] << std::endl;
    }

    std::cout << std::endl << std::endl;
    std::cout << "exactResolution: 2^-" << exactResolution << std::endl;
    std::cout << "discError: " << discError << std::endl;

    std::cout << std::endl;
    std::cout << "Total time: " << totalTimer.elapsed() << " seconds." << std::endl;

    std::cout.rdbuf(coutbuf); //reset to standard output again

    std::cout << "Problem " << problemCount << " done!" << std::endl;
    problemCount++;
    }
} catch (Dune::Exception e) {

    std::cout << e << std::endl;

}
}