diff --git a/triangle_gcm_initial_analysis.m b/triangle_gcm_initial_analysis.m
new file mode 100644
index 0000000000000000000000000000000000000000..3144fde4c5d8d35253d9c60ae270b2e62bcb80b1
--- /dev/null
+++ b/triangle_gcm_initial_analysis.m
@@ -0,0 +1,83 @@
+(* ::Package:: *)
+
+Needs["ErrorBarPlots`"]
+
+
+(* ::Section:: *)
+(*Data import*)
+
+
+gsraw=Import[NotebookDirectory[]<>"data/graphdata_initialtris.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. *) *)
+
+
+gdata=GatherBy[gsraw,{#[[1,2]]&,#[[1,1]]&}];
+(* Data format: *)
+(* gdata[[ tau index, n index, run index , datatype index ]] *)
+(* datatype index:
+1: {n,tau}
+2: avgtriangles when mixed
+3: {GCM1 starting triangle counts summed, GCM1 number of successes}
+4: {GCM2 starting triangle counts summed, GCM2 number of successes}
+*)
+nlabels=Map["n = "<>ToString[#]&,gdata[[1,All,1,1,1]]];
+taulabels=Map["tau = "<>ToString[#]&,gdata[[All,1,1,1,2]]];
+
+
+(* ::Section:: *)
+(*Greedy configuration model*)
+
+
+(* ::Subsection:: *)
+(*Distribution of initial #triangles for GCM1,GCM2,EG compared to average #triangles.*)
+
+
+ (* Consider all runs *)
+getIt[x_,avg_]:=If[x[[2]]>=5, x[[1]]/x[[2]]/avg,-0.5]
+getAverage[run_]:=Module[{avg},
+    avg=20*run[[2]];
+    If[avg>0,
+    { getIt[run[[3]],avg], getIt[run[[4]],avg] }
+    , {3,3}]
+]
+getTotalStats[runs_]:=Transpose[Map[getAverage,runs]];
+totalStats=Map[getTotalStats,gdata,{2}];
+
+
+(* Yellow: GCM1 (take new highest everytime *)
+(* Blue: GCM2 (finish highest first, more similar to EG) *)
+histogramsTotal=Map[Histogram[#,{0.05},PlotRange->{{-0.5,2},Automatic},ImageSize->300,AxesOrigin->{0,0}]&,totalStats,{2}];
+
+
+TableForm[histogramsTotal,TableHeadings->{taulabels,nlabels}]
+
+
+(* ::Subsubsection:: *)
+(*Exporting plots*)
+
+
+tauIndices={2,5,8};
+nIndices={2};
+makeHistogram[datasets_,n_,tau_]:=Histogram[datasets,{0.05},PlotRange->{{-0.5,1.5},Automatic},ImageSize->300,AxesOrigin->{0,0}];
+histogramsTotal=Map[makeHistogram[#,1000,tau]&,totalStats[[tauIndices,nIndices]],{2}];
+TableForm[histogramsTotal,TableHeadings->{taulabels[[tauIndices]],nlabels[[nIndices]]}]
+
+
+(* ::Subsection:: *)
+(*GCM1 vs GCM2 success rates*)
+
+
+(* gdata[[ tau index, n index, run index , {ntau, #tris, ds, greedyTriangles} ]] *)
+successrates=Map[{#[[3,2]],#[[4,2]]}&,gdata,{3}];
+successrates=Map[Transpose,successrates,{2}];
+successratesDelta=Map[#[[3,2]]-#[[4,2]]&,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}]*)
+
+
+