#include "exports.hpp"

#include "graph.hpp"

#include "graph_ccm.hpp"

#include "graph_powerlaw.hpp"

#include "switchchain.hpp"

#include <algorithm>

#include <fstream>

#include <iostream>

#include <numeric>

#include <random>

#include <vector>

int main(int argc, char* argv[]) {

    // Simulation parameters

    float tauValues[] = {2.1f, 2.2f, 2.3f, 2.4f, 2.5f, 2.6f, 2.7f, 2.8f, 2.9f};

    //float tauValues[] = {2.1f, 2.3f, 2.5f, 2.7f, 2.9f};

    const int totalDegreeSamples = 200;

    auto getMixingTime = [](int n, float tau) {

        return int(50.0f * (50.0f - 30.0f * (tau - 2.0f)) * n);

    };

    auto getMeasurements = [](int n, float tau) {

        (void)n;

        (void)tau;

        return 100;

    };

    auto getMeasureSkip = [](int n, float tau) {

        (void)tau;

        return 10 * n; // Take a sample every ... steps

    };

    // Output file

    std::ofstream outfile;

    if (argc >= 2)

        outfile.open(argv[1]);

    else

    // Mathematica does not accept normal scientific notation

    outfile << std::fixed;

    outfile << '{';

    bool outputComma = false;

    std::mt19937 rng(std::random_device{}());

    Graph g;

    for (int numVertices = numVerticesMin; numVertices <= numVerticesMax;

         numVertices += numVerticesStep) {

        for (float tau : tauValues) {

            // For a single n,tau take samples over several instances of

            // the degree distribution.

            for (int degreeSample = 0; degreeSample < totalDegreeSamples;

                 ++degreeSample) {

                DegreeSequence ds;

                generatePowerlawGraph(numVertices, tau, g, ds, rng);

                std::cout << "Running (n,tau) = (" << numVertices << ',' << tau

                          << "). " << std::flush;

                //

                // Test the GCM1 and GCM2 success rate

                //

                int successrate1 = 0;

                int successrate2 = 0;

                for (int i = 0; i < 100; ++i) {

                    Graph gtemp;

                    // Take new highest degree every time

                    if (constrainedConfigurationModel(ds, gtemp, rng, false)) {

                        ++successrate1;

                    }

                    // Finish all pairings of highest degree first

                    if (constrainedConfigurationModel(ds, gtemp, rng, true)) {

                        ++successrate2;

                    }

                }

                SwitchChain chain;

                if (!chain.initialize(g)) {

                    std::cerr << "Could not initialize Markov chain.\n";

                    return 1;

                }

                long long trianglesTotal = 0;

                std::cout << " Finished CCM generation." << std::flush;

                int mixingTime = getMixingTime(numVertices, tau);

                for (int i = 0; i < mixingTime; ++i) {

                    chain.doMove();

                }

                chain.g.getTrackedTriangles() = chain.g.countTriangles();

                int measurements = getMeasurements(numVertices, tau);

                int measureSkip = getMeasureSkip(numVertices, tau);

                for (int i = 0; i < measurements; ++i) {

                    for (int j = 0; j < measureSkip; ++j)

                        chain.doMove(true);

                    trianglesTotal += chain.g.getTrackedTriangles();

                }