TARGETS += switchchain_ccm_cputime

template <class T1, class T2>

std::ostream &operator<<(std::ostream &s, const std::pair<T1, T2>& p) {

    s << '{' << p.first << ',' << p.second << '}';

    return s;

}

template< class RandomIt, class Compare >

void insertionSort( RandomIt first, RandomIt last, Compare comp ) {

    for (RandomIt next = first;;) {

        RandomIt a = next;

        next++;

        if (next == last)

            break;

        RandomIt b = next;

        auto newvalue = *next;

        //                next

        // 1 2 3 4  5  6   4

        //             a   b

        //          a  b

        //       a  b

        // a b

        while (b != first && comp(newvalue, *a)) { // if newvalue < *a

            *b = *a;

            a--;

            b--;

        }

        *b = newvalue;

    }

}

            insertionSort(degrees.begin(), degrees.end(),

#include "exports.hpp"

#include "graph.hpp"

#include "graph_ccm.hpp"

#include "graph_powerlaw.hpp"

#include "switchchain.hpp"

#include <algorithm>

#include <chrono>

#include <fstream>

#include <iostream>

#include <numeric>

#include <random>

#include <vector>

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

    // Simulation parameters

    const int numVerticesMin = 1000;

    const int numVerticesMax = 1000;

    const int numVerticesStep = 1000;

    //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 = 10;

    const int totalDegreeSamples = 1;

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

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

    };

    // Output file

    std::ofstream outfile;

    if (argc >= 2)

        outfile.open(argv[1]);

    else

        outfile.open("graphdata_ccm_cputime.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 << "degreeSamples: " << totalDegreeSamples << std::endl;

    outfile << "canonical ds" << std::endl;

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

    outfile << "measurements: full time evol\n";

    outfile << "data:\n";

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

    outfile << "2: edges\n";

    outfile << "3: HH timed triangle seq\n";

    outfile << "4: {ccm1 failed attempts, timed triangle seq}\n";

    outfile << "5: {ccm2 failed attempts, timed triangle seq}\n";

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

    // Mathematica does not accept normal scientific notation

    outfile << std::fixed;

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

    bool outputComma = false;

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

    Graph g;

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

         numVertices += numVerticesStep) {

        for (float tau : tauValues) {

            int mixingTime = getMixingTime(numVertices, tau);

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

                generateCanonicalPowerlawGraph(numVertices, tau, g, ds);

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

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

                SwitchChain chain;

                if (!chain.initialize(g, true)) {

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

                    return 1;

                }

                std::vector<std::pair<double, unsigned int>> triangleSeq(mixingTime);

                {

                    // record start time

                    auto start = std::chrono::high_resolution_clock::now();

                    // Incorporate Havel-Hakimi time

                    g.createFromDegreeSequence(ds);

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

                        auto now = std::chrono::high_resolution_clock::now();

                        std::chrono::duration<double> dt = now - start;

                        triangleSeq[i] =

                            std::make_pair(dt.count(), chain.g.getTrackedTriangles());

                        chain.doMove(true);

                    }

                }

                std::cout << " Finished timed HH time evol." << std::flush;

                if (outputComma)

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

                outputComma = true;

                outfile << '{';

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

                outfile << ',' << g.edgeCount();

                outfile << ',' << triangleSeq;

                for (int ccmType = 1; ccmType <= 2; ++ccmType) {

                    bool ccmMethod = (ccmType == 1 ? false : true);

                    // record start time

                    auto start = std::chrono::high_resolution_clock::now();

                    bool failed = true;

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

                        Graph gtemp;

                        if (constrainedConfigurationModel(ds, gtemp, rng,

                                                          ccmMethod)) {

                            chain.initialize(gtemp, true);

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

                                auto now =

                                    std::chrono::high_resolution_clock::now();

                                std::chrono::duration<double> dt = now - start;

                                triangleSeq[i] = std::make_pair(

                                    dt.count(), chain.g.getTrackedTriangles());

                                chain.doMove(true);

                            }

                            outfile << ',' << '{' << i << ',' << triangleSeq << '}';

                            failed = false;

                            break;

                        }

                    }

                    if (failed)

                        outfile << ",{1000,{}}";

                }

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

                std::cout << " Finished timed CCM time evols." << std::flush;

                std::cout << std::endl;

            }

        }

    }

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

    return 0;

}

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

(* ::Subsection:: *)

(*New code with CPU time*)

getCombinedTimeData[run_]:=Module[{maxTime,skipPts,hhData,ccm1Data,ccm2Data},

maxTime=Length[run[[3]]];

skipPts=Max[1,Round[maxTime/500]];

hhData=  {{0,run[[3,1,2]]}}~Join~run[[3,1;;-1;;skipPts]];

ccm1Data={{0,run[[4,2,1,2]]}}~Join~run[[4,2,1;;-1;;skipPts]];

ccm2Data={{0,run[[5,2,1,2]]}}~Join~run[[5,2,1;;-1;;skipPts]];

{Legended[hhData,"\[Tau] = "<>ToString[run[[1,2]]]],ccm1Data,ccm2Data}

]

dataSets=Map[getCombinedTimeData,gdata,{3}];

dataSetsFlattened=Flatten[dataSets,3];

colorList=Table[ColorData[97,"ColorList"][[1+Floor[i/3]]],{i,0,Length[dataSetsFlattened]-1}];

ListPlot[dataSetsFlattened[[1;;3]],Joined->True,PlotRange->All]

z2

plot1=ListPlot[dataSetsFlattened,Joined->True,PlotRange->{All,All},PlotStyle->colorList,ImageSize->300,PlotLabel->nlabels[[1]],Frame->True,FrameLabel->{"seconds","number of triangles"}]

Export[NotebookDirectory[]<>"plots/timeevol_ccm.pdf",plot1]

Export[NotebookDirectory[]<>"plots/timeevol_ccm_log.pdf",plot1log]