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

sample8=d[SelectFirst[#n==8&]];
sample12=d[SelectFirst[#n==12&]];

c8[x_]=Cf[sample8][x];
c12[x_]=Cf[sample12][x];

LogLogPlot[{(1-3/π)x,c8[x],c12[x]},{x,0.0001,25},AxesLabel->{β,"c(β) / β"},PlotLegends->{"High-T expansion",n==8,n==12}]

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

fit8=ff[sample8,0.5,0.01];
fit12=ff[sample12,0.5,0.01];

Plot[Evaluate[{s8[x],fit8[x],s12[x],fit12[x]}],{x,0,25},AxesLabel->{β,(s[0]-s[β])}]

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

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