diff --git a/cpp/Makefile b/cpp/Makefile index 3da471776c31b2c52446224ba5263e8379e04c3d..0a5064bf92997966049497679c05f956bc3607c8 100644 --- a/cpp/Makefile +++ b/cpp/Makefile @@ -12,6 +12,7 @@ CXXFLAGS += -Wno-int-in-bool-context TARGETS += switchchain TARGETS += switchchain_canonical_creationfreqs +TARGETS += switchchain_canonical_mixingtime TARGETS += switchchain_canonical_properties TARGETS += switchchain_ccm_constructionrate TARGETS += switchchain_ccm_cputime diff --git a/cpp/switchchain_canonical_mixingtime.cpp b/cpp/switchchain_canonical_mixingtime.cpp new file mode 100644 index 0000000000000000000000000000000000000000..cfd5fe08e37e6db0bba97084220234f33facc37b --- /dev/null +++ b/cpp/switchchain_canonical_mixingtime.cpp @@ -0,0 +1,144 @@ +#include "exports.hpp" +#include "graph.hpp" +#include "graph_powerlaw.hpp" +#include "graph_spectrum.hpp" +#include "switchchain.hpp" +#include +#include +#include +#include +#include +#include + +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 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>> samples; + for (int maxTime = numVertices; maxTime <= 20 * numVertices; + maxTime += numVertices) { + samples.push_back(make_pair(maxTime, std::vector())); + } + + 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 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; +} diff --git a/triangle_canonical_mixingtime.m b/triangle_canonical_mixingtime.m new file mode 100644 index 0000000000000000000000000000000000000000..cb8ff6ad2e8391563f88f11de4691def6141c1e9 --- /dev/null +++ b/triangle_canonical_mixingtime.m @@ -0,0 +1,57 @@ +(* ::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] + + +