AENC/switchchain Changeset - d0883e1df741 · Centrum Wiskunde & Informatica (CWI)

Changeset - d0883e1df741

Parent rev.

Child rev.

[Not reviewed]

0 6 0

Tom Bannink - 8 years ago 2017-08-17 16:59:00
tom.bannink@cwi.nl

Add proper cputime plots

6 files changed with 125 insertions and 44 deletions:

cpp/graph.hpp

cpp/graph_ccm.hpp

cpp/switchchain_ccm_cputime.cpp

plots/timeevol_ccm.pdf

bin+mod

plots/timeevol_ccm_log.pdf

bin+mod

triangle_ccm_timeevol_plots.m

0 comments (0 inline, 0 general)

cpp/graph.hpp

➞

Show inline comments

@@ @@ -23,41 +23,44 @@ class StoredEdge { @@
 class DiDegree {
   public:
     unsigned int in;
     unsigned int out;
 };
 typedef std::vector<unsigned int> DegreeSequence;
 typedef std::vector<DiDegree> DiDegreeSequence;
 template< class RandomIt, class Compare >
 void insertionSort( RandomIt first, RandomIt last, Compare comp ) {
     if (first == last)
         return;
     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
         //             next
         // 1 2 3 4 5 6  4
         //           a  b
         // 1 2 3 4 5 x  6
         //         a b
         // 1 2 3 4 x 5  6
         //       a b
         // 1 2 3 4 4 5  6
         while (b != first && comp(newvalue, *a)) { // if newvalue < *a
             *b = *a;
             a--;
             b--;
             b = a--;
+        }
         *b = newvalue;
+    }
+}
 class Graph {
   public:
     Graph() {}
     Graph(unsigned int n) { reset(n); }
     ~Graph() {}
@@ @@ -77,62 +80,81 @@ class Graph { @@
+    }
     unsigned int edgeCount() const { return edges.size(); }
     const Edge &getEdge(unsigned int i) const { return edges[i].e; }
     const auto& getAdj() const { return adj; }
     const auto& getBooleanAdj() const { return badj; }
     // When the degree sequence is not graphics, the Graph can be
     // in any state, it is not neccesarily empty
     bool createFromDegreeSequence(const DegreeSequence &d) {
+    bool createFromDegreeSequence(const DegreeSequence &d, bool slowSort = false) {
         // Havel-Hakimi algorithm
         // Based on Erdos-Gallai theorem
         unsigned int n = d.size();
         // degree, vertex index
         std::vector<std::pair<unsigned int, unsigned int>> degrees(n);
         for (unsigned int i = 0; i < n; ++i) {
             degrees[i].first = d[i];
             degrees[i].second = i;
+        }
         // Clear the graph
         reset(n);
         // First use general sorting algorithm
         std::sort(
             degrees.begin(), degrees.end(),
             [](const auto &p1, const auto &p2) { return p1.first < p2.first; });
         while (!degrees.empty()) {
             insertionSort(degrees.begin(), degrees.end(),
                           [](const auto& p1, const auto& p2) {
                               return p1.first < p2.first;
                           });
             // Highest degree is at back of the vector
             // Take it out
             unsigned int degree = degrees.back().first;
             unsigned int u = degrees.back().second;
             degrees.pop_back();
             // In the end, all remaining degrees are zero
             if (degree == 0) {
                 break;
+            }
             if (degree > degrees.size()) {
                 return false;
+            }
             // Now loop over the last 'degree' entries of degrees
             auto rit = degrees.rbegin();
             for (unsigned int i = 0; i < degree; ++i) {
                 if (rit->first == 0 || !addEdge({u, rit->second})) {
                     return false;
+                }
                 rit->first--;
                 ++rit;
+            }
             // To compare influence of sorting methods
             if (slowSort) {
                 // Re-sort using general sort
                 std::sort(degrees.begin(), degrees.end(),
                               [](const auto &p1, const auto &p2) {
                                   return p1.first < p2.first;
                               });
             } else {
                 // Re-sort using insertion sort
                 insertionSort(degrees.begin(), degrees.end(),
                               [](const auto &p1, const auto &p2) {
                                   return p1.first < p2.first;
                               });
+            }
+        }
         return true;
+    }
     DegreeSequence getDegreeSequence() const {
         DegreeSequence d(adj.size());
         std::transform(adj.begin(), adj.end(), d.begin(),
                        [](const auto &u) { return u.size(); });
         return d;
+    }
     // Assumes valid vertex indices

cpp/graph_ccm.hpp

➞

Show inline comments

                                    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();
     // degree, vertex index
     std::vector<std::pair<unsigned int, unsigned int>> degrees(n);
     for (unsigned int i = 0; i < n; ++i) {
         degrees[i].first = ds[i];
         degrees[i].second = i;
+    }
     // This should make the random pairing faster
     // More likely to pair up with something at the end of an array
     // So erasing that element from the array then requires less moves
     std::sort(
             degrees.begin(), degrees.end(),
             [](const auto &p1, const auto &p2) { return p1.first < p2.first; });
     // remaining half-edges , iterator into `degrees`
     std::vector<std::pair<unsigned int, decltype(degrees.begin())>> available;
     available.reserve(n);
     // Clear the graph
     g.reset(n);
     while (!degrees.empty()) {
         // Get the highest degree:
         // If there are multiple highest ones, we pick the first one
         unsigned int dmax = 0;
         auto uIter = degrees.begin();
             return false;
         if (method2) {
             // Now assign randomly to the remaining vertices
             // weighted by the number of half-edges they have left
             while (dmax--) {
                 std::uniform_int_distribution<> distr(1, availableEdges);
                 // Go from edge index to number. We should really have a proper
                 // data structure like some cumulative tree
                 unsigned int halfEdge = distr(rng);
                 unsigned int cumulative = 0;
                 auto vIter = uIter;
                 for (auto iter = available.begin(); iter != available.end();
                 // Reverse has higher probability to be fast because high
                 // degrees are near the end
                 for (auto iter = available.rbegin(); iter != available.rend();
                      ++iter) {
                     cumulative += iter->first;
                     if (halfEdge <= cumulative) {
                         vIter = iter->second;
                         availableEdges -= iter->first;
                         available.erase(iter);
                         // Reverse iterators are tricky
                         available.erase(std::next(iter).base());
                         break;
+                    }
+                }
                 if (vIter == uIter)
                     std::cerr << "ERROR 1 in CCM2.\n";
                 if (vIter->first == 0)
                     std::cerr << "ERROR 2 in CCM2.\n";
                 if (!g.addEdge({u, vIter->second}))
                     std::cerr << "ERROR 3 in CCM2.\n";
                 uIter->first--;
                 vIter->first--;
+            }
             for (auto iter = degrees.begin(); iter != degrees.end();)
                 if (iter->first == 0)
                     iter = degrees.erase(iter);
                 else
                     ++iter;
         } else {
             std::uniform_int_distribution<> distr(1, availableEdges);
             // Go from edge index to number. We should really have a proper
             // data structure like some cumulative tree
             auto vIter = uIter;
             unsigned int halfEdge = distr(rng);
             unsigned int cumulative = 0;
             for (auto iter = available.begin(); iter != available.end();
             // Reverse has higher probability to be fast because high degrees
             // are near the end
             for (auto iter = available.rbegin(); iter != available.rend();
                  ++iter) {
                 cumulative += iter->first;
                 if (halfEdge <= cumulative) {
                     vIter = iter->second;
                     break;
+                }
+            }
             if (vIter == uIter)
                 std::cerr << "ERROR 4 in CCM2.\n";
             // pair u to v
             if (vIter->first == 0)
                 std::cerr << "ERROR 5 in CCM1.\n";

cpp/switchchain_ccm_cputime.cpp

➞

Show inline comments

@@ @@ -4,30 +4,30 @@ @@
 #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 numVerticesMin = 10000;
     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};
-    float tauValues[] = {2.1f, 2.3f, 2.5f, 2.7f, 2.9f};
     float tauValues[] = {2.1f, 2.5f, 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]);
@@ @@ -44,24 +44,25 @@ int main(int argc, char* argv[]) { @@
             << " 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 << "6: slow-sort-HH 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) {
@@ @@ -70,35 +71,34 @@ int main(int argc, char* argv[]) { @@
             // 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);
                     if (!chain.initialize(g, true)) {
                         std::cerr << "Could not initialize Markov chain.\n";
                         return 1;
+                    }
                     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)
@@ @@ -130,23 +130,43 @@ int main(int argc, char* argv[]) { @@
                                     dt.count(), chain.g.getTrackedTriangles());
                                 chain.doMove(true);
+                            }
                             outfile << ',' << '{' << i << ',' << triangleSeq << '}';
                             failed = false;
                             break;
+                        }
+                    }
                     if (failed)
                         outfile << ",{1000,{}}";
+                }
                 // Slow sort method
+                {
                     // record start time
                     auto start = std::chrono::high_resolution_clock::now();
                     // Incorporate Havel-Hakimi time
                     g.createFromDegreeSequence(ds, true);
                     if (!chain.initialize(g, true)) {
                         std::cerr << "Could not initialize Markov chain.\n";
                         return 1;
+                    }
                     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 << ',' << triangleSeq;
                 outfile << '}' << std::flush;
                 std::cout << " Finished timed CCM time evols." << std::flush;
                 std::cout << std::endl;
+            }
+        }
+    }
     outfile << '\n' << '}';
     return 0;
+}

plots/timeevol_ccm.pdf

➞

Show inline comments

binary diff not shown

plots/timeevol_ccm_log.pdf

➞

Show inline comments

binary diff not shown

triangle_ccm_timeevol_plots.m

➞

Show inline comments

 plot1=ListPlot[dataSetsFlattened,Joined->True,PlotRange->{{0,40000},{0,6000}},PlotStyle->colorList,ImageSize->300,PlotLabel->nlabels[[1]],Frame->True,FrameLabel->{"timesteps","number of triangles"}]
 plot1log=ListLogPlot[dataSetsFlattened,Joined->True,PlotRange->{{0,40000},{1,15000}},PlotStyle->colorList,ImageSize->300,PlotLabel->nlabels[[1]],Frame->True,FrameLabel->{"timesteps","number of triangles"}]
 Export[NotebookDirectory[]<>"plots/timeevol_ccm.pdf",plot1]
 Export[NotebookDirectory[]<>"plots/timeevol_ccm_log.pdf",plot1log]
 (* ::Subsection:: *)
 (*New code with CPU time*)
 getCombinedTimeData[run_]:=Module[{maxTime,skipPts,hhData,ccm1Data,ccm2Data},
+getCombinedTimeData[run_]:=Module[{maxTime,skipPts,avg,hhData,ccm1Data,ccm2Data,hhSlowData,line1,line2,line3,line4,lineAvg},
 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}
 avg=Mean[run[[3,Floor[(1/2)*Length[run[[3]]]];;-1,2]]];
 hhData=  run[[3,1;;-1;;skipPts]];
 (*ccm1Data=run[[4,2,1;;-1;;skipPts]];*)
 ccm2Data=run[[5,2,1;;-1;;skipPts]];
 hhSlowData=run[[6,1;;-1;;skipPts]];
 line1={{0,run[[3,1,2]]},run[[3,1]]};
 line2={{0,run[[4,2,1,2]]},run[[4,2,1]]};
 line3={{0,run[[5,2,1,2]]},run[[5,2,1]]};
 line4={{0,run[[6,1,2]]},run[[6,1]]};
 lineAvg={{0,avg},{3*run[[3,-1,1]],avg}};
 {Legended[hhData,"\[Tau] = "<>ToString[run[[1,2]]]],
 (*ccm1Data,*)
 ccm2Data,
 hhSlowData,
 line1,
 (*line2,*)
 line3,
 line4,
 lineAvg}
+]
 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]
 getStartPoints[run_]:={run[[3,1]],run[[5,2,1]],run[[6,1]]};
 getLogScaleStartPoints[run_]:={{run[[3,1,1]],Log[run[[3,1,2]]]},{run[[5,2,1,1]],Log[run[[5,2,1,2]]]},{run[[6,1,1]],Log[run[[6,1,2]]]}};
 z2
 plot1=ListPlot[dataSetsFlattened,Joined->True,PlotRange->{All,All},PlotStyle->colorList,ImageSize->300,PlotLabel->nlabels[[1]],Frame->True,FrameLabel->{"seconds","number of triangles"}]
 dataSets=Map[getCombinedTimeData,gdata,{3}];
 startPoints=Map[getStartPoints,gdata,{3}];
 logScaleStartPoints=Map[getLogScaleStartPoints,gdata,{3}];
 dataSetsFlattened=Flatten[dataSets,2];
 styleList=Table[{
 Directive[Dashing[None],ColorData[97,"ColorList"][[i]]],
 Directive[Dashing[None],ColorData[97,"ColorList"][[i]]],
 Directive[Dashing[None],ColorData[97,"ColorList"][[i]]],
 Directive[Dashed,ColorData[97,"ColorList"][[i]]],
 Directive[Dashed,ColorData[97,"ColorList"][[i]]],
 Directive[Dashed,ColorData[97,"ColorList"][[i]]],
 Directive[Dashing[None],Black,Thin]
 },{i,1,Length[dataSetsFlattened]}];
 dataSetsFlattened=Flatten[dataSetsFlattened,1];
 styleList=Flatten[styleList,1];
 pointList={Red,PointSize[Medium],Point[Flatten[startPoints,3]]};
 logPointList={Red,PointSize[Medium],Point[Flatten[logScaleStartPoints,3]]};
 plot1=ListPlot[dataSetsFlattened,Joined->True,PlotRange->{{0,0.55},All},PlotStyle->styleList,Epilog->pointList,ImageSize->300,PlotLabel->nlabels[[1]],Frame->True,FrameLabel->{"seconds","number of triangles"}]
 plot1log=ListLogPlot[dataSetsFlattened,Joined->True,PlotRange->{{0,0.55},All},PlotStyle->styleList,Epilog->logPointList,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]

0 comments (0 inline, 0 general)