TARGETS += switchchain_canonical_mixingtime

#include "exports.hpp"

#include "graph.hpp"

#include "graph_powerlaw.hpp"

#include "graph_spectrum.hpp"

#include "switchchain.hpp"

#include <algorithm>

#include <fstream>

#include <iostream>

#include <numeric>

#include <random>

#include <vector>

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

    // Simulation parameters

    const int numVerticesMin = 200;

    const int numVerticesMax = 1000;

    const int numVerticesStep = 200;

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

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

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

    };

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

        (void)n;

        (void)tau;

        return 10000;

    };

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

        (void)tau;

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

    };

    // Output file

    std::ofstream outfile;

    if (argc >= 2)

        outfile.open(argv[1]);

    else

        outfile.open("graphdata_canonical_mixingtime.m");

    if (!outfile.is_open()) {

        std::cout << "ERROR: Could not open output file.\n";

        return 1;

    }

    // Output Mathematica-style comment to indicate file contents

    outfile << "(*\n";

    outfile << "n from " << numVerticesMin << " to " << numVerticesMax

            << " step " << numVerticesStep << std::endl;

    outfile << "tauValues: " << tauValues << std::endl;

    outfile << "Canonical degree sequence.\n";

    outfile << "max time: 20 n\n";

    outfile << "samples per time: 1000\n";

    outfile << "time skip: n\n";

    outfile << "For uniform samples:\n";

    outfile << "mixingTime: 50 * (50 - 5 (tau - 2)) n\n";

    outfile << "measurements: 10000\n";

    outfile << "measureSkip: 30 n\n";

    outfile << "data:\n";

    outfile << "1: {n,tau}\n";

    outfile << "2: { {timestamp 1, {samples}}, {timestamp 2, {samples}} }\n";

    outfile << "3: {uniform samples}}\n";

    outfile << "*)" << std::endl;

    // Mathematica does not accept normal scientific notation

    outfile << std::fixed;

    outfile << '{' << '\n';

    bool outputComma = false;

    SwitchChain chain;

    Graph g;

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

         numVertices += numVerticesStep) {

        for (float tau : tauValues) {

            DegreeSequence ds;

            generateCanonicalPowerlawGraph(numVertices, tau, g, ds);

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

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

            if (outputComma)

                outfile << ',' << '\n';

            outputComma = true;

            outfile << '{' << '{' << numVertices << ',' << tau << '}';

#if 0

            std::vector<int> samples;

            outfile << ',' << '{';

            for (int maxTime = numVertices; maxTime <= 20 * numVertices;

                 maxTime += numVertices) {

                samples.clear();

                for (int sample = 0; sample < 1000; ++sample) {

                    chain.initialize(g, true);

                    for (int i = 0; i < maxTime; ++i)

                        chain.doMove(true);

                    samples.push_back(chain.g.getTrackedTriangles());

                }

                if (maxTime != numVertices)

                    outfile << ',';

                outfile << '{' << maxTime << ',' << samples << '}';

                std::cout << "t=" << maxTime << ' ' << std::flush;

            }

            outfile << '}';

#else

            std::vector<std::pair<int, std::vector<int>>> samples;

            for (int maxTime = numVertices; maxTime <= 20 * numVertices;

                 maxTime += numVertices) {

                samples.push_back(make_pair(maxTime, std::vector<int>()));

            }

            for (int sample = 0; sample < 5000; ++sample) {

                chain.initialize(g, true);

                int curTime = 0;

                for (auto &piv : samples) {

                    for (; curTime < piv.first; ++curTime)

                        chain.doMove(true);

                    piv.second.push_back(chain.g.getTrackedTriangles());

                }

            }

            outfile << ',' << samples;

#endif

            std::cout << "\nTaking uniform samples." << std::flush;

            // Uniform samples

            chain.initialize(g);

            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);

            std::vector<int> usamples;

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

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

                    chain.doMove(true);

                usamples.push_back(chain.g.getTrackedTriangles());

            }

            outfile << ',' << usamples;

            outfile << '}' << std::flush;

            std::cout << std::endl;

        }

    }

    outfile << '\n' << '}';

    return 0;

}

(* ::Package:: *)

gsraw=Import[NotebookDirectory[]<>"data/graphdata_canonical_mixingtime.m"];

gdata=GatherBy[gsraw,{#[[1,2]]&}];

(* Data format: *)

(* gdata[[ tau index, n index, datatype index ]] *)

(* datatype index:

1: {n,tau}

2: { {time1, {samples1}}, {time2, {samples2}} , ... }

3: {uniform samples}

*)

choices={1,2,3,6,8,12,14};

(*Histogram[histogramdata[[All,2]],Automatic,"Probability",ChartLegends\[Rule]histogramdata[[All,1]]]*)

makeHistogram[run_,choices_]:=Module[{gridsize,labels},

gridsize=Max[0.5,Mean[run[[3]]]/200];

labels="t = "<>ToString[#/run[[1,1]]]<>"\[CenterDot]n"&/@run[[2,choices,1]];

Show[

SmoothHistogram[run[[2,choices,2]],gridsize,PlotRange->{All,Automatic},PlotLegends->labels,PlotLabel->"n = "<>ToString[run[[1,1]]],ImageSize->500],

SmoothHistogram[run[[3]],gridsize,PlotLegends->{"uniform"},PlotStyle->{Thick,Black}]]]

makeHistogram[gdata[[1,-1]],choices]

makeHistogram[gdata[[-1,-1]],choices]

(* Get some sort of total variation distance between to sets of samples *)

getTVDistance[samples1_,samples2_]:=Module[{max,probs1,probs2},

max=Max[Max[samples1],Max[samples2]];

probs1=BinCounts[samples1,{0,max+1,1}];

probs1=probs1/Total[probs1];

probs2=BinCounts[samples2,{0,max+1,1}];

probs2=probs2/Total[probs2];

Total[Abs[probs1-probs2]]/2

]

getRatios[run_]:=Module[{avg,sd,scalefactor},

avg=Mean[run[[3]]];

sd=StandardDeviation[run[[3]]];

scalefactor=1/(run[[1,1]]*Log[run[[1,1]]]^0);

{"n,tau = "<>ToString[run[[1]]],

Map[{#[[1]]*scalefactor,Mean[#[[2]]]/avg}&,run[[2]]],

Map[{#[[1]]*scalefactor,(Mean[#[[2]]]-avg)/sd}&,run[[2]]],

Map[{#[[1]]*scalefactor,getTVDistance[#[[2]],run[[3]]]}&,run[[2]]]

}

]

ratios=Map[getRatios,gdata,{2}];

Map[ListPlot[#[[All,2]],Joined->True,PlotMarkers->Automatic,PlotRange->(1+{-0.15,+0.15}),PlotLegends->#[[All,1]],ImageSize->300]&,ratios]

Map[ListPlot[#[[All,3]],Joined->True,PlotMarkers->Automatic,PlotRange->(0+{-1,+1}),PlotLegends->#[[All,1]],ImageSize->300]&,ratios]

Map[ListPlot[#[[All,4]],Joined->True,PlotMarkers->Automatic,PlotRange->({0,1}),PlotLegends->#[[All,1]],ImageSize->300]&,ratios]