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[[1]])^2)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_][β_]:=logN[n]+n(n-1)/4Log[((2π E)/(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]]}]]

Cf[dat_][x_]:=Interpolation[Normal@Prepend[Thread[{Normal@dat["βs"],Normal@dat["Cs"]/Normal@dat["βs"]/((dat["n"])^2)}],{0,0}]][x]

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

Sf[dat_][x_]:=Integrate[Cf[dat][x],x]

ff[dat_,Δx_,δx_]:=Module[{x},Module[{s=Sf[dat][x],n=dat["n"],β,xmax=dat["β1"]},
NonlinearModelFit[Table[{x,s},{x,xmax-Δx,xmax,δx}],(logZ0[n]-logZ[n][β])/(n^2),logN[n],β]
]]

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

d[All,{#n,#N,#β1}&]

sample8=d[SelectFirst[#n==8&&#N==128&]];
sample12=d[Select[#n==12&&#N==128&&#β1==50&]][[2]];
sample16=d[SelectFirst[#n==16&&#N==64&]];
sample20=d[SelectFirst[#n==20&&#N==64&]];

c8[x_]=Cf[sample8][x];
c12[x_]=Cf[sample12][x];
c16[x_]=Cf[sample16][x];
c20[x_]=Cf[sample20][x];

Plot[{x c8[x],x c12[x],x c16[x]},{x,0.0001,10},AxesLabel->{β,"c(β)"},PlotRange->All,PlotLegends->{n==8,n==12,n==16}]

Show[LogLogPlot[{(1-3/π)x,c8[x],c12[x],c16[x]},{x,0.0001,25},AxesLabel->{β,"c(β) / β"},PlotLegends->{"High-T expansion",n==8,n==12,n==16},PlotStyle->Prepend[ColorData[97]/@Range[3],Black]]]

s8[x_]=Sf[sample8][x];
s12[x_]=Sf[sample12][x];
s16[x_]=Sf[sample16][x];

fit8=ff[sample8,1,0.01];
fit12=ff[sample12,1,0.01];
fit16=ff[sample16,1,0.01];

LogLinearPlot[Evaluate[{s8[x],fit8[x],s12[x],fit12[x],s16[x],fit16[x]}],{x,1,25},AxesLabel->{β,(s[0]-s[β])},PlotStyle->ColorData[97]/@Riffle[Range[3],Range[3]],PlotLegends->{n==8,None,n==12,None,n==16,None}]

ns=Import["https://oeis.org/A206711/b206711.txt","Table"][[2;;,2]]/2;

(Exp[logN[8]]-ns[[2]])/ns[[2]]/.fit8["BestFitParameters"]

(Exp[logN[12]]-ns[[3]])/ns[[3]]/.fit12["BestFitParameters"]

(Exp[logN[16]]-ns[[4]])/ns[[4]]/.fit16["BestFitParameters"]

ρ8=sample8[Thread[{#βs,5000#ρs[[All,5001]]/Sqrt@#n}]&];
ρ12=sample12[Thread[{#βs,5000#ρs[[All,5001]]/Sqrt@#n}]&];
ρ16=sample16[Thread[{#βs,5000#ρs[[All,5001]]/Sqrt@#n}]&];
ρ20=sample20[Thread[{#βs,5000#ρs[[All,5001]]/Sqrt@#n}]&];

ListLogPlot[{ρ8,ρ12,ρ16},PlotRange->{{0,10},{0.0001,0.5}},Joined->True,PlotLegends->{n==8,n==12,n==16},AxesLabel->{β,(Subscript[ρ, 0])[β]}]

E8=sample8[Thread[{#βs,#Es/((#n)^2)}]&];
E12=sample12[Thread[{#βs,#Es/((#n)^2)}]&];
E16=sample16[Thread[{#βs,#Es/((#n)^2)}]&];
E20=sample20[Thread[{#βs,#Es/((#n)^2)}]&];

ListPlot[{E8,E12,E16},PlotRange->{{0,10},All},Joined->True,PlotLegends->{n==8,n==12,n==16},AxesLabel->{β,e[β]}]

ListPlot[d[All,Thread[{#βs,Range[#N]/Sqrt[#N]}]&],Joined->True]

ListLogPlot[d[SelectFirst[#n==8&],10000/2#ρs[[-4]]&],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"]

UnitConvert[Quantity[415330.76,"Seconds"],"Hours"]

UnitConvert[Quantity[528330.62,"Seconds"],"Hours"]