(* ::Package:: *)

gsraw=Import[NotebookDirectory[]<>"data/graphdata_canonical_mixingtime_histogram_merged.m"];
gdata=GatherBy[gsraw,{#[[1,2]]&}];
(* Data format: *)
(* gdata[[ tau index, n index, datatype index ]] *)
(* datatype index:
1: {n,tau}
2: { {time1, {samples1}}, {time2, {samples2}} , ... }  samples can also be a HistogramList kind of list
3: {uniform samples}
*)


(*Histogram[histogramdata[[All,2]],Automatic,"Probability",ChartLegends\[Rule]histogramdata[[All,1]]]*)

makeHistogram[run_,plotrange_,choices_]:=Module[{gridsize,labels,wd,uni},
uni=WeightedData[run[[3,All,1]],run[[3,All,2]]];(* Only in case of histogram-type file input *)
gridsize=Max[1,Mean[uni]/100];
labels="t = "<>ToString[NumberForm[N[#/run[[1,1]]],{3,1}]]<>" n"&/@run[[2,choices,1]];
wd=run[[2,choices,2]];
wd=Map[WeightedData[#[[All,1]],#[[All,2]]]&,wd];(* Only in case of histogram-type file input *)
Show[
SmoothHistogram[wd,gridsize,
PlotRange->plotrange,
PlotLegends->Placed[labels,Scaled[{0.7,0.7}]],
Epilog->Text["n = "<>ToString[run[[1,1]]]<>", \[Tau] = "<>ToString[run[[1,2]]],Scaled[{0.15,0.95}]],
ImageSize->300,
Frame->True,
FrameLabel->{"triangles","probability"}
],
SmoothHistogram[uni,gridsize,PlotRange->plotrange,PlotLegends->Placed[{"uniform"},Scaled[{0.7,0.7}]],PlotStyle->{Thick,Black}]]]

normalize[h_]:=Transpose[{h[[All,1]],h[[All,2]]/Total[h[[All,2]]]}];

makeHistogram2[run_,plotrange_,choices_]:=Module[{gridsize,labels,uni,probs},
uni=WeightedData[run[[3,All,1]],run[[3,All,2]]];(* Only in case of histogram-type file input *)
labels="t = "<>ToString[NumberForm[N[#/run[[1,1]]],{3,1}]]<>" n"&/@run[[2,choices,1]];
probs=run[[2,choices,2]];
probs=Map[normalize,probs];
Show[
ListPlot[probs,Joined->True,
PlotRange->plotrange,
PlotLegends->Placed[labels,Scaled[{0.7,0.7}]],
Epilog->Text["n = "<>ToString[run[[1,1]]]<>", \[Tau] = "<>ToString[run[[1,2]]],Scaled[{0.15,0.95}]],
ImageSize->300,
Frame->True,
FrameLabel->{"triangles","probability"}
],
ListPlot[normalize[run[[3]]],Joined->True,PlotRange->plotrange,PlotLegends->Placed[{"uniform"},Scaled[{0.7,0.7}]],PlotStyle->{Thick,Black}]]]

plot1=makeHistogram[gdata[[1,5]],{{3500,7800},{0,0.005}},{10,50,100}]
plot2=makeHistogram2[gdata[[3,5]],{{50,300},{0,0.04}},{10,20,30}]
plot3=makeHistogram2[gdata[[5,5]],{{0,70},{0,0.18}},{5,10,20}]
plot4=makeHistogram[gdata[[1,-1]],Automatic,{50,100,150,200}]


Export[NotebookDirectory[]<>"plots/triangle_distributions_over_time.pdf",plot2]


(* ::Section:: *)
(*Total variation distance*)


(* Get some sort of total variation distance between to sets of samples *)
getTVDistance1[samples1_,samples2_]:=Module[{max,probs1,probs2,binsize},
max=Max[Max[samples1],Max[samples2]];
binsize=Max[1,Floor[Mean[samples2]/500]];
binsize=1;
probs1=BinCounts[samples1,{0,max+binsize,binsize}];
probs1=probs1/Total[probs1];
probs2=BinCounts[samples2,{0,max+binsize,binsize}];
probs2=probs2/Total[probs2];
Total[Abs[probs1-probs2]]/2
]


(* Get some sort of total variation distance between to HistogramList-like objects *)
getTVDistance2[hist1_,hist2_,filterradius_]:=Module[{m,M,probs,probs1,probs2},
m=Min[hist1[[1,1]],hist2[[1,1]]];
M=Max[hist1[[-1,1]],hist2[[-1,1]]];
probs=ConstantArray[0,M-m+1];
probs1=hist1[[All,2]];
If[filterradius>1,probs1=GaussianFilter[probs1,10]]; (* test *) 
probs1=probs1/Total[probs1];
probs2=hist2[[All,2]];
If[filterradius>1,probs2=GaussianFilter[probs2,10]]; (* test *) 
probs2=probs2/Total[probs2];
probs[[1+hist1[[1,1]]-m;;1+hist1[[-1,1]]-m]]+=probs1;
probs[[1+hist2[[1,1]]-m;;1+hist2[[-1,1]]-m]]-=probs2;
Total[Abs[probs]]/2
]


(* ::Subsection:: *)
(*Plots*)


(*getRatios[run_]:=Module[{avg,sd,scalefactor},
avg=Mean[run[[3]]];
sd=StandardDeviation[run[[3]]];
scalefactor=1/(run[[1,1]]*Log[2,run[[1,1]]]^(1*(4-run[[1,2]])));
{"n,tau = "<>ToString[run[[1]]],
Map[{#[[1]]*scalefactor,Mean[#[[2]]]/avg}&,run[[2]]],
Map[{#[[1]]*scalefactor,(Mean[#[[2]]]-avg)/sd}&,run[[2]]],
}
]
ratios=Map[getRatios,gdata,{2}];
(*Map[ListPlot[#[[All,2]],Joined->True,PlotMarkers->Automatic,PlotRange->(1+{-0.15,+0.15}),PlotLegends->#[[All,1]],ImageSize->300]&,ratios]*)
(*Map[ListPlot[#[[All,3]],Joined->True,PlotMarkers->Automatic,PlotRange->(0+{-1,+1}),PlotLegends->#[[All,1]],ImageSize->300]&,ratios]*)
*)
getTVDs[run_]:=Module[{scaling},
scaling[step_,n_,tau_]:=step/(n*Log[n]^(4-tau));
{"n = "<>ToString[run[[1,1]]],"\[Tau] = "<>ToString[run[[1,2]]],
Map[{scaling[#[[1]],run[[1,1]],run[[1,2]]],getTVDistance2[#[[2]],run[[3]],30]}&,run[[2]]]}
]
TVDs=Map[getTVDs,gdata,{2}];


Map[Show[ListPlot[#[[All,3]],
Joined->True,(*PlotMarkers->Automatic,*)
PlotLegends->Placed[#[[All,1]],Center],ImageSize->500,
PlotRange->{0,1},Frame->True,FrameLabel->{"steps / (n (\!\(\*SubscriptBox[\(log\), \(2\)]\)n\!\(\*SuperscriptBox[\()\), \(4 - \[Tau]\)]\))","TV distance"},
PlotLabel->#[[1,2]]
],
Plot[{0.1,0.05},{x,0,20},PlotStyle->Directive[Black,Dashed]]
]&,TVDs]


(* ::Subsection:: *)
(*Mixing time from TV distance*)


getMixingTime[run_]:=Module[{i=1,timestamp=0,prevtimestamp=0},
While[i<=Length[run[[2]]],
    timestamp=run[[2,i,1]];
    If[getTVDistance2[run[[2,i,2]],run[[3]],20]<=0.1,
        Break[];
    ];
    prevtimestamp=timestamp;
    i++;
];
{run[[1,1]],run[[1,2]],prevtimestamp, timestamp} (* Range in which mixingtime lies *)
]
mixingtimes=Map[getMixingTime,gdata,{2}];


scaling[n_,tau_]:=n*Log[2,n]^(4-tau);
scaling[n_,tau_]:=n;
scaling[n_,tau_]:=n*Log[2,n]^1.5;

plotData=Map[{#[[1]],#[[4]]/scaling[#[[1]],#[[2]]]}&,mixingtimes,{2}]~Join~Map[{#[[1]],#[[3]]/scaling[#[[1]],#[[2]]]}&,mixingtimes,{2}];
fillings={1->{6},2->{7},3->{8},4->{9},5->{10}};
taulabels=Map["\[Tau] = "<>ToString[#[[1,2]]]&,mixingtimes];
ListPlot[plotData,
Joined->True,PlotMarkers->Automatic,ImageSize->500,
PlotRange->All,
PlotLegends->taulabels,
Filling->fillings,PlotStyle->{Automatic,Automatic,Automatic,Automatic,Automatic,None,None,None,None,None},
Frame->True,FrameLabel->{"n","ETMT / n"}]