INCLUDES += -I.

all: switchchain switchchain_exponent switchchain_initialtris

switchchain:

switchchain_exponent:

switchchain_initialtris:

    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;

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

                }

                // std::cout << std::endl;

                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;

}

    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;

(*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

(* 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;

(*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}];

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}]];