Changeset - 9df034849ada
Parent rev.
Child rev.
[Not reviewed]
0 2 0
Tom Bannink - 8 years ago 2017-04-06 17:09:45
tom.bannink@cwi.nl
Write some conclusions to top of notebook
2 files changed with 140 insertions and 30 deletions:
cpp/switchchain.cpp
108
19
showgraphs.m
32
11
0 comments (0 inline, 0 general)
cpp/switchchain.cpp
➞
Show inline comments
 
@@ -78,6 +78,9 @@ bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, auto& rng, bool meth
 
        degrees[i].second = i;
 
    }
 

			




	
	
	
		
 
    std::vector<decltype(degrees.begin())> available;

			




	
	
	
		
 
    available.reserve(n);

			




	
	
	
		
 

			




	
	
	
		
 
    // Clear the graph

			




	
	
	
		
 
    g.reset(n);

			




	
	
	
		
 

			




	
	
		
 
@@ -99,6 +102,10 @@ bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, auto& rng, bool meth

			




	
	
	
		
 
        if (dmax > degrees.size() - 1)

			




	
	
	
		
 
            return false;

			




	
	
	
		
 

			




	
	
	
		
 
        if (dmax == 0) {

			




	
	
	
		
 
            std::cerr << "ERROR 1 in GCM.\n";

			




	
	
	
		
 
        }

			




	
	
	
		
 

			




	
	
	
		
 
        unsigned int u = uIter->second;

			




	
	
	
		
 

			




	
	
	
		
 
        if (method2) {

			




	
	
		
 
@@ -120,21 +127,28 @@ bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, auto& rng, bool meth

			




	
	
	
		
 
                    vIter++;

			




	
	
	
		
 
            }

			




	
	
	
		
 
        } else {

			




	
	
	
		
 
            // Pair with a random vertex that is not u itself

			




	
	
	
		
 
            std::uniform_int_distribution<> distr(0, degrees.size() - 2);

			




	
	
	
		
 
            auto vIter = degrees.begin() + distr(rng);

			




	
	
	
		
 
            if (vIter == uIter)

			




	
	
	
		
 
                vIter++;

			




	
	
	
		
 
            // Pair with a random vertex that is not u itself and to which

			




	
	
	
		
 
            // u has not paired yet!

			




	
	
	
		
 
            available.clear();

			




	
	
	
		
 
            for (auto iter = degrees.begin(); iter != degrees.end(); ++iter) {

			




	
	
	
		
 
                if (iter->second != u && !g.hasEdge({u, iter->second}))

			




	
	
	
		
 
                    available.push_back(iter);

			




	
	
	
		
 
            }

			




	
	
	
		
 
            if (available.empty())

			




	
	
	
		
 
                return false;

			




	
	
	
		
 
            std::uniform_int_distribution<> distr(0, available.size() - 1);

			




	
	
	
		
 
            auto vIter = available[distr(rng)];

			




	
	
	
		
 
            // pair u to v

			




	
	
	
		
 
            if (vIter->first == 0)

			




	
	
	
		
 
                std::cerr << "ERROR in GCM1.\n";

			




	
	
	
		
 
            if (!g.addEdge({uIter->second, vIter->second}))

			




	
	
	
		
 
                std::cerr << "ERROR 2 in GCM1.\n";

			




	
	
	
		
 
            if (!g.addEdge({u, vIter->second}))

			




	
	
	
		
 
                std::cerr << "ERROR. Could not add edge in GCM1.\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

			




	
	
	
		
 
            if (dmax == 1) {

			




	
	
	
		
 
                // Remove both

			




	
	
	
		
 
                // Erasure invalidates all iterators AFTER the erased one

			




	
	
	
		
 
                if (vIter > uIter) {

			




	
	
	
		
 
                    degrees.erase(vIter);

			




	
	
	
		
 
                    degrees.erase(uIter);

			




	
	
		
 
@@ -144,6 +158,7 @@ bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, auto& rng, bool meth

			




	
	
	
		
 
                }

			




	
	
	
		
 
            } else {

			




	
	
	
		
 
                // Remove only v if it reaches zero

			




	
	
	
		
 
                uIter->first--;

			




	
	
	
		
 
                vIter->first--;

			




	
	
	
		
 
                if (vIter->first == 0)

			




	
	
	
		
 
                    degrees.erase(vIter);

			




	
	
		
 
@@ -167,6 +182,7 @@ int main() {

			




	
	
	
		
 
    float tauValues[] = {2.1f, 2.2f, 2.3f, 2.4f, 2.5f, 2.6f, 2.7f, 2.8f};

			




	
	
	
		
 

			




	
	
	
		
 
    Graph g;

			




	
	
	
		
 
    Graph g1;

			




	
	
	
		
 
    Graph g2;

			




	
	
	
		
 

			




	
	
	
		
 
    std::ofstream outfile("graphdata.m");

			




	
	
		
 
@@ -181,8 +197,9 @@ int main() {

			




	
	
	
		
 
            //std::poisson_distribution<> degDist(12);

			




	
	
	
		
 

			




	
	
	
		
 
            // For a single n,tau take samples over several instances of

			




	
	
	
		
 
            // the degree distribution

			




	
	
	
		
 
            for (int degreeSample = 0; degreeSample < 500; ++degreeSample) {

			




	
	
	
		
 
            // the degree distribution.

			




	
	
	
		
 
            // 500 samples seems to give reasonable results

			




	
	
	
		
 
            for (int degreeSample = 0; degreeSample < 200; ++degreeSample) {

			




	
	
	
		
 
                // Generate a graph

			




	
	
	
		
 
                // might require multiple tries

			




	
	
	
		
 
                for (int i = 1; ; ++i) {

			




	
	
		
 
@@ -207,14 +224,23 @@ int main() {

			




	
	
	
		
 
                }

			




	
	
	
		
 

			




	
	
	
		
 
                //

			




	
	
	
		
 
                // Test the greedy configuration model success rate

			




	
	
	
		
 
                // Test the GCM1 and GCM2 success rate

			




	
	
	
		
 
                //

			




	
	
	
		
 
                std::vector<int> greedyTriangles;

			




	
	
	
		
 
                int successrate = 0;

			




	
	
	
		
 
                std::vector<int> greedyTriangles1;

			




	
	
	
		
 
                std::vector<int> greedyTriangles2;

			




	
	
	
		
 
                int successrate1 = 0;

			




	
	
	
		
 
                int successrate2 = 0;

			




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

			




	
	
	
		
 
                    if (greedyConfigurationModel(ds, g2, rng, true)) {

			




	
	
	
		
 
                        ++successrate;

			




	
	
	
		
 
                        greedyTriangles.push_back(g2.countTriangles());

			




	
	
	
		
 
                    Graph gtemp;

			




	
	
	
		
 
                    if (greedyConfigurationModel(ds, gtemp, rng, false)) {

			




	
	
	
		
 
                        ++successrate1;

			




	
	
	
		
 
                        greedyTriangles1.push_back(gtemp.countTriangles());

			




	
	
	
		
 
                        g1 = gtemp;

			




	
	
	
		
 
                    }

			




	
	
	
		
 
                    if (greedyConfigurationModel(ds, gtemp, rng, true)) {

			




	
	
	
		
 
                        ++successrate2;

			




	
	
	
		
 
                        greedyTriangles2.push_back(gtemp.countTriangles());

			




	
	
	
		
 
                        g2 = gtemp;

			




	
	
	
		
 
                    }

			




	
	
	
		
 
                }

			




	
	
	
		
 

			




	
	
		
 
@@ -224,6 +250,17 @@ int main() {

			




	
	
	
		
 
                    return 1;

			




	
	
	
		
 
                }

			




	
	
	
		
 

			




	
	
	
		
 
                SwitchChain chain1, chain2;

			




	
	
	
		
 
                bool do1 = true, do2 = true;

			




	
	
	
		
 
                if (!chain1.initialize(g1)) {

			




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

			




	
	
	
		
 
                    do1 = false;

			




	
	
	
		
 
                }

			




	
	
	
		
 
                if (!chain2.initialize(g2)) {

			




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

			




	
	
	
		
 
                    do2 = false;

			




	
	
	
		
 
                }

			




	
	
	
		
 

			




	
	
	
		
 
                std::cout << "Running n = " << numVertices << ", tau = " << tau

			




	
	
	
		
 
                          << ". \t" << std::flush;

			




	
	
	
		
 

			




	
	
		
 
@@ -232,7 +269,7 @@ int main() {

			




	
	
	
		
 
                //constexpr int measureSkip =

			




	
	
	
		
 
                //    200; // Take a sample every ... steps

			




	
	
	
		
 
                int mixingTime = 0;

			




	
	
	
		
 
                constexpr int measurements = 10000;

			




	
	
	
		
 
                constexpr int measurements = 5000;

			




	
	
	
		
 
                constexpr int measureSkip = 1;

			




	
	
	
		
 

			




	
	
	
		
 

			




	
	
		
 
@@ -255,7 +292,8 @@ int main() {

			




	
	
	
		
 
                          << " moves succeeded ("

			




	
	
	
		
 
                          << 100.0f * float(movesDone) /

			




	
	
	
		
 
                                 float(mixingTime + measurements * measureSkip)

			




	
	
	
		
 
                          << "%)." << std::endl;

			




	
	
	
		
 
                          << "%).";

			




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

			




	
	
	
		
 

			




	
	
	
		
 
                if (outputComma)

			




	
	
	
		
 
                    outfile << ',';

			




	
	
		
 
@@ -263,8 +301,59 @@ int main() {

			




	
	
	
		
 

			




	
	
	
		
 
                std::sort(ds.begin(), ds.end());

			




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

			




	
	
	
		
 
                outfile << ',' << triangles << ',' << ds;

			




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

			




	
	
	
		
 
                outfile << ',' << triangles;

			




	
	
	
		
 
                outfile << ',' << ds;

			




	
	
	
		
 
                outfile << ',' << greedyTriangles1;

			




	
	
	
		
 
                outfile << ',' << greedyTriangles2;

			




	
	
	
		
 

			




	
	
	
		
 
                if (do1) {

			




	
	
	
		
 
                    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;

			




	
	
	
		
 
                }

			




	
	
	
		
 
                if (do2) {

			




	
	
	
		
 
                    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;

			




	
	
	
		
 
            }

			




	
	
	
		
 
        }

			




	
	
	
		
 
    }

			




        

    


    
    

    

        

                

                    showgraphs.m
                

                

                  
                      
                        

                      ➞

                      
                        
                  

                  
                      
                  
                      
                  
                      
                  
                      
                  
                  
                

                

                    Show inline comments
                    
                

        

        

            



	
	
		
 
@@ -14,6 +14,8 @@ Needs["ErrorBarPlots`"]

			




	
	
	
		
 
(*   (b) How often does this procedure work/fail. Might still be faster to do switchings from Erdos-Gallai.*)

			




	
	
	
		
 
(*   (c) Compare two greedy ways: (c1) first take highest and finish all its pairings  (c2) take new highest after every single pairing*)

			




	
	
	
		
 
(*   (d) Time evolution for GCM on top of Erdos-Gallai time evolution.*)

			




	
	
	
		
 
(* The initial #triangles in GCM2 is somewhere between 0 and 5 standard deviations below the average #triangles, whereas the #triangles in Erdos-Gallai can be as high as 100 standard deviations above the average.*)

			




	
	
	
		
 
(* TODO: Not only compare number of standard deviations but also percentage above/below average.*)

			




	
	
	
		
 
(**)

			




	
	
	
		
 
(*- Count #moves that result in +-k triangles (one move could create many triangles at once!)*)

			




	
	
	
		
 
(**)

			




	
	
		
 
@@ -37,13 +39,20 @@ Needs["ErrorBarPlots`"]

			




	
	
	
		
 
(*- Compute  Sum over i<j<k  of  (1-Exp[- d_i d_j / (2E)]) * (1 - Exp[-d_j d_k / (2E)]) * (1 - Exp[-d_k d_i / (2E)]) .*)

			




	
	
	
		
 
(*  Computing the f(i,j) = (1-Exp[..]) terms is fine, but then computing Sum[ f(i,j) f(j,k) f(i,k) ) ] over n^3 entries is very slow.*)

			




	
	
	
		
 
(*  *)

			




	
	
	
		
 
(*  - Improve runtime*)

			




	
	
	
		
 
(*- Improve runtime*)

			




	
	
	
		
 
(*   (a) Don't remove/add edges from the std::vector. Simply replace them. Done, is way faster for large n.*)

			




	
	
	
		
 
(*   (b) Do not choose the three permutations with 1/3 probability: choose the "staying" one with zero probability. Should still be a valid switch chain?*)

			




	
	
	
		
 
(*   *)

			




	
	
	
		
 
(*   - Experimental mixing time as function of n. At (n,tau)=(1000,2.5) it seems to be between 10.000 and 20.000 steps.*)

			




	
	
	
		
 
(*     Done. Seems to be something like  (1/2)(32-26(tau-2))n  so we run it for that time without the factor (1/2).*)

			




	
	
	
		
 
(*  *)

			




	
	
	
		
 
(*- Experimental mixing time as function of n. At (n,tau)=(1000,2.5) it seems to be between 10.000 and 20.000 steps.*)

			




	
	
	
		
 
(*   Done. Seems to be something like  (1/2)(32-26(tau-2))n  so we run it for that time without the factor (1/2).*)

			




	
	
	
		
 
(**)

			




	
	
	
		
 
(*- GCM1: Greedy Configuration Model 1: take highest degree and completely do all its pairings (at random)*)

			




	
	
	
		
 
(*- GCM2: Greedy Configuration Model 2: take highest degree and do a single pairing, then take new highest degree. So this matters mostly if there are multiple high degree nodes*)

			




	
	
	
		
 
(*The success rates, conditioned on the degree sequence being graphical, is almost always higher using GCM2. For certain degree sequences the success rate of GCM2 can be 0.9 higher than that of GCM1. (i.e. amost always works vs almost always fails).*)

			




	
	
	
		
 
(*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%.*)

			




	
	
	
		
 
(**)

			




	
	
	
		
 
(**)

			




	
	
	
		
 
(*  *)

			




	
	
	
		
 

			




	
	
	
		
 

			




	
	
		
 
@@ -65,11 +74,7 @@ Grid[Partition[gs,10],Frame->All]

			




	
	
	
		
 

			




	
	
	
		
 

			




	
	
	
		
 
(* ::Section:: *)

			




	
	
	
		
 
(*Plot triangle counts*)

			




	
	
	
		
 

			




	
	
	
		
 

			




	
	
	
		
 
(* ::Subsection:: *)

			




	
	
	
		
 
(*Data import and data merge*)

			




	
	
	
		
 
(*Data import*)

			




	
	
	
		
 

			




	
	
	
		
 

			




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

			




	
	
		
 
@@ -82,11 +87,19 @@ nlabels=Map["n = "<>ToString[#]&,gdata[[1,All,1,1,1]]];

			




	
	
	
		
 
taulabels=Map["tau = "<>ToString[#]&,gdata[[All,1,1,1,2]]];

			




	
	
	
		
 

			




	
	
	
		
 

			




	
	
	
		
 
(* ::Subsection:: *)

			




	
	
	
		
 
(*Merge data*)

			




	
	
	
		
 

			




	
	
	
		
 

			




	
	
	
		
 
newData=Import[NotebookDirectory[]<>"data/graphdata_3.m"];

			




	
	
	
		
 
mergedData=Import[NotebookDirectory[]<>"data/graphdata_merged.m"];

			




	
	
	
		
 
Export[NotebookDirectory[]<>"data/graphdata_merged_new.m",Join[mergedData,newData]]

			




	
	
	
		
 

			




	
	
	
		
 

			




	
	
	
		
 
(* ::Section:: *)

			




	
	
	
		
 
(*Plot triangle counts*)

			




	
	
	
		
 

			




	
	
	
		
 

			




	
	
	
		
 
(* ::Subsection:: *)

			




	
	
	
		
 
(*Plot empirical distribution of maximum degree*)

			




	
	
	
		
 

			




	
	
		
 
@@ -241,7 +254,7 @@ getStats[run_]:=Module[{avg,stddev},

			




	
	
	
		
 
stats=Map[getStats,gdata,{3}];

			




	
	
	
		
 

			




	
	
	
		
 

			




	
	
	
		
 
histograms=Map[Histogram[#[[1,3]],PlotRange->{{-8,8},Automatic},PlotLabel->"ErdosGallai deviation: "<>ToString[#[[1,2]]]]&,stats,{2}];

			




	
	
	
		
 
histograms=Map[Histogram[#[[1,3]],PlotRange->{{-5,5},Automatic},PlotLabel->"ErdosGallai deviation: "<>ToString[#[[1,2]]]]&,stats,{2}];

			




	
	
	
		
 

			




	
	
	
		
 

			




	
	
	
		
 
TableForm[histograms,TableHeadings->{taulabels,nlabels}]

			




	
	
		
 
@@ -252,8 +265,16 @@ TableForm[histograms,TableHeadings->{taulabels,nlabels}]

			




	
	
	
		
 

			




	
	
	
		
 

			




	
	
	
		
 
(* gdata[[ tau index, n index, run index , {ntau, #tris, ds, greedyTriangles} ]] *)

			




	
	
	
		
 
successrates=Map[Length[#[[4]]]&,gdata,{3}];

			




	
	
	
		
 
successrates=Map[{Length[#[[4]]],Length[#[[5]]]}&,gdata,{3}];

			




	
	
	
		
 
successrates=Map[Transpose,successrates,{2}];

			




	
	
	
		
 
successratesDelta=Map[Length[#[[5]]]-Length[#[[4]]]&,gdata,{3}];

			




	
	
	
		
 

			




	
	
	
		
 
rateHistograms=Map[Histogram[#,{10},PlotRange->{{0,100},Automatic}]&,successrates,{2}];

			




	
	
	
		
 
TableForm[rateHistograms,TableHeadings->{taulabels,nlabels}]

			




	
	
	
		
 

			




	
	
	
		
 
rateHistograms=Map[Histogram[#,{10},PlotRange->{{-100,100},Automatic}]&,successratesDelta,{2}];

			




	
	
	
		
 
TableForm[rateHistograms,TableHeadings->{taulabels,nlabels}]

			




	
	
	
		
 
(*TableForm[Transpose[rateHistograms],TableHeadings->{nlabels,taulabels}]*)

			




	
	
	
		
 

			




	
	
	
		
 

			




	
	
	
		
 

			




        

    



    


    



    



      

      





    
    0 comments (0 inline, 0 general)
    





    


  

  







  


    




    








    
        Server instance: kalli-2997587
    
    
        This site is powered by
            Kallithea,
        which is
        © 2010–2022 by various authors & licensed under GPLv3.
            – Support