From 2c86fdd4152ed34533baf517172afa13e9b5cd5c 2017-07-10 18:38:29 From: Tom Bannink Date: 2017-07-10 18:38:29 Subject: [PATCH] Add start of ccm time evol --- diff --git a/cpp/Makefile b/cpp/Makefile index bed50277fc98818587ad053622dab1ef4091c136..e4cb5af5c37838f2f7d8882049d02064b6334bbb 100644 --- a/cpp/Makefile +++ b/cpp/Makefile @@ -12,10 +12,11 @@ CXXFLAGS += -Wno-int-in-bool-context TARGETS += switchchain TARGETS += switchchain_canonical_properties +TARGETS += switchchain_ccm_initialtris +TARGETS += switchchain_ccm_timeevol TARGETS += switchchain_properties TARGETS += switchchain_exponent TARGETS += switchchain_exponent_new -TARGETS += switchchain_initialtris TARGETS += switchchain_mixingtime TARGETS += switchchain_spectrum TARGETS += switchchain_successrates diff --git a/cpp/graph_gcm.hpp b/cpp/graph_ccm.hpp similarity index 91% rename from cpp/graph_gcm.hpp rename to cpp/graph_ccm.hpp index b53a8acb6625f4f23ac6b1c5fc4ad33ee8c464a3..c3033c2136b8e5af284533d20fa6298a5aa42848 100644 --- a/cpp/graph_gcm.hpp +++ b/cpp/graph_ccm.hpp @@ -10,7 +10,7 @@ // method2 = true -> take highest degree and finish its pairing completely // method2 = false -> take new highest degree after every pairing template -bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, RNG& rng, bool method2) { +bool constrainedConfigurationModel(DegreeSequence& ds, Graph& g, RNG& rng, bool method2) { // 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(); @@ -47,7 +47,7 @@ bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, RNG& rng, bool metho return false; if (dmax == 0) { - std::cerr << "ERROR 1 in GCM.\n"; + std::cerr << "ERROR 1 in CCM.\n"; } unsigned int u = uIter->second; @@ -61,9 +61,9 @@ bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, RNG& rng, bool metho auto vIter = degrees.begin(); while (dmax--) { if (vIter->first == 0) - std::cerr << "ERROR in GCM2.\n"; + std::cerr << "ERROR in CCM2.\n"; if (!g.addEdge({u, vIter->second})) - std::cerr << "ERROR. Could not add edge in GCM2.\n"; + std::cerr << "ERROR. Could not add edge in CCM2.\n"; vIter->first--; if (vIter->first == 0) vIter = degrees.erase(vIter); @@ -84,9 +84,9 @@ bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, RNG& rng, bool metho auto vIter = available[distr(rng)]; // pair u to v if (vIter->first == 0) - std::cerr << "ERROR 2 in GCM1.\n"; + std::cerr << "ERROR 2 in CCM1.\n"; if (!g.addEdge({u, vIter->second})) - std::cerr << "ERROR. Could not add edge in GCM1.\n"; + std::cerr << "ERROR. Could not add edge in CCM1.\n"; // Purge anything with degree zero // Be careful with invalidating the other iterator! // Degree of u is always greater or equal to the degree of v diff --git a/cpp/switchchain.cpp b/cpp/switchchain.cpp index 930fedebbcfacbe96e1f3c7406b4f1238297b37e..9305d4286be94616747a376926624a08f049c430 100644 --- a/cpp/switchchain.cpp +++ b/cpp/switchchain.cpp @@ -1,7 +1,7 @@ #include "switchchain.hpp" #include "exports.hpp" #include "graph.hpp" -#include "graph_gcm.hpp" +#include "graph_ccm.hpp" #include "graph_powerlaw.hpp" #include "graph_spectrum.hpp" #include diff --git a/cpp/switchchain_ccm_initialtris.cpp b/cpp/switchchain_ccm_initialtris.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d62f9258837e25112454508dab2ba7b357731599 --- /dev/null +++ b/cpp/switchchain_ccm_initialtris.cpp @@ -0,0 +1,150 @@ +#include "exports.hpp" +#include "graph.hpp" +#include "graph_ccm.hpp" +#include "graph_powerlaw.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 = 2000; + const int numVerticesStep = 400; + + 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 = 200; + + auto getMixingTime = [](int n, float tau) { + return int(50.0f * (50.0f - 30.0f * (tau - 2.0f)) * n); + }; + auto getMeasurements = [](int n, float tau) { + (void)n; + (void)tau; + return 100; + }; + auto getMeasureSkip = [](int n, float tau) { + (void)tau; + return 10 * n; // Take a sample every ... steps + }; + + // Output file + std::ofstream outfile; + if (argc >= 2) + outfile.open(argv[1]); + else + outfile.open("graphdata_ccm_initialtris.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 << "measurements: 100\n"; + outfile << "measureSkip: 10 n\n"; + outfile << "data:\n"; + outfile << "1: {n,tau}\n"; + outfile << "2: avgTriangles\n"; + outfile << "3: {ccmTris1, ccmsrate1} \n"; + outfile << "4: {ccmTris2, ccmsrate2} \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); + + std::cout << "Running (n,tau) = (" << numVertices << ',' << tau + << "). " << std::flush; + + // + // Test the GCM1 and GCM2 success rate + // + long long gcmTris1 = 0; + long long gcmTris2 = 0; + int successrate1 = 0; + int successrate2 = 0; + for (int i = 0; i < 100; ++i) { + Graph gtemp; + // Take new highest degree every time + if (constrainedConfigurationModel(ds, gtemp, rng, false)) { + ++successrate1; + gcmTris1 += gtemp.countTriangles(); + } + // Finish all pairings of highest degree first + if (constrainedConfigurationModel(ds, gtemp, rng, true)) { + ++successrate2; + gcmTris2 += gtemp.countTriangles(); + } + } + + SwitchChain chain; + if (!chain.initialize(g)) { + std::cerr << "Could not initialize Markov chain.\n"; + return 1; + } + + + long long trianglesTotal = 0; + + std::cout << " Finished CCM generation." << std::flush; + + 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); + for (int i = 0; i < measurements; ++i) { + for (int j = 0; j < measureSkip; ++j) + chain.doMove(true); + trianglesTotal += chain.g.getTrackedTriangles(); + } + + std::cout << " Finished mixing and measurements." << std::flush; + + if (outputComma) + outfile << ',' << '\n'; + outputComma = true; + + float avgTriangles = + float(trianglesTotal) / float(measurements); + outfile << '{'; + outfile << '{' << numVertices << ',' << tau << '}'; + outfile << ',' << avgTriangles; + outfile << ',' << '{' << gcmTris1 << ',' << successrate1 << '}'; + outfile << ',' << '{' << gcmTris2 << ',' << successrate2 << '}'; + outfile << '}' << std::flush; + + std::cout << std::endl; + } + } + } + outfile << '}'; + return 0; +} diff --git a/cpp/switchchain_ccm_timeevol.cpp b/cpp/switchchain_ccm_timeevol.cpp new file mode 100644 index 0000000000000000000000000000000000000000..22db8bb0eea8b3bab2703031c368a214fd8109c3 --- /dev/null +++ b/cpp/switchchain_ccm_timeevol.cpp @@ -0,0 +1,154 @@ +#include "exports.hpp" +#include "graph.hpp" +#include "graph_ccm.hpp" +#include "graph_powerlaw.hpp" +#include "switchchain.hpp" +#include +#include +#include +#include +#include +#include + +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; + + 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_timeevol.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 triangle seq\n"; + outfile << "4: {ccm1 failed attempts, triangle seq}\n"; + outfile << "5: {ccm2 failed attempts, 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 triangleSeq(mixingTime); + for (int i = 0; i < mixingTime; ++i) { + chain.doMove(true); + triangleSeq[i] = chain.g.getTrackedTriangles(); + } + + std::cout << " Finished 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); +#if 0 + outfile << ',' << '{'; + for (int i = 0; i < 10; ++i) { + if (i != 0) + outfile << ','; + Graph gtemp; + if (constrainedConfigurationModel(ds, gtemp, rng, + ccmMethod)) { + chain.initialize(gtemp, true); + for (int i = 0; i < mixingTime; ++i) { + chain.doMove(true); + triangleSeq[i] = chain.g.getTrackedTriangles(); + } + outfile << triangleSeq; + } else { + outfile << '{' << '}'; + } + } + outfile << '}'; +#endif + bool failed = true; + for (int i = 0; i < 50; ++i) { + Graph gtemp; + if (constrainedConfigurationModel(ds, gtemp, rng, + ccmMethod)) { + chain.initialize(gtemp, true); + for (int i = 0; i < mixingTime; ++i) { + chain.doMove(true); + triangleSeq[i] = chain.g.getTrackedTriangles(); + } + outfile << ',' << '{' << i << ',' << triangleSeq << '}'; + failed = false; + break; + } + } + if (failed) + outfile << ",{50,{}}"; + } + + outfile << '}' << std::flush; + + std::cout << " Finished CCM time evols." << std::flush; + + std::cout << std::endl; + } + } + } + outfile << '\n' << '}'; + return 0; +} diff --git a/cpp/switchchain_initialtris.cpp b/cpp/switchchain_initialtris.cpp deleted file mode 100644 index 158e26550a13b9db01e08a9cb6d917f7be2afed4..0000000000000000000000000000000000000000 --- a/cpp/switchchain_initialtris.cpp +++ /dev/null @@ -1,115 +0,0 @@ -#include "exports.hpp" -#include "graph.hpp" -#include "graph_gcm.hpp" -#include "graph_powerlaw.hpp" -#include "switchchain.hpp" -#include -#include -#include -#include -#include -#include - -int main() { - // 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; - - std::ofstream outfile("graphdata_initialtris.m"); - outfile << '{'; - bool outputComma = false; - - for (int numVertices = 200; numVertices <= 2000; numVertices += 400) { - for (float tau : tauValues) { - // For a single n,tau take samples over several instances of - // the degree distribution. - for (int degreeSample = 0; degreeSample < 200; ++degreeSample) { - DegreeSequence ds; - generatePowerlawGraph(numVertices, tau, g, ds, rng); - - std::cout << "Running n = " << numVertices << ", tau = " << tau - << "." << std::flush; - - // - // Test the GCM1 and GCM2 success rate - // - long long gcmTris1 = 0; - long long gcmTris2 = 0; - 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; - gcmTris1 += gtemp.countTriangles(); - } - // Finish all pairings of highest degree first - if (greedyConfigurationModel(ds, gtemp, rng, true)) { - ++successrate2; - gcmTris2 += gtemp.countTriangles(); - } - } - - SwitchChain chain; - if (!chain.initialize(g)) { - std::cerr << "Could not initialize Markov chain.\n"; - return 1; - } - - int mixingTime = (32.0f - 20.0f * (tau - 2.0f)) * numVertices; - constexpr int measurements = 20; - constexpr int measureSkip = 200; - - int movesDone = 0; - - long long trianglesTotal = 0; - - std::cout << " .. \t" << std::flush; - - for (int i = 0; i < mixingTime; ++i) { - if (chain.doMove()) - ++movesDone; - } - for (int i = 0; i < measurements; ++i) { - for (int j = 0; j < measureSkip; ++j) - if (chain.doMove()) - ++movesDone; - trianglesTotal += chain.g.countTriangles(); - } - - std::cout << movesDone << '/' << mixingTime + measurements * measureSkip - << " moves succeeded (" - << 100.0f * float(movesDone) / - float(mixingTime + measurements * measureSkip) - << "%)."; - //std::cout << std::endl; - - if (outputComma) - outfile << ',' << '\n'; - outputComma = true; - - float avgTriangles = - float(trianglesTotal) / float(measurements); - outfile << '{'; - outfile << '{' << numVertices << ',' << tau << '}'; - outfile << ',' << avgTriangles; - outfile << ',' << '{' << gcmTris1 << ',' << successrate1 << '}'; - outfile << ',' << '{' << gcmTris2 << ',' << successrate2 << '}'; - outfile << '}' << std::flush; - - std::cout << std::endl; - } - } - } - outfile << '}'; - return 0; -} diff --git a/cpp/switchchain_spectrum.cpp b/cpp/switchchain_spectrum.cpp index 65dd3fb144ba2f158149bec11fc4388542bc36ec..de4ab7ca70bdd9ac54ded6577669578240b03de8 100644 --- a/cpp/switchchain_spectrum.cpp +++ b/cpp/switchchain_spectrum.cpp @@ -1,7 +1,7 @@ #include "switchchain.hpp" #include "exports.hpp" #include "graph.hpp" -#include "graph_gcm.hpp" +#include "graph_ccm.hpp" #include "graph_spectrum.hpp" #include "graph_powerlaw.hpp" #include diff --git a/plots/triangle_exponent.pdf b/plots/triangle_exponent.pdf index 539fb6a60ec4679e91bd31ba5829b3283faa1379..ebdb0f4fbb796b478a361cdf2ef07cb0f725df4a 100644 GIT binary patch delta 6075 zcmZWrby(Ehwx%5hl$OQ;r4bp18HRL_E`4k^BjQGzbF-NDZAz3(`F_l1fV> zAl>2hJLjI~-1~j&pS9MzSM7K2=h?d@7{4zX|D75zV5nNmT1PdP@%*c*&kLVJe-w^` zl7P@a%^YXG+a|<@X`Spr-s%||(#wt05kFDBGeT~1TJUQUouG3xo_)00xh0)o+G~Am zt2IdVu%fW;eER_~)pTK6qw)L8Hz~*8b+y|U5#RmmE>GDQDzDu4>LaEzJ|`U;&RT2F zHu6scSGv15>_wTk)3v9hohQbmP6y5_E;N-EmeP0ozLa${D93G?#5Frie3|u8%+h+P zEmXc&A0-k~X?D`(G``a{EGcOE-HG4&0p}<(SL&i-Maj7dwLhu)+yvn>z}UyUIpD*0 zQ|q|w(IUri((zbOzrounrHKQn`&WUnkL)S07%%*B8Du0B75n|UU_1?VT>c)iPiD*U@r1t*O zRU0FSitXl77?AztSkbeP7PNHGp^HG$5X90-!rHRJfbjeS_8WZ557=CUZ>Y}EY&oS_ z?f2~Vf%^zqp^op$iJpi0x@3i&f!2}%l(>7BLtDc^rebRF%Ke4=MFDfSOD*`Ii60 zDpr8UOA9+@+#bjuc-O~}eqmR(+j5R3nAGrJ=b62_F^K=%MN)$vSsT+&hzaHIvTY5t zi{Xt-QM|mBvpAJff46NA4N$4A96U+NA!fQ8qaK-ZB=FvxP#h%t%EZ=YeT*+e_NY$FlSJg7> zI34jst^1X+#BkZyXK6NwN5d@-cx)@mo;)|1a_lo(^)4NL5k&tnsR{tyR`XFaTa{Ef zNYuKCPy1;=37o4*)^$ISc$7^$No;6FM;MhC&o1L0pG~OVx=kO?`yW!^w+8wU>S3KV zL#g>D*ftJ>C)>(VT5j;{*Xpl>jl@TZ%arF0$i>FN6^3QzQ*KN-)=VTeD-Lgo9G(RV zuuiOg()5ZI0RBy)`S1io*y3zkpNXP(P%1k3)qPC{FEr{Zi1Jyi8am+_X z)T&G>k|zacM<@OYK-v#&q=JN9d1*-2lISp1D2wItuqOwS+FgEERaZY|iO&^JGiJ%J zqwrICx?kT}o*NuX89n}Wy+1qTxq7Q@dv@skJiutNwKIR^dbU(gAKp0%E`8H$UXNDj zibAr0eQ$VKN`Bl<9-os5=Jk|EQ@(>}4K5k;zS=O0(kK>-DT)9wson2L4NYxbbH{w) z)@@RRzx{Mik0mV-dnmSQL@V8RDP4r{%?papvid37L_+U0e9^(xA`>MeP;inJ#`0XH z5Gb5MjpciAM59$Z*=CSdXymd(U~@phkX;DNQ#yZYs<_vR39*j`k4Y%5qYf?82ap_4 z71Giy975abD74Lq=Q3qti18}n4LI8bKcSXEEGvZRrtoYV8bf#@aH zpi0#C5^#dEJ{{4X{FJhV9IGcse-z~byu8~t3+8s?w_sKecuE$!e(i<-NI+-Ic}Nr? zs&mW*a$OW}dVKDSd7Gp6)6aWnCv#3mYMO;*|HIFd@bmfkgg?>Y39HXkq{}nmw57X_l@i0$v9tVE5HDAMeyY?{O*;f0qJrO2bw`>BI{S72KuqCQ zNbe`~S^BWHREe;MbY#A>S^OJ0dT``wS+7maON!$DMvApads12Di512B7armn^XAeY z#Q-ItF&Y=#TR#X){Nw3H+OMT>?fn>r0rrPGoeBG|j2047GNxZV{D?@v6RUq$``$VV zzDKen-FyJ(yD2tyFlmz;sVjv7!n4fcC8u_GAkmwBa2cT=>Ya z%f5T5a^N7iPs-H2RThi5i)L;zq<1BoRO*!k-b61bjtUa&ci-GCzq#)i*5`-@8+FU4 zThu6e>pQlznbVEnLNSNIBpVu=vZi`>Jt8=)@|VqXpSq>?=!|42y=9@_4D0;b%lC7} z5V_7Z@s=h29O+KV<1z%K?GrkkpL4(osq`M@ENXCSUhOBhgY>l0rvP7VP*wCD9PVX% z58p_5V-*)Fwm2&NiJ#8;E$Uf;q6%X}x&`0Kg7Gaxi(bJrQ4wgyF-s&26dXrxG<~HM z?;F6y5$!xlR7Ry=2bYmP8`Z0d=OzC(yKC~*j&6C?^v)>6k_JFxdL6ULRH;KuYu!~# z9X(2nct=)w`eswp^ z{*H|B(axDTUCmhg@Q-+XiPHX79@%5&w*0W|t^~8jGP)dE(lX+Lo>;dCJcU^MxiUX> z#&;yuH4luc{PuvtWE9VUQW!DLQC7kOWI9x%D|}jW6pg!q(VFBb^AoduTt^I4773uN zwIRF6d+9;wU0U2%J|lWkU(mpTWoi#MyDj2pOOpHfa!KNy0JG6a`r1kxZV%+$-c5+h z_G6~KG7J2MU5+??A<`|Wy++DZuA>6>P)DnK1p=T;FCdmEV8SWKStBrhA>mz>5r4Gc z%!p6q7{@&b_u*Iw4+manS?05c0WCCFt>un&p+N#IPhVJ4Yf7(pTWJKcq@c3Ms5amh z96SE2i2d<B-h7XspeqBD!T(TGTK9*H!t!QKNk=e@p>o={D z`IVf__yElZlk6#pQI2`$+6^O8U%!yhU6k39mGH5YPNW2ikli8}w=H~x!s)Z;%rS!) zn=`5XdP(12u81YfRI9koyEdqGp~V`hncoQxGn`L=W6~J@JkPWUn!ADYw3W-Z#puv` z53>zwCn9Y#obdZk!cd=-(a~X40bRR}o>5c*K$p_3hft1J%A!K7i}cT)`aJ&ZUiv{6 zRMbmpRibf9k%!B$Jtwqba~%M4GGwn3S(@ANX+$+s8Gy zn#*NEm8wFIxfO_*+y*S4?Z!EQ6tI@&SNR?X6g^A$+EiX8?#H145Tdto>c)y58mKnl zVBy(nkLQz86mM02M|Y1Vl@LpotA-Jtur5d?qA2B27dUo03KxL8db^A=-IUVHSGQnn zDuJ#L<7^Z$GWIMP!wgUfH?@Fumafa zhx~Kh59|jGocXS?&%SdHr94&tjQ_$+&_Un;UM%zI)wX>OlgCK8ezJfOFudYZsR%oDnNP+EF=v`y?$Z&dA7=c1} zLFgaX@{q1L^D7j1^YcP7e6!n$ykZA3a7;M9;s2Uzo}wReiLVv;Q1i|HMH_sBVl49~ z^!?(`047T*&G6D2=;0$2ZxdUA%s0YZ2A`?N`CH}Q@>4NMscP`T<0t0%Y%)=|R=(V| zcpxaXaO8aFTD{@%Zhh~(Mj=H+Acl#BB()>v@^p22vKcsK-`_as{_H9EUi8lIiy^@G z>~v*%`edgmKWVpp1^5bEvBKnoSczc5P@zg3C>0+-2qVGL(qK1tS1U6|@T=5Py;)}# zJ-TL?M5PHJv8QOLq5R9bforz0-lNG_)Pu#d%U%Nhu6ijbnPJ`GWh0aw>%4)(_!}SQEUq-qTRbY`!aZ5bQ?Zr0r#O z)L$Q_i03u$7p>JHvJ}M4F`0`Nv|P>a+RLYrM28593$mI<%7= zUhAsE9|Ui|zYQ;4xxfAi^*-YLD*eOcq)S%d@Fep$YCTl7BXNmNZW2Ph(|uEY{d<2a zYY%JTa`dvR+RU?w&XiES>$8ap1F95c<2C<5!(Qw+&c$MhB|}8EM>XT|T&n-+aOjKO z7N=Rsj1P^MT}O4rK%VgM=H#l$71`6ny}C_jUmt|%?@EnN>2=rVR*-XAy9nl~obqm< zHtcMv_=_o!H2Ljd=T4`|2$w}fpcrfu5SfR6R7g!Cd2G?|d-Bn9aug`qY4%ZnnC-Lp zBjfqfiTYlimn!HWI+ohp5qr@4+-4wVatHPBi@%Lyv(pJVSEn zD|^K8nT))yhF4C0c(YPy5KqEZrLHib+{XQ5%7Q<1;AX0TUAU!YPuxdKgOP0l#w#1$ zV>@-HrudBHgSeL?Mc#2>;s?svO`LY2U*GNYl_~Hg+TBs}P#?jIp%6NF*@^*W8OD$4 z>=NkTi-AL#EJYaf19mQM+cD06R;F$6njg|qQ$B@i$rD+AA&ya-#@22_d^4&5i0a!4 zs`mmRs|JeuN{+5*S7>U!zBp+ye{ltP=*R2iqQ^cmt$c&0>Dv|#hdCkm5i~_ZDV*9* z-h42SF_O0mldWcDUa3yw*s$TxNmxkjcTZ%BU&i||JEAU|ny?k1K$yYXs!wbBs8mFR zOmrs)^C8$t=Y+nga8bcXpca3k>=!qj$pJ^p?UxhXLaq+!4h80k~Y6RY}gQM zRmkcjo)xF0lNsbyIB?Sv{%-7TmQ|sw`hYj~3ERO%fybOl@5n7Ud|N#YC}AwUx0M#_ zmk*Y&_vs~@o-ThSA2J1-$+R`%cK8va>-%ovO+sic*aX)pTl1JMqHi=0gJ7F((REC{ zBQsywsmMhW3J^?Oy#Gl>#sL0mmu)B*7E7l-npOQ%hlAGJeyn*f5hK~7p711oP?Kp|9t&u5FW>ETqz`eB&cEg- zyo$@+T4l18${ES=ppffXsnK1z#^%RMvUG~!MJEE^#A=3WFcr@ZFD1nzx)rvtyS$~$ z$yNpGV1k^GTJO|PB+&;PcD^L4TfoD++_(;sguZ^ zvS$=%Lr?nEhd|`roEx(-&7qH?oUEHtDsx#rwoNWHnk5&XK?tQ=<7juS+0u@6Js2?) zMRRvHk65J_sF6I|p^h;Y_H`7$Fs~_+r+7cs%Y4sEG7_-R_l&sh(sHB!%@8dWwIH<* zi@lA-+ZW5kR09OSw3_*3v@l^keYNOUrBzvzAF(hMc^P1#Qd6Ss*FZ2QkYW!pgI;XH z%R-){U~y`s4otk=+ob7{dYj1wmZ_-5Ec1!J9y4QamKU=;Vu?|4-P*&^&W*bEWYy2s zCiJ>ikEOI?T4wcAUhmA_(4NunwMn?|W?cOu^-H_IU`IU|F}IsnS2pf@=zm<3N+llx zn--VKvjv2-J3D3sLXXXEC2}e@V;RoC(s5rupY8r} zP$%rLwEnGPPPv&ayS?)pZhqww^Rf}|ND`sVCfhq(JrGvRd{G8dyeOwXh$)L%DjEMG zzJa$JlqE}vF|L;U$=Cft!*Oh*cJXD7a?9FyI3+NIuZC|+-ju0Zh1-{Zjoas3MdEQM zwP&`lkc{@53ALwEy~g~}eJmvTN$m?FZVLEK{VXJGZVadWW-f2V1Y20#7<_}p$eVgu zNS<2^LSNG-wL?sAz2R_UNPT*EXq~KDGF*vT0o^1VP`n~q;fW8=!m1yO>2e%D3QD>< z0zlBChAbOI;k+qlmPsh3_5+2`!g(EcOUyY3Ka~gck(wXxBeljHUFMs9O9QwbP4jef zzhKTLjf%i#GC1kp@!n+}{Zzxj~3Glw# zZLfDq2KwExl>hSlJ(0lNmNgUSuNR56E}M!UriPEv#c0md<<7rdQdl5a zyPX&Aq81tKH06ED-CM6oHmhP zp!F>-TN(WVRonqqhAI}%r+exHpo%1oMzgux?T(sFA#~_rwOlF z+{mn!OhY+e@2$J2*Liu zV(?||0Wg@*KNb-PQRLqeprS%T(7$0qFcCQHAK0JdzhiJwF`<88Na$a=!2g>H0Tupx zV<7|-{*T585!gQpB9J0~9|npL7K8t-F%$uV{<(sG$zNNID`(lC%N5sZN|jGyUwmb>YfjkwY+` z(91$_-Aw$fta_S#Xucge--=m%=Mml9{mYAVcXna$CV^Hyd};8uS_5ZJdzOV8Zt`fu z&2d&HlL>J9-hBkcxpW61=A>DJB;=SrnemVUrI}A@KlFs9-P58W@9jI8qjVl>+PLPF zPUmWCiIV-q@})SbHJ&B>om}ih)x%)VI-Li)Xp(I1e7A&;57+=_O>{rmn37Vn7@3W~`H0Z8ZtoLb z>l+z!E>>n#5^A@ro*do3rTw;|$+~8YQnGegt9vXkz994p4SmuRBt1$bpoUJ2SqMb__2)#LH@a)G$s4XO@2Gb+i3=Ps7MboW@U7?a&?d{LBW=KvoCbK$*S-&Q zNB^^vzv^2#Kon?^&;EAJ^h^9Zd6zIY@0}R_E|NM(cuHjpfNncJ9-du@9XuT&q1P0JodSJx_dlv*ByO)BS1~yWk^YhNFu>_1dOv8-C7e z6NTN!cSYsUbdvQt2)*3?xmp3>=6iFO4Z3;jEgC<+b@|oi|E*=pV_8pe@KtZf^x~om z{O7@yZO!%zShs8B>pFz)SDHVqz8$#rZ#Pz}Ucc`vausTQ*qv@d@gd$k3~;ej?T?t! zcO)y<5A@zb%FWks9MAmiSHbX`|= znu;IH7R2ViC_c0t+zt$x_gBDJT1t(DO4G7)F z(Sg~rAfNi+2nq?Kye-T7n%-ZT6T-1vAT=V3?>rcGoe1`BmKVAORUuca&+mCg5B`Fd zRM%j^z6+R<0(FVw)k|P%iBU6Tu8`Wt;5vux*1|}&lmg%c4k~m9afv)Zb*>ju7btH+ zElIo5R9E|6>+6EAim{ljeR^TP|BDj*CS83UN(iKo4dNjK^IszI*#{*>M#EZc7`p<0 zP<>7*CoI^m8z~ipDM^-aWFXo7 z&&F>&q5;@t5(QN*Mg?l7P_u9llDfmg*Vp1klWh?-uo}E!n;I=p;5;r~ca{ACKikYT z|GPOi`Y-NnBRAPxgBA(bX<89yU+2pba?i=wAgf%n!YkBi?}zRa=3DO2%)@l%HIZ0; z<#2LJb9uf^v{rWjhdGt0AQ`+4*`3sA+q>O-6{rj8lO0`DR4=}MChonS_KDs(dup^DH zkjfv@X#=Ge`Yh=2SYb#yOWr=t;)F>r{53?`R4}e&VeV<|)KlW!7?cWVDA~#nLkc_+ zXxCe|^ebyUOkb9oA%)e3%w*$1k6Bl;TvZqT3x4 z^IkKt0_(77Pn+Mub5E>srB=xL0LVN6@73!i6WG`IWmI);H@dsnAClM@z^D%vOL_T;fu{y#Mq>putFuFoS}5gaoy0#v(<5 zq~m9pdMtP7s^Lg|-|Fu;u!v=D;N*Op{ni~n$?8LuZ)XO~?I#gVqvODA4?OAqOZpi7 zUS@~tQSqiNUkB^77URnTD{=7YSU86tze9DpuzuL@?4-yh&c~anKtn$qWoAlf!rJfiB9%S`#1!NL#koxEPKD%H*tcgD8Z{UJqO?Zo#}Z1%l(Y(+>`*@0N%nl#0E-2vJgTgyJ~l56lI{7$*KRtt z?{2rPzF}h~)<}d;Z8EOU^rD@&VfXJ4{=B=wQF}P$q$P3e#Fs_uE5AnNsjW_nbqVe; zuAUS~PFuV=mGaB!71L?AnLhqN-wc0OtBP5xeBV~M#R?EAosMfH4XjPG{r)~H1z%-S zbef)&ifft(F?zoqU$7SHC97A{ynjyLR1NqaO0P?eewUKBI9Y~}=>~lEB=VdsN29Fe z&NpXY&~nVT<1aPh{ml3T-a}B1bFr%JW)%? z^{rVPD}Q;8RD$$UW9n;S<5EP{)h|^K@syl}d}-HZV8<8LO>Q_F{m|$9G1WRM?k5-t z;##@Lt`J<-m9M~~=?@A6QRd8NSxW0Jh8iv;%C=*+Jcec0W&q6hT;9hqW=`q069z43 z6;k`d1TU^7P?cR{dYlQ+vJ=fnWRCgx(xP?eo7qaD$+Xx3_pg^$W+EnUfM4(AiW1%s z6RpE-wiTx7_5cNr*FO7`tT+zh_a{=_)*JnT>2RH~kT0M0cASrUAZ~uG)RwnQ1gA4h z1;5->ItSRZV&ZHo`&f!x3MYQ{y;8g{hkDJA=TDtWwTd5YH_`c6VE--O<0YT7 zJ@I%bpI?L^UlN>VF`t9jkt}{G=+1E?laQd>rIrXqFh09xF+hl|&W;JCY1+cvFliqz4UCnjwPlZD-+ncRLsP2@Ho37fA{@h8jxLd+33bWO1iH2v39j= zV-$3zqp!Zu^kSn-j+xPwd10@hqqiSziid=aV|U}MKI)qL-u5bAF7+R4zMM$OvKxFw zLVe{1M^rgbF2J7rgWhIN*dL-eF%bJn(+U)qQY)cIxX4CZ)o~ac7mxk^Nt@T>!{IJ; z91^ma1qW-s;!R|1dZ+P8@kowUF>mRkszn70Uj~$yJHD@FsM|By?gOH#x}s7_CYt++ zcJS~eO=JJLLst&F)CaOucT@B*Q!@OB=2tDq1|>rvgU9ew?H#L60`)D1cE#W|Z0~A| ziXdgpEN&)YJcWT)=Gu`G)5accv&+$AkH`k^J4GrYL6aXgUjf)-u~!kpV7c zBQC8Qef>&F=6Ae>?4G?~F=&uq^w8&?{n`>6tE zl+#l#y~V85YJG1h-Xlj#dWfIH{(FBGF-WMna@-(0acv{>ZuZ4dx&F?&%V_eAWxr@Q zBa`d_)k@|X4VJAuf=$?k6N84u!Sl+WJVBBsJ#J7o=_%=NC_A5$tLDdVfE3`PJ5xax z!LQ#yy`{ZY(#zM51+3*-PdunOMIJ^99X9K_{ ztVZZfh*~H=l%&m*)PgPnxKwi;N@C5B5PI{HT8Bvl2J+cqfr27}yFds^qLrQhK?I|6 zIfSN+rO>Njr7XVy8=*fg(=^m}U9+q0GN%oIR8qoM9e;4!m{xvLX2QiYQAy9RF^$ZB z|GbS$Q`z!Kl3qJKXQ9>jCe919In>+d<(y@AtCZC{e_5lgy@CD|IGwi4diI8Ju6nzq zXmjD=&C{>VT=KvAHXpN=wpSmT7w9yB#kG?=m1qn*KYyY4*4zm@8me8+3UO-xc0^33|q1$j6n@ss;;5UKpC3wdc$|6$T46jt})K?YsCrjJH|H zAy#9j`w3~I$#^%{kY`ITm&UTWPH3;LjZvjOWnv>0lSCp@aZKt!lA?2$D3!{gTx5oL zTxY7nLr%TtIYJ)=bOJ$y41-(Ttlg%|2S!-YNQi6F8L(=8GKKfK4FB>lp@~2xmt==s zddpG9Or{lEF+XtUKy-XIX^S^G%t+F-;@6h*)45qwrEh7EEMf#Rdr?1UwG9g{wZE2|O39+~zbQkOCpUzUa7i(^C`w#h zP!!HDDk|{j^K(^V++mP^3sE>6DUteA zg$IDc#l`E5wjU"data/graphdata_exponent_mix32.m"]; +gsraw=Import[NotebookDirectory[]<>"data/graphdata_exponent_hightau.m"]; gsraw=SortBy[gsraw,#[[1,1]]&]; (* Sort by n *) averagesGrouped=GatherBy[gsraw,{#[[1,2]]&,#[[1,1]]&}]; +gsraw2=Import[NotebookDirectory[]<>"data/graphdata_canonical_properties2.m"]; +gsraw2=SortBy[gsraw2,#[[1,1]]&]; (* Sort by n *) +averagesGrouped2=GatherBy[gsraw2,{#[[1,2]]&,#[[1,1]]&}]; +canonicalDatapoints=Map[{#[[1,1,1]],Mean[#[[All,2]]]}&,averagesGrouped2,{2}]; + + (* 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]]]; @@ -61,6 +67,10 @@ mediansFits=Map[Fit[#,{1,logn},logn]&,mediansLoglogdata]; mediansFitsExtra=Map[LinearModelFit[#,logn,logn]&,mediansLoglogdata]; +canonicalloglogdata=Log[canonicalDatapoints[[All,nRange]]]; +canonicalFitsExtra=Map[LinearModelFit[#,logn,logn]&,canonicalloglogdata]; + + averagesFitsExtra[[1]]["ParameterConfidenceIntervalTable"] averagesFitsExtra[[1]]["BestFitParameters"] averagesFitsExtra[[1]]["ParameterErrors"] @@ -91,10 +101,18 @@ Show[ListPlot[{averagesExponents,mediansExponents},Joined->True,PlotMarkers->Aut (* For visual, shift the tau values slightly left or right to distinguish the two datasets *) tauValues=averagesGrouped[[All,1,1,1,2]]; averagesExponentsErrorBars=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&, -Transpose[{tauValues-0.001,averagesFitsExtra}]]; +Transpose[{tauValues-0.003,averagesFitsExtra}]]; mediansExponentsErrorBars=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&, -Transpose[{tauValues+0.001,mediansFitsExtra}]]; -plot2=Show[ErrorListPlot[{averagesExponentsErrorBars,mediansExponentsErrorBars},Joined->True,PlotMarkers->Automatic,Frame->True,FrameLabel->{"tau","triangle exponent"},PlotRange->{{2,3},{0,1.6}},ImageSize->300],Plot[3/2(3-tau),{tau,2,3},PlotStyle->{Dashed}]] +Transpose[{tauValues+0.003,mediansFitsExtra}]]; +canonicalExponentsErrorBars=Map[{{#[[1]],#[[2]]["BestFitParameters"][[2]]},ErrorBar[#[[2]]["ParameterConfidenceIntervals"][[2]]-#[[2]]["BestFitParameters"][[2]]]}&, +Transpose[{tauValues+0.000,canonicalFitsExtra}]]; +plot2=Show[ +ErrorListPlot[{averagesExponentsErrorBars,mediansExponentsErrorBars,canonicalExponentsErrorBars}, +Joined->True,PlotMarkers->Automatic, +Frame->True,FrameLabel->{"tau","triangle exponent"}, +PlotRange->{{2,3},{0,1.6}}, +ImageSize->300], +Plot[3/2(3-tau),{tau,2,3},PlotStyle->{Black,Dashed}]] Export[NotebookDirectory[]<>"plots/triangle_exponent.pdf",plot2]