getPData[dir_]:=With[{data=Import[#],metadata=ToExpression/@Rest[StringSplit[FileBaseName[#],"_"]]},<|"n"->metadata[[1]],"β0"->metadata[[2]],"βf"->metadata[[3]],"nβ"->metadata[[4]],"As"->MapThread[MapThread[Around,{#1,#2}]&,{Rest[data],Sqrt[Rest[data]]}]/First[First[data]]|>]&/@FileNames[FileNameJoin[{dir,"probs_*.dat"}]]

getEData[dir_]:=With[{data=Import[#],metadata=ToExpression/@Rest[StringSplit[FileBaseName[#],"_"]]},<|"n"->metadata[[1]],"β0"->metadata[[2]],"βf"->metadata[[3]],"nβ"->metadata[[4]],"Es"->First@Rest[data]/First[First[data]]|>]&/@FileNames[FileNameJoin[{dir,"energies_*.dat"}]]

getρData[dir_]:=With[{data=Import[#],metadata=ToExpression/@Rest[StringSplit[FileBaseName[#],"_"]]},<|"n"->metadata[[1]],"β0"->metadata[[2]],"βf"->metadata[[3]],"nβ"->metadata[[4]],"ρs"->MapThread[MapThread[Around,{#1,#2}/Total[#1]]&,{data,Sqrt[data]}]|>]&/@FileNames[FileNameJoin[{dir,"rhos_*.dat"}]]

getData[dir_]:=Dataset[Merge[First]/@GatherBy[Join[getPData[dir],getEData[dir],getρData[dir]], {#n,#β0,#βf,#nβ}&]]

c[βi_,βj_][A_]:=2(((((βi)^2)+((βj)^2)))/(((βi-βj))^2))((InverseErfc[A])^2)

Ts[dat_]:=Range[dat["β0"],dat["βf"],(dat["βf"]-dat["β0"])/(dat["nβ"]-1)]

Cstmp[dat_]:=With[{ts=Ts[dat]},{{#,#+ts[[2]]-ts[[1]]}&/@ts,Last/@dat["As"]}]

Cs[dat_]:=With[{ts=Ts[dat]},MapThread[{Mean[{#1[[1]],#1[[2]]}],c[#1[[1]],#1[[2]]][#2]}&,{{#,#+ts[[2]]-ts[[1]]}&/@Most[ts],Last/@dat["As"]}]]

Z0[n_]:=(2^((n-1)(n/4+1)))Product[(π^(k/2))Gamma[k/2],{k,2,n}]

ΔSs[dat_]:=With[{ts=1/Ts[dat],cs=Cs[dat]},{Mean[1/#1],((#2(#1[[2]]-#1[[1]]))/(Mean[#1]))}&@@@Thread[{Rest@Rest[FoldList[{#1[[2]],#2}&,{0,0},ts]],Cs[dat][[All,2]]}]]

Ss[dat_]:=With[{ss=ΔSs[dat]},Thread[{ss[[All,1]],Accumulate[ss[[All,2]]]}]]

d=getData["~/doc/research/hadamard/code/"];

ListPlot[d[SortBy["n"],With[{L=((#n)^2)},{#1,#2/L}&@@@Cs[#]]&],Joined->True,PlotRange->All,PlotLegends->Normal@d[SortBy["n"],ToString[#n]&],AxesLabel->{β,c}]

ListPlot[d[SortBy["n"],With[{L=((#n)^2)},{#1,#2/L}&@@@Ss[#]]&],Joined->True,PlotRange->All,PlotLegends->Normal@d[SortBy["n"],ToString[#n]&],AxesLabel->{β,s[0]-s[β]}]

ListPlot[d[SortBy["n"],Thread[{Ts[#],#Es/((#n)^2)}]&],Joined->True,PlotLegends->Normal@d[SortBy["n"],ToString[#n]&],AxesLabel->{β,e}]

ListPlot[d[SortBy["n"],Thread[{Ts[#],10000/2#["ρs"][[All,5001]]}]&,FailureAction->None],Joined->True,PlotLegends->Normal@d[SortBy["n"],ToString[#n]&],AxesLabel->{β,(Subscript[ρ, 0])}]

ListPlot[d[SelectFirst[#n==4&],10000/2#ρs[[-1]]&],Joined->True]

tData={{4,Quantity[13,"Minutes"]+Quantity[32.95,"Seconds"]},{8,Quantity[42,"Minutes"]+Quantity[21.57,"Seconds"]},{16,Quantity[2,"Hours"]+Quantity[55,"Minutes"]+Quantity[5.77,"Seconds"]}};

tFit=LinearModelFit[Log2[{#1,#2/Quantity["Days"]}&@@@tData],x,x]

Show[ListPlot[Log2[{#1,#2/Quantity["Days"]}&@@@tData],Joined->True],Plot[tFit[x],{x,0,5},PlotStyle->ColorData[97][2]]]

2^tFit[5]