diff --git a/cpp/Makefile b/cpp/Makefile index 452b979d4981bbf23fd0d533c5221afbc9680a93..bed50277fc98818587ad053622dab1ef4091c136 100644 --- a/cpp/Makefile +++ b/cpp/Makefile @@ -11,8 +11,10 @@ CXXFLAGS += -DNDEBUG CXXFLAGS += -Wno-int-in-bool-context TARGETS += switchchain +TARGETS += switchchain_canonical_properties TARGETS += switchchain_properties TARGETS += switchchain_exponent +TARGETS += switchchain_exponent_new TARGETS += switchchain_initialtris TARGETS += switchchain_mixingtime TARGETS += switchchain_spectrum diff --git a/cpp/graph_powerlaw.hpp b/cpp/graph_powerlaw.hpp index 9d6981eaef64d669f6d06bdf6f604698a532eaa7..bd6eaee6790916f5583bac67db0ce83135ba6d93 100644 --- a/cpp/graph_powerlaw.hpp +++ b/cpp/graph_powerlaw.hpp @@ -33,3 +33,22 @@ void generatePowerlawGraph(int n, float tau, Graph& g, DegreeSequence& ds, } } } + +void generateCanonicalPowerlawGraph(int n, float tau, Graph &g, + DegreeSequence &ds) { + ds.resize(n); + powerlaw_distribution degDist(tau, 1, n); + degDist.getFixedDistribution(n, ds.begin()); + + unsigned int sum = std::accumulate(ds.begin(), ds.end(), 0); + if (sum % 2) { + //TODO: Can we do this: ?? + ds.back()++; + } + + if (g.createFromDegreeSequence(ds)) + return; + g.reset(n); + std::cerr << "ERROR: Canonical ds is not graphical" << std::endl; +} + diff --git a/cpp/powerlaw.hpp b/cpp/powerlaw.hpp index 5dd0194ce349f77bf3f82f642f43d0e0bd2656eb..b325a2ab534a2dbeb29cde138d8d310b5e9a9817 100644 --- a/cpp/powerlaw.hpp +++ b/cpp/powerlaw.hpp @@ -1,10 +1,23 @@ #pragma once #include #include +#include #include // Discrete powerlaw distribution // obtained by rounding a continuous powerlaw distribution +// Let m = minimum and M = maximum and a = alpha +// Then for continuous version we have +// PDF: f(x) = C x^{-a} +// with C = (1-a)/(M^{1-a} - m^{1-a}) +// This gives CDF: F(x) = (x^{1-a} - m^{1-a}) / (M^{1-a} - m^{1-a}) +// with inversion F^{-1}(y) = ( y M^{1-a} + (1-y) m^{1-a} )^{1/(1-a)} +// i.e. linear interpolate between M^{1-a} < m^{1-a} +// For M = infty this becomes +// F'^{-1}(y) = m (1-y)^{1/(1-a)} +// = m (y')^{-1/(a-1)} +// where higher y' means lower x + class powerlaw_distribution { public: // xmin and xmax are inclusive @@ -29,6 +42,30 @@ class powerlaw_distribution { return x; } + // Generate non-random ``canonical'' powerlaw distribution + // Same as above operator but where the random number between + // 0 and 1 is replaced by stepping from 0 to 1 in fixed steps + template + void getFixedDistribution(int n, FwdIterator output) const { + // go in linear steps over interval + // [ M^{1-a} , m^{1-a} ] + double minVal = std::pow(double(xmax), 1.0f - alpha); + double maxVal = std::pow(double(xmin), 1.0f - alpha); + unsigned int x; + // Now it outputs in increasing order + for (int i = 0; i < n; ++i) { + double r1 = double(i) / double(n); + double r2 = r1 * minVal + (1.0 - r1) * maxVal; + x = std::round(std::pow(r2, _exponent)); + if (x > xmax || x < xmin) { + std::cerr << "ERROR: x not in [xmin,xmax] : " << x + << " not in [" << xmin << ',' << xmax + << "]; i = " << i << std::endl; + } + *output++ = x; + } + } + private: float alpha; double _exponent; diff --git a/cpp/switchchain_canonical_properties.cpp b/cpp/switchchain_canonical_properties.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d58097b4b869d31eb64dcaee078f62557297f1d6 --- /dev/null +++ b/cpp/switchchain_canonical_properties.cpp @@ -0,0 +1,183 @@ +#include "exports.hpp" +#include "graph.hpp" +#include "graph_powerlaw.hpp" +#include "graph_spectrum.hpp" +#include "switchchain.hpp" +#include +#include +#include +#include +#include +#include + +double getDSTN(const DegreeSequence& ds) { + std::vector> vals(ds.size()); + for (auto& v : vals) { + v.resize(ds.size(), 0); + } + + auto D = 0u; + for (auto d : ds) + D += d; + + double factor = 1.0 / double(D); + + for (auto i = 0u; i < ds.size(); ++i) { + for (auto j = i + 1; j < ds.size(); ++j) { + vals[i][j] = 1.0 - std::exp(-(ds[i] * ds[j] * factor)); + } + } + + double result = 0.0; + for (auto i = 0u; i < ds.size(); ++i) { + for (auto j = i + 1; j < ds.size(); ++j) { + for (auto k = j + 1; k < ds.size(); ++k) { + result += vals[i][j] * vals[j][k] * vals[i][k]; + } + } + } + return result; +} + +int main(int argc, char* argv[]) { + // Simulation parameters + const int numVerticesMin = 1000; + const int numVerticesMax = 10000; + const int numVerticesStep = 1000; + + float tauValues[] = {2.1f, 2.2f, 2.3f, 2.4f, 2.5f, 2.6f, 2.7f, 2.8f, 2.9f}; + + //const int totalDegreeSamples = 5000; + + auto getMixingTime = [](int n, float tau) { + return int(50.0f * (50.0f - 30.0f * (tau - 2.0f)) * n); + }; + constexpr int measurements = 10; + constexpr int measureSkip = 1000; // Take a sample every ... steps + + // Output file + std::ofstream outfile; + if (argc >= 2) + outfile.open(argv[1]); + else + outfile.open("graphdata_canonical_properties.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 << "mixingTime: 50 * (50 - 30 (tau - 2)) n\n"; + outfile << "data:\n"; + outfile << "1: {n,tau}\n"; + outfile << "2: avgTriangles\n"; + outfile << "3: edges\n"; + outfile << "4: dstn\n"; + outfile << "5: { HH A, HH L, average A, average L } where for each there is (average of) {lambda1 , lambda1 - lambda2, lambda1/lambda2}\n"; + outfile << "6: switching successrate after mixing\n"; + outfile << "7: initial HH triangles\n"; + outfile << "*)" << std::endl; + + // Mathematica does not accept normal scientific notation + outfile << std::fixed; + outfile << '{'; + bool outputComma = false; + + Graph g; + for (int numVertices = numVerticesMin; numVertices <= numVerticesMax; + numVertices += numVerticesStep) { + for (float tau : tauValues) { + DegreeSequence ds; + generateCanonicalPowerlawGraph(numVertices, tau, g, ds); + + SwitchChain chain; + if (!chain.initialize(g)) { + std::cerr << "Could not initialize Markov chain.\n"; + return 1; + } + + std::cout << "Running (n,tau) = (" << numVertices << ',' << tau + << "). " << std::flush; + + // Mix + int mixingTime = getMixingTime(numVertices, tau); + for (int i = 0; i < mixingTime; ++i) { + chain.doMove(); + } + + std::cout << "Mixing done. " << std::flush; + + std::array HHAspectrum; + std::array HHLspectrum; + std::array avgAspectrum; + std::array avgLspectrum; + + auto getSpectralValues = + [](const std::vector &s) -> std::array { + auto l1 = s[s.size() - 1]; + auto l2 = s[s.size() - 2]; + return {l1, l1 - l2, l1 / l2}; + }; + + GraphSpectrum gs_start(g); + GraphSpectrum gs(chain.g); + + HHAspectrum = + getSpectralValues(gs_start.computeAdjacencySpectrum()); + HHLspectrum = + getSpectralValues(gs_start.computeLaplacianSpectrum()); + + long long trianglesTotal = 0; + int movesDone = 0; + avgAspectrum.fill(0); + avgLspectrum.fill(0); + for (int i = 0; i < measurements; ++i) { + for (int j = 0; j < measureSkip; ++j) + if (chain.doMove()) + ++movesDone; + trianglesTotal += chain.g.countTriangles(); + auto sA = getSpectralValues(gs.computeAdjacencySpectrum()); + auto sL = getSpectralValues(gs.computeLaplacianSpectrum()); + for (auto i = 0u; i < 3; ++i) { + avgAspectrum[i] += sA[i]; + avgLspectrum[i] += sL[i]; + } + } + float avgTriangles = float(trianglesTotal) / float(measurements); + float successrate = + float(movesDone) / float(measurements * measureSkip); + for (auto &f : avgAspectrum) + f /= float(measurements); + for (auto &f : avgLspectrum) + f /= float(measurements); + + std::cout << "Measuring done." << std::flush; + + if (outputComma) + outfile << ',' << '\n'; + outputComma = true; + + outfile << '{' << '{' << numVertices << ',' << tau << '}'; + outfile << ',' << avgTriangles; + outfile << ',' << g.edgeCount(); + outfile << ',' << getDSTN(ds); + outfile << ',' << '{' << HHAspectrum; + outfile << ',' << HHLspectrum; + outfile << ',' << avgAspectrum; + outfile << ',' << avgLspectrum; + outfile << '}'; + outfile << ',' << successrate; + outfile << ',' << g.countTriangles(); + outfile << '}' << std::flush; + + std::cout << "Output done." << std::endl; + } + } + outfile << '}'; + return 0; +} diff --git a/cpp/switchchain_exponent_new.cpp b/cpp/switchchain_exponent_new.cpp new file mode 100644 index 0000000000000000000000000000000000000000..3366283f49ba09c0512e994f03e0422951159d27 --- /dev/null +++ b/cpp/switchchain_exponent_new.cpp @@ -0,0 +1,112 @@ +#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 = 20000; + const int numVerticesMax = 50000; + const int numVerticesStep = 10000; + + //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 = 1000; + + auto getMixingTime = [](int n, float tau) { + return int(50.0f * (50.0f - 30.0f * (tau - 2.0f)) * n); + }; + constexpr int measurements = 10; + constexpr int measureSkip = 1000; // Take a sample every ... steps + + // Output file + std::ofstream outfile; + if (argc >= 2) + outfile.open(argv[1]); + else + outfile.open("graphdata_exponent_new.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 << "mixingTime: 50 * (50 - 30 (tau - 2)) n\n"; + outfile << "data:\n"; + outfile << "1: {n,tau}\n"; + outfile << "2: avgTriangles\n"; + outfile << "*)" << std::endl; + + // Mathematica does not accept normal scientific notation + outfile << std::fixed; + outfile << '{'; + bool outputComma = false; + + std::mt19937 rng(std::random_device{}()); + Graph g; + for (int numVertices = numVerticesMin; numVertices <= numVerticesMax; + numVertices += numVerticesStep) { + for (float tau : tauValues) { + // 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); + + SwitchChain chain; + if (!chain.initialize(g)) { + std::cerr << "Could not initialize Markov chain.\n"; + return 1; + } + + std::cout << "Running (n,tau) = (" << numVertices << ',' + << tau << "). " << std::flush; + + // Mix + int mixingTime = getMixingTime(numVertices, tau); + for (int i = 0; i < mixingTime; ++i) { + chain.doMove(); + } + + std::cout << "Mixing done. " << std::flush; + + long long trianglesTotal = 0; + for (int i = 0; i < measurements; ++i) { + for (int j = 0; j < measureSkip; ++j) + chain.doMove(); + trianglesTotal += chain.g.countTriangles(); + } + float avgTriangles = + float(trianglesTotal) / float(measurements); + + std::cout << "Measuring done." << std::flush; + + if (outputComma) + outfile << ',' << '\n'; + outputComma = true; + + outfile << '{' << '{' << numVertices << ',' << tau << '}'; + outfile << ',' << avgTriangles; + outfile << '}' << std::flush; + + std::cout << "Output done." << std::endl; + } + } + } + outfile << '}'; + return 0; +} diff --git a/triangle_exponent_plots.m b/triangle_exponent_plots.m index b88398c79d7059660338a64c3be02f5d1493305d..230af3640f86d2ac3426dabf6366eaff7a62f246 100644 --- a/triangle_exponent_plots.m +++ b/triangle_exponent_plots.m @@ -15,7 +15,10 @@ Needs["ErrorBarPlots`"] (*#triangles = const\[Cross]n^(T(\[Tau])) where T(\[Tau])=3/2 (3-\[Tau])*) -(* When importing from exponent-only-data file *) +(* When importing from exponent-only-data file or property-data file *) +(* graphdata_exponent_mix32.m *) +(* graphdata_exponent_highN.m *) +(* graphdata_properties2.m *) gsraw=Import[NotebookDirectory[]<>"data/graphdata_exponent_mix32.m"]; gsraw=SortBy[gsraw,#[[1,1]]&]; (* Sort by n *) averagesGrouped=GatherBy[gsraw,{#[[1,2]]&,#[[1,1]]&}]; @@ -40,7 +43,9 @@ ListLogLogPlot[averagesErrorBars[[All,All,1]],Joined->True,PlotMarkers->Automati (*Fitting the log-log-plot*) -loglogdata=Log[averagesErrorBars[[All,All,1]]]; +nRange=5;;-1; +nRange=All; +loglogdata=Log[averagesErrorBars[[All,nRange,1]]]; fits=Map[Fit[#,{1,logn},logn]&,loglogdata]; fitsExtra=Map[LinearModelFit[#,logn,logn]&,loglogdata]; diff --git a/triangle_spectrum_plots.m b/triangle_spectrum_plots.m index 2fbbd2afd660a1d7fb4b82a8b5c7240057f93cf4..e525f1aa7a49ff7c4cda6e15d7c26d53a00a96b5 100644 --- a/triangle_spectrum_plots.m +++ b/triangle_spectrum_plots.m @@ -1,6 +1,6 @@ (* ::Package:: *) -gsraw=Import[NotebookDirectory[]<>"cpp/graphdata_spectrum.m"]; +gsraw=Import[NotebookDirectory[]<>"data/graphdata_properties.m"]; (* 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. *) *) @@ -9,19 +9,20 @@ gdata=GatherBy[gsraw,{#[[1,2]]&,#[[1,1]]&}]; (* gdata[[ tau index, n index, run index , datatype index ]] *) (* datatype index: 1: {n,tau} -2: etmt +2: avgTriangles +3: edges +4: dstn +5: { HH A, HH L, average A, average L } where for each there is (average of) {lambda1 , lambda1 - lambda2, lambda1/lambda2} +6: switching successrate after mixing +7: initial HH triangles *) tauvalues=gdata[[All,1,1,1,2]]; nlabels=Map["n = "<>ToString[#]&,gdata[[1,All,1,1,1]]]; taulabels=Map["\[Tau] = "<>ToString[#]&,gdata[[All,1,1,1,2]]]; -ListPlot[gdata[[1,1,1,{4,6}]]] -ListPlot[gdata[[1,1,1,{5,7}]]] - - -Histogram[gdata[[1,1,1,{4,6}]],50] -Histogram[gdata[[1,1,1,{5,7}]],50] +Histogram[gdata[[1,1,All,5,3,1]]] +Histogram[gdata[[1,1,All,5,3,2]]]