getPData[dir_]:=With[{data=Import[#],metadata=ToExpression/@Rest[StringSplit[FileBaseName[#],"_"]]},<|"n"->metadata[[1]],"rs"->Partition[metadata[[2;;]],3],"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]],"rs"->Partition[metadata[[2;;]],3],"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]],"rs"->Partition[metadata[[2;;]],3],"ρs"->(#/Total[#]&)/@data(*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,#rs}&]]

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

Ts[dat_]:=Flatten[Rest[Range[#1,#2,(#2-#1)/#3]]&@@@dat["rs"]]

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]}&,{MapThread[#1+{0,#2}&,{Most[ts],Differences[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/"];

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

Export["/tmp/c.png",plotC];

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

Export["/tmp/s.png",plotS];

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

Export["/tmp/e.png",plotE];

plotR=ListLogPlot[Select[Length[#rs]==3&][d][SortBy["n"],Thread[{Ts[#],10000/2#["ρs"][[All,5001]]/((#n)^(1/2))}]&,FailureAction->None],Joined->True,PlotLegends->Normal@Select[Length[#rs]==1&][d][SortBy["n"],ToString[#n]&],AxesLabel->{β,(Subscript[ρ, 0])}]

Export["/tmp/ρ.png",plotR];

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

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

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"]},{32,Quantity[14,"Hours"]+Quantity[31,"Minutes"]+Quantity[15.87,"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[7]

UnitConvert[(Quantity[2,"Hours"]+Quantity[55,"Minutes"]+Quantity[5.77,"Seconds"])16 2,"Hours"]