INCLUDES += -I.

all: switchchain switchchain_exponent switchchain_initialtris

switchchain:

switchchain_exponent:

switchchain_initialtris:

# target : dep1 dep2 dep3

# 	$@ = target

# 	$< = dep1

# 	$^ = dep1 dep2 dep3

        // Consider one of the three possible permutations

        // 1) e1.u - e1.v and e2.u - e2.v (original)

        // 2) e1.u - e2.u and e1.v - e2.v

        // 3) e1.u - e2.v and e1.v - e2.u

        bool switchType = permutationDistribution(mt);

        return g.exchangeEdges(e1index, e2index, switchType);

    }

    Graph g;

    std::mt19937 mt;

    std::uniform_int_distribution<> edgeDistribution;

    //std::uniform_int_distribution<> permutationDistribution;

    std::bernoulli_distribution permutationDistribution;

};

void getTriangleDegrees(const Graph& g) {

    std::vector<std::array<std::size_t,3>> trids;

    const auto& adj = g.getAdj();

    int triangles = 0;

    for (auto& v : adj) {

        for (unsigned int i = 0; i < v.size(); ++i) {

            for (unsigned int j = i + 1; j < v.size(); ++j) {

                if (g.hasEdge({v[i], v[j]})) {

                    ++triangles;

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

                    trids.push_back(ds);

                }

            }

        }

    }

    assert(triangles % 3 == 0);

    return;

}

//

// Assumes degree sequence does NOT contain any zeroes!

//

// method2 = true  -> take highest degree and finish its pairing completely

// method2 = false -> take new highest degree after every pairing

    // Similar to Havel-Hakimi but instead of pairing up with the highest ones

    // that remain, simply pair up with random ones

    unsigned int n = ds.size();

    // degree, vertex index

    std::vector<std::pair<unsigned int, unsigned int>> degrees(n);

    for (unsigned int i = 0; i < n; ++i) {

        degrees[i].first = ds[i];

        degrees[i].second = i;

    }

    std::vector<decltype(degrees.begin())> available;

    available.reserve(n);

    // Clear the graph

    g.reset(n);

    while (!degrees.empty()) {

        std::shuffle(degrees.begin(), degrees.end(), rng);

        // Get the highest degree:

        // If there are multiple highest ones, we pick a random one,

        // ensured by the shuffle.

        // The shuffle is needed anyway for the remaining part

        unsigned int dmax = 0;

            if (dmax == 1) {

                // Remove both

                // Erasure invalidates all iterators AFTER the erased one

                if (vIter > uIter) {

                    degrees.erase(vIter);

                    degrees.erase(uIter);

                } else {

                    degrees.erase(uIter);

                    degrees.erase(vIter);

                }

            } else {

                // Remove only v if it reaches zero

                uIter->first--;

                vIter->first--;

                if (vIter->first == 0)

                    degrees.erase(vIter);

            }

            //degrees.erase(std::remove_if(degrees.begin(), degrees.end(),

            //                             [](auto x) { return x.first == 0; }));

        }

    }

    return true;

}

    // Generate a random degree sequence

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

    // Goal:

    // Degrees follow a power-law distribution with some parameter tau

    // Expect:  #tri = const * n^{ something }

    // The goal is to find the 'something' by finding the number of triangles

    // for different values of n and tau

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

    Graph g;

    Graph g1;

    Graph g2;

    outfile << '{';

    bool outputComma = false;

        for (float tau : tauValues) {

            DegreeSequence ds(numVertices);

            powerlaw_distribution degDist(tau, 1, numVertices);

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

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

            // the degree distribution.

            // 500 samples seems to give reasonable results

                // Generate a graph

                // might require multiple tries

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

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

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

                    // First make the sum even

                    unsigned int sum = std::accumulate(ds.begin(), ds.end(), 0);

                    if (sum % 2) {

                        continue;

                        // Can we do this: ??

                        ds.back()++;

                    }

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

                    }

                }

                //

                // Test the GCM1 and GCM2 success rate

                //

                std::vector<int> greedyTriangles1;

                std::vector<int> greedyTriangles2;

                int successrate1 = 0;

                int successrate2 = 0;

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

                    Graph gtemp;

                    // Take new highest degree every time

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

                        ++successrate1;

                        greedyTriangles1.push_back(gtemp.countTriangles());

                        g1 = gtemp;

                    }

                    // Finish all pairings of highest degree first

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

                        ++successrate2;

                        greedyTriangles2.push_back(gtemp.countTriangles());

                        g2 = gtemp;

                    }

                }

                SwitchChain chain;

                if (!chain.initialize(g)) {

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

                    return 1;

                }

                std::cout << "Running n = " << numVertices << ", tau = " << tau

                          << ". \t" << std::flush;

                //int mixingTime = (32.0f - 26.0f*(tau - 2.0f)) * numVertices; //40000;

                //constexpr int measurements = 50;

                //constexpr int measureSkip =

                //    200; // Take a sample every ... steps

                int mixingTime = 0;

                constexpr int measurements = 50000;

                constexpr int measureSkip = 1;

                int triangles[measurements];

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

                }

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

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

                }

                if (outputComma)

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

                outputComma = true;

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

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

                outfile << ',' << triangles;

                outfile << ',' << ds;

                outfile << ',' << greedyTriangles1;

                outfile << ',' << greedyTriangles2;

                    movesDone = 0;

                    SwitchChain& c = chain1;

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

                        if (c.doMove())

                            ++movesDone;

                    }

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

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

                            if (c.doMove())

                                ++movesDone;

                        triangles[i] = c.g.countTriangles();

                    }

                    std::cout << movesDone << '/' << mixingTime + measurements * measureSkip

                        << " moves succeeded ("

                        << 100.0f * float(movesDone) /

                        float(mixingTime + measurements * measureSkip)

                        << "%).";

                    outfile << ',' << triangles;

                }

                    movesDone = 0;

                    SwitchChain& c = chain2;

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

                        if (c.doMove())

                            ++movesDone;

                    }

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

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

                            if (c.doMove())

                                ++movesDone;

                        triangles[i] = c.g.countTriangles();

                    }

                    std::cout << movesDone << '/' << mixingTime + measurements * measureSkip

                        << " moves succeeded ("

                        << 100.0f * float(movesDone) /

                        float(mixingTime + measurements * measureSkip)

                        << "%).";

                    outfile << ',' << triangles;

                }

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

                std::cout << std::endl;

            }

        }

    }

    outfile << '}';

    return 0;

}

                          << " random edges but they keep conflicting.\n";

            }

        } while (edgeConflicts(g.getEdge(e1index), g.getEdge(e2index)));

        // Consider one of the three possible permutations

        // 1) e1.u - e1.v and e2.u - e2.v (original)

        // 2) e1.u - e2.u and e1.v - e2.v

        // 3) e1.u - e2.v and e1.v - e2.u

        bool switchType = permutationDistribution(mt);

        return g.exchangeEdges(e1index, e2index, switchType);

    }

    Graph g;

    std::mt19937 mt;

    std::uniform_int_distribution<> edgeDistribution;

    //std::uniform_int_distribution<> permutationDistribution;

    std::bernoulli_distribution permutationDistribution;

};

//

// Assumes degree sequence does NOT contain any zeroes!

//

// method2 = true  -> take highest degree and finish its pairing completely

// method2 = false -> take new highest degree after every pairing

    // Similar to Havel-Hakimi but instead of pairing up with the highest ones

    // that remain, simply pair up with random ones

    unsigned int n = ds.size();

    // degree, vertex index

    std::vector<std::pair<unsigned int, unsigned int>> degrees(n);

    for (unsigned int i = 0; i < n; ++i) {

        degrees[i].first = ds[i];

        degrees[i].second = i;

    }

    std::vector<decltype(degrees.begin())> available;

    available.reserve(n);

    // Clear the graph

    g.reset(n);

    while (!degrees.empty()) {

        std::shuffle(degrees.begin(), degrees.end(), rng);

        // Get the highest degree:

        // If there are multiple highest ones, we pick a random one,

        // ensured by the shuffle.

        // The shuffle is needed anyway for the remaining part

        unsigned int dmax = 0;

(* ::Subsection:: *)

(*Observations on `Greedy Configuration Model'*)

(* ::Text:: *)

(*- GCM1: Greedy Configuration Model 1: take highest degree and do a single pairing, then take new highest degree*)

(*- GCM2: Greedy Configuration Model 2: take highest degree and completely do all its pairings (at random)*)

(*The difference does not matter if one node is by far the highest.*)

(**)

(*The success rates, conditioned on the degree sequence being graphical, is almost always higher using GCM2. For certain degree sequences the success rate of GCM2 can be 0.9 higher than that of GCM1. (i.e. amost always works vs almost always fails).*)

(*For tau > ~2.3 the success rate of GCM2 seems to be higher than 80% for most sequences.*)

(*For tau < ~2.3 the success rate of GCM2 can drop to less than 10% for some sequences but for many sequences it is still larger than 80%.*)

(**)

(*Success rate of GCM seems to be correlated with mixing time from Erdos-Gallai: higher success rate implies lower mixing time.  *)

(**)

(*Initial #triangles in both GCM1 and GCM2 is always below the average #triangles whereas Erdos-Gallai is usually many times higher than average.*)

(* ::Section:: *)

(*Data import*)

(* gsraw=SortBy[gsraw,{#[[1,1]]&,#[[1,2]]&}]; (* Sort by n and then by tau. The {} forces a *stable* sort because otherwise Mathematica sorts also on triangle count and other things. *) *)

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

(* Data format: *)

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

(* datatype index:

1: {n,tau}

2: #triangles time sequence

3: degree sequence

4: GCM1 starting triangle counts

5: GCM2 starting triangle counts

6: GCM1 time sequence

7: GCM2 time sequence

*)

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

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

(* ::Subsection:: *)

(*Merge data*)

newData=Import[NotebookDirectory[]<>"data/graphdata_3.m"];

getProperty[ds1_]:=Module[{ds,n=Length[ds1],tmp=ConstantArray[0,{Length[ds1],Length[ds1]}]},

ds=N[ds1/Sqrt[N[Total[ds1]]]]; (* scale degrees by 1/Sqrt[total] *)

(* The next table contains unneeded entries, but its faster to have a square table for the sum *)

tmp=Table[1.-Exp[-ds[[i]]ds[[j]]],{i,1,n},{j,1,n}];

Sum[tmp[[i,j]]*tmp[[j,k]]*tmp[[i,k]],{i,3,n},{j,2,i-1},{k,1,j-1}] (* somehow i>j>k is about 60x faster than doing i<j<k !!! *)

(* This sparser table is slower

tmp=Table[1.-Exp[-ds[[i]]ds[[j]]],{i,1,n-1},{j,i+1,n}];

(* tmp[[a,b]] is now with  ds[[a]]*ds[[a+b]] *)

Sum[tmp[[i,j-i]]*tmp[[j,k-j]]*tmp[[i,k-i]],{i,1,n-2},{j,i+1,n-1},{k,j+1,n}]

*)

];

(* gdata[[ tau index, n index, run index , {ntau, #tris, ds} ]] *)

avgAndProp=ParallelMap[{getProperty[#[[3]]],Mean[#[[2,1;;-1]]]}&,gdata[[2,2,1;;100]]];

Show[ListPlot[avgAndProp,AxesOrigin->{0,0},AxesLabel->{"degree-sequence-property","<#triangles>"},AspectRatio->1],Plot[x,{x,1,1000}]]

(* ::Subsection:: *)

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

measureSkip=1;

minCount=Min[Map[Min[#[[2]]]&,selectedData]];

maxCount=Max[Map[Max[#[[2]]]&,selectedData]];

maxTime=Max[Map[Length[#[[2]]]&,selectedData]];

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

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

(* ::Subsection:: *)

(*Compute 'mixing time'*)

(* Compute average of last part and check when the value drops below that for the first time *)

getMixingTime[values_]:=Module[{avg},

    (* average over the last 20 percent *)

    avg=Mean[values[[-Round[Length[values]/5];;-1]]];

    FirstPosition[values,_?(#<=avg&)][[1]]

]

(* gdata[[ tau index, n index, run index , {ntau, #tris, ds} ]] *)

measureSkip=1;

mixingTimes=Map[{#[[1,1]],(1/#[[1,1]])measureSkip * getMixingTime[#[[2]]]}&,gdata,{3}];

mixingTimesBars=Map[

    {{#[[1,1]],Mean[#[[All,2]]]},ErrorBar[StandardDeviation[#[[All,2]]](*/Sqrt[Length[#]]*)]}&

,mixingTimes,{2}];

ErrorListPlot[mixingTimesBars,Joined->True,PlotMarkers->Automatic,AxesLabel->{"n","~~mixing time divided by n"},PlotLegends->taulabels]

(* For n fixed, look at function of tau *)

measureSkip=1;

mixingTimes=Map[{#[[1,2]],(1/#[[1,1]])measureSkip * getMixingTime[#[[2]]]}&,gdata,{3}];

mixingTimesBars=Map[

    {{#[[1,1]],Mean[#[[All,2]]]},ErrorBar[StandardDeviation[#[[All,2]]]]}&

(* ::Package:: *)

Needs["ErrorBarPlots`"]

(* ::Section:: *)

(*Triangle exponent*)

(* ::Text:: *)

(*Expected number of triangles is the following powerlaw*)

(* ::DisplayFormula:: *)

(*#triangles = const\[Cross]n^(T(\[Tau]))      where    T(\[Tau])=3/2 (3-\[Tau])*)

(* When importing from exponent-only-data file *)

gsraw=SortBy[gsraw,#[[1,1]]&]; (* Sort by n *)

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

(* averagesGrouped[[ tau index, n index, run index , {ntau, avgtri} ]] *)

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

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

averagesErrorBars=Map[

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

ErrorBar[StandardDeviation[#[[All,2]]]]

}&,averagesGrouped,{2}];

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

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

(* ::Subsection:: *)

(*Fitting the log-log-plot*)

loglogdata=Log[averagesErrorBars[[All,All,1]]];

fits=Map[Fit[#,{1,logn},logn]&,loglogdata];

fitsExtra=Map[LinearModelFit[#,logn,logn]&,loglogdata];

fitsExtra[[1]]["ParameterConfidenceIntervalTable"]

fitsExtra[[1]]["BestFitParameters"]

fitsExtra[[1]]["ParameterErrors"]

fitsExtra[[1]]["ParameterConfidenceIntervals"]

(* ::Subsection:: *)

(*Plot of T(\[Tau])*)

tauValues=averagesGrouped[[All,1,1,1,2]];

exponents=Transpose[{tauValues,fits[[All,2,1]]}];

(* ::Subsection:: *)

(*T(\[Tau]) including error bars*)

tauValues=averagesGrouped[[All,1,1,1,2]];

exponentsErrorBars=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&,

Transpose[{tauValues,fitsExtra}]];