diff --git a/cpp/showgraphs.m b/cpp/showgraphs.m
index f9efe4e7088cfebf80fddea00414d73ed856b55b..beec01993c5f52bbb41157fe018ab21c4d0297b0 100644
--- a/cpp/showgraphs.m
+++ b/cpp/showgraphs.m
@@ -24,18 +24,73 @@ Grid[Partition[gs,10],Frame->All]
 (*Plot triangle counts*)
 
 
-gsraw=Import[NotebookDirectory[]<>"graphdata.m"];
+(* ::Subsection:: *)
+(*Data import and data merge*)
+
+
+gsraw=Import[NotebookDirectory[]<>"graphdata_merged.m"];
+
+
+newData=Import[NotebookDirectory[]<>"graphdata_tau_multi3.m"];
+mergedData=Import[NotebookDirectory[]<>"graphdata_merged.m"];
+Export[NotebookDirectory[]<>"graphdata_merged_new.m",Join[mergedData,newData]]
+
+
+(* ::Subsection:: *)
+(*Plot empirical distribution of maximum degree*)
+
+
+maxDegrees=Map[{#[[1]],Max[#[[3]]]}&,gsraw];
+maxDegrees=GatherBy[maxDegrees,{#[[1,2]]&,#[[1,1]]&}];
+
+
+Histogram[maxDegrees[[1,-1,All,2]],PlotRange->{{0,1100},{0,200}},AxesLabel->{"d_max","frequency"}]
+Histogram[maxDegrees[[2,-1,All,2]],PlotRange->{{0,1100},{0,200}},AxesLabel->{"d_max","frequency"}]
+Histogram[maxDegrees[[3,-1,All,2]],PlotRange->{{0,1100},{0,200}},AxesLabel->{"d_max","frequency"}]
 
 
-Map[ListPlot[#[[2]],Joined->True,PlotRange->All]&,gsraw[[1;;3]]]
+(* ::Subsection:: *)
+(*Plot triangle count over "time" in Markov chain instances*)
 
 
-averages=Map[{#[[1]],Mean[#[[2,-1000;;-1]]]}&,gsraw];
-averagesGrouped=GatherBy[averages,#[[1]]&];
+numPlots=10;
+minCount=Min[Map[Min[#[[2]]]&,gsraw[[1;;numPlots]]]];
+maxCount=Max[Map[Max[#[[2]]]&,gsraw[[1;;numPlots]]]];
+maxTime=Max[Map[Length[#[[2]]]&,gsraw[[1;;numPlots]]]];
+skipPts=Round[maxTime/100]; (* Plotting every point is slow. Plot only once per `skipPts` timesteps *)
+coarseData=Map[#[[2,1;;-1;;skipPts]]&,gsraw[[1;;numPlots]]];
+labels=Map["{n,tau} = "<>ToString[#[[1]]]&,gsraw[[1;;numPlots]]];
+ListPlot[coarseData,Joined->True,PlotRange->{minCount,maxCount},DataRange->{0,maxTime},PlotLegends->labels]
+(* Map[ListPlot[#,Joined->True,PlotRange\[Rule]{minCount,maxCount},DataRange\[Rule]{0,maxTime}]&,coarseData] *)
+
+
+(* ::Subsection:: *)
+(*Plot average #triangles vs n*)
+
+
+averages=Map[{#[[1]],Mean[#[[2,1;;-1]]]}&,gsraw];
+averagesGrouped=GatherBy[averages,{#[[1,1]]&,#[[1,2]]&}];
+(* averagesGroupes[[ n index, tau index , run index, {ntau, tri, ds} ]] *)
+averagesGrouped=Transpose[averagesGrouped];
+(* averagesGrouped[[ tau index, n index, run index ]] *)
+nlabels=Map["n = "<>ToString[#]&,averagesGrouped[[1,All,1,1,1]]];
+taulabels=Map["tau = "<>ToString[#]&,averagesGrouped[[All,1,1,1,2]]];
 averagesErrorBars=Map[
-{{#[[1,1]],Mean[#[[All,2]]]},
+{{#[[1,1,1]],Mean[#[[All,2]]]},
 ErrorBar[StandardDeviation[#[[All,2]]]/Sqrt[Length[#]]]
-}&,averagesGrouped];
+}&,averagesGrouped,{2}];
+
+
+ErrorListPlot[averagesErrorBars,Joined->True,PlotMarkers->Automatic,AxesLabel->{"n","\[LeftAngleBracket]triangles\[RightAngleBracket]"},PlotLegends->taulabels]
+
+
+(* ::Subsection:: *)
+(*Plot #triangles distribution for specific (n,tau)*)
+
+
+histograms=Map[Histogram[#[[All,2]]]&,averagesGrouped,{2}];
+nlabels=Map["n = "<>ToString[#]&,averagesGrouped[[1,All,1,1,1]]];
+taulabels=Map["tau = "<>ToString[#]&,averagesGrouped[[All,1,1,1,2]]];
 
 
-ErrorListPlot[averagesErrorBars,Joined->True,PlotMarkers->Automatic,AxesLabel->{"n","\[LeftAngleBracket]triangles\[RightAngleBracket]"}]
+TableForm[histograms,TableHeadings->{taulabels,nlabels}]