#include <vector>

template <typename T>

std::ostream& operator<<(std::ostream& s, const std::vector<T>& v) {

    if (v.empty()) {

        s << '{' << '}';

        return s;

    }

    auto iter = v.begin();

    s << '{' << *iter++;

    for (; iter != v.end(); ++iter) {

        s << ',' << *iter;

    }

    s << '}';

    return s;

}

// Mathematica style export for arrays

// SFINAE to disable char arrays since it will mess up

// things of the form  s << "some string";

template <

    typename T, std::size_t N,

    typename = typename std::enable_if<!std::is_same<T, char>::value, T>::type>

std::ostream& operator<<(std::ostream& s, const T (&v)[N]) {

    if (N == 0) {

        s << '{' << '}';

        return s;

    }

    s << '{' << v[0];

    for (size_t i = 1; i < N; ++i)

        s << ',' << v[i];

    s << '}';

    return s;

}

#include <algorithm>

#pragma once

#include <cassert>

#include <cmath>

#include <random>

// Discrete powerlaw distribution

// obtained by rounding a continuous powerlaw distribution

class powerlaw_distribution {

  public:

    // xmin and xmax are inclusive

    powerlaw_distribution(float a, unsigned int _xmin, unsigned int _xmax)

        : alpha(a), xmin(_xmin), xmax(_xmax) {

        assert(xmin > 0);

        assert(xmin <= xmax);

        assert(alpha > 1.0f);

        _exponent = -1.0f / (alpha - 1.0f);

    }

    ~powerlaw_distribution() {}

    template <class Generator>

    unsigned int operator()(Generator& rng) const {

        unsigned int x = xmax + 1;

        while (x > xmax) {

            // Generate uniform value in [0,1)

            double r = std::generate_canonical<

                double, std::numeric_limits<double>::digits>(rng);

            x = std::round(std::pow(r, _exponent) * xmin);

        }

        return x;

    }

  private:

    float alpha;

    double _exponent;

    unsigned int xmin, xmax;

};

(* ::Subsection:: *)

(*Data import and data merge*)

gsraw=Import[NotebookDirectory[]<>"graphdata_merged.m"];

newData=Import[NotebookDirectory[]<>"graphdata_tau_multi3.m"];

mergedData=Import[NotebookDirectory[]<>"graphdata_merged.m"];

Export[NotebookDirectory[]<>"graphdata_merged_new.m",Join[mergedData,newData]]

(* ::Subsection:: *)

(*Plot empirical distribution of maximum degree*)

maxDegrees=Map[{#[[1]],Max[#[[3]]]}&,gsraw];

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

Histogram[maxDegrees[[1,-1,All,2]],PlotRange->{{0,1100},{0,200}},AxesLabel->{"d_max","frequency"}]

Histogram[maxDegrees[[2,-1,All,2]],PlotRange->{{0,1100},{0,200}},AxesLabel->{"d_max","frequency"}]

Histogram[maxDegrees[[3,-1,All,2]],PlotRange->{{0,1100},{0,200}},AxesLabel->{"d_max","frequency"}]

(* ::Subsection:: *)

(*Plot triangle count over "time" in Markov chain instances*)

numPlots=10;

minCount=Min[Map[Min[#[[2]]]&,gsraw[[1;;numPlots]]]];

maxCount=Max[Map[Max[#[[2]]]&,gsraw[[1;;numPlots]]]];

maxTime=Max[Map[Length[#[[2]]]&,gsraw[[1;;numPlots]]]];

skipPts=Round[maxTime/100]; (* Plotting every point is slow. Plot only once per `skipPts` timesteps *)

coarseData=Map[#[[2,1;;-1;;skipPts]]&,gsraw[[1;;numPlots]]];

labels=Map["{n,tau} = "<>ToString[#[[1]]]&,gsraw[[1;;numPlots]]];

ListPlot[coarseData,Joined->True,PlotRange->{minCount,maxCount},DataRange->{0,maxTime},PlotLegends->labels]

(* Map[ListPlot[#,Joined->True,PlotRange\[Rule]{minCount,maxCount},DataRange\[Rule]{0,maxTime}]&,coarseData] *)

(* ::Subsection:: *)

(*Plot average #triangles vs n*)

averages=Map[{#[[1]],Mean[#[[2,1;;-1]]]}&,gsraw];

averagesGrouped=GatherBy[averages,{#[[1,1]]&,#[[1,2]]&}];

(* averagesGroupes[[ n index, tau index , run index, {ntau, tri, ds} ]] *)

averagesGrouped=Transpose[averagesGrouped];

(* averagesGrouped[[ tau index, n index, run index ]] *)

nlabels=Map["n = "<>ToString[#]&,averagesGrouped[[1,All,1,1,1]]];

taulabels=Map["tau = "<>ToString[#]&,averagesGrouped[[All,1,1,1,2]]];

{{#[[1,1,1]],Mean[#[[All,2]]]},

}&,averagesGrouped,{2}];

ErrorListPlot[averagesErrorBars,Joined->True,PlotMarkers->Automatic,AxesLabel->{"n","\[LeftAngleBracket]triangles\[RightAngleBracket]"},PlotLegends->taulabels]

(* ::Subsection:: *)

(*Plot #triangles distribution for specific (n,tau)*)

histograms=Map[Histogram[#[[All,2]]]&,averagesGrouped,{2}];

nlabels=Map["n = "<>ToString[#]&,averagesGrouped[[1,All,1,1,1]]];

taulabels=Map["tau = "<>ToString[#]&,averagesGrouped[[All,1,1,1,2]]];

TableForm[histograms,TableHeadings->{taulabels,nlabels}]

#include "powerlaw.hpp"

    // Generate a random degree sequence

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

    float tauValues[] = {2.2f, 2.5f, 2.8f};

    for (int numVertices = 200; numVertices <= 2000; numVertices += 200) {

        for (float tau : tauValues) {

            DegreeSequence ds(numVertices);

            powerlaw_distribution degDist(tau, 2, numVertices);

            //std::poisson_distribution<> degDist(12);

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

            // the degree distribution

            for (int degreeSample = 0; degreeSample < 200; ++degreeSample) {

                // Generate a graph

                // might require multiple tries

                for (int i = 1; ; ++i) {

                    std::generate(ds.begin(), ds.end(),

                                  [&degDist, &rng] { return degDist(rng); });

                    if (g.createFromDegreeSequence(ds))

                        break;

                    // When 10 tries have not worked, output a warning

                    if (i % 10 == 0) {

                        std::cerr << "Warning: could not create graph from "

                                     "degree sequence. Trying again...\n";

                    }

                SwitchChain chain;

                if (!chain.initialize(g)) {

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

                    return 1;

                }

                std::cout << "Starting switch Markov chain with n = "

                          << numVertices << ", tau = " << tau << ". \t"

                          << std::flush;

                constexpr int mixingTime = 15000;

                constexpr int measureTime = 10000;

                constexpr int measureSkip =

                    100; // Take a sample every 100 steps

                constexpr int measurements =

                    (measureTime - 1) / measureSkip + 1;

                int movesDone = 0;

                int triangles[measurements];

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

                    if (chain.doMove())

                        ++movesDone;

                }

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

                    if (chain.doMove())

                        ++movesDone;

                    if (i % measureSkip == 0)

                        triangles[i / measureSkip] = chain.g.countTriangles();

                }

                std::cout << movesDone << '/' << mixingTime + measureTime

                          << " moves succeeded." << std::endl;

                if (outputComma)

                    outfile << ',';

                outputComma = true;

                std::sort(ds.begin(), ds.end());

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

                outfile << ',' << triangles << ',' << ds << '}' << std::flush;

            }