getData[dir_]:=Dataset[With[{data=Import[#],metadata=ToExpression/@Rest[StringSplit[FileBaseName[#],"_"]]},<|"n"->metadata[[1]],"β0"->metadata[[2]],"β1"->metadata[[3]],"N"->metadata[[4]],"βs"->data[[1]],"f"->data[[2]],"ps"->data[[3]],"Es"->data[[4]],"Cs"->data[[5]],"ρs"->(#/Total[#]&)/@data[[6;;]]|>]&/@FileNames[FileNameJoin[{dir,"hmm_*.dat"}]]]

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

Ts[dat_]:=dat["βs"]

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

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]}],dat["ps"]}]]

logZ0[n_]:=(n-1)(n/4+1)Log[2]+Sum[(k/2)Log[π]-Log[Gamma[k/2]],{k,2,n}]

logZ[n_][β_]:=lognum[n]+β (n^2)+n(n-1)/4Log[((2π)/(β n))]

ΔSs[dat_]:=With[{ts=Ts[dat],cs=Cs[dat]},{Mean[#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]]]}]]

ss8=FullSimplify[(logZ[8][β]-logZ0[8])/(8^2)-β,Assumptions->{β>0}]

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

c8=Evaluate[d[3,Interpolation[Normal@Prepend[With[{L=((#n)^2)},{#1,#2/L/#1}&@@@Thread[{#βs,#Cs ((#βs)^2)}]],{0,0}]]&]]

d[SortBy["n"],#Cs&]//Normal

ListLogLogPlot[{d[2,Cs],d[2,Thread[{#βs,#Cs((#βs)^2)}]&]},Joined->True]

c16=Evaluate[d[SelectFirst[#n==16&&#N==256&],Interpolation[Normal@Prepend[With[{L=((#n)^2)},{#1,#2/L/#1}&@@@Cs[#]],{0,0}]]&]]

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

Show[LogLogPlot[{c8[x],(1-3/π)x},{x,0.0001,30},PlotRange->All],ListLogLogPlot[d[3,With[{L=((#n)^2)},{#1,#2/L/#1}&@@@Thread[{#βs,#Cs ((#βs)^2)}]]&],PlotRange->All]]

s8=Table[{X,-NIntegrate[c8[x],{x,0,X},WorkingPrecision->40]},{X,0,19.0,0.25}];

s8fit=NonlinearModelFit[Select[s8,#[[1]]>15&],ss8,{lognum[8]},β,MaxIterations->1000]

Show[ListPlot[s8],Plot[s8fit[x],{x,0,20}]]

Exp[lognum[8]]/.s8fit["BestFitParameters"]

fits=Normal@d[SortBy["n"],With[{L=((#n)^2)},NonlinearModelFit[Select[{#1,#2/L}&@@@Ss[#],#[[1]]>15&],0.48029-(1-1/#n) Log[2π E/(β #n)]/4-logN/L,{a,logN},β]]&]

Exp[logN]/.s16fit["BestFitParameters"]

4954521600./2

88526812916367202104587059200000./2

3776127947893930552689423154306445475840000000.

ListPlot[d[SortBy["n"],Thread[{#βs,Range[#N]}]&],Joined->True]

ListPlot[d[SortBy["n"],#f&],Joined->True]

ListPlot[d[SortBy["n"],Thread[{Most@#βs,#ps}]&],Joined->True,PlotRange->All]

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

plotS=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[β]-s[0]}]

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

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

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

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

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

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

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"]