BeginPackage["Schofield`"]

$Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0}

β[D_:2] := Piecewise[
  {
    {1/8, D == 2},
    {0.326419, D == 3},
    {1/2, D == 4},
    {β, True}
  }
]

δ[D_:2] := Piecewise[
  {
    {15, D == 2},
    {4.78984, D == 3},
    {3, D == 4},
    {δ, True}
  }
]

α[D_:2] := Piecewise[
  {
    {0, D == 2},
    {0.11008, D == 3},
    {0, D == 4},
    {α, True}
  }
]

Δ[D_:2] := β[D] δ[D]

os := 2^(1/12) Exp[-1/8] Glaisher^(3/2)

t[θ_] := ((θ/1)^2 - 1)
h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i] θ^(2*i+1), {i, 0, n}]
η[g_][θ_] := t[θ] / (g[θ] / I)^(1 / Δ[2])

fLow[B_, θc_][θ_] := (θc Exp[1/(B θc)] ExpIntegralEi[-1/(B θc)] + (θ - θc) Exp[-1/(B (θ - θc))] ExpIntegralEi[1/(B (θ - θc))]) / π

RFLow[B_, θc_][θ_] := fLow[B, θc][θ] + fLow[B, θc][-θ]
RFHigh[ξ0_][ξ_] := (-I ξ+ ξ0)^(5/6) + (I ξ+ ξ0)^(5/6)  - 2 ξ0^(5/6)
RFReg[n_][θ_] := AL fLow[B, θc][-θ] + AH RFHigh[θ0][θ] + Sum[A[i] θ^(2 i), {i, 1, n}]
RF[n_][θ_] := RFReg[n][θ] + AL fLow[B, θc][θ]
dfLow[AL_, B_][m_] := Piecewise[{{AL m! Gamma[m - 1] B^(m - 1) / π, m > 1}, {AL θc Exp[1/(B θc)] ExpIntegralEi[-1/(B θc)] / π, m == 0}, {0, True}}]
dRFc[n_][m_] := dfLow[AL, B][m] + D[RFReg[n][θ], {θ, m}] /. θ -> θc

RFC[n_][θ_] := RF[n][θ] + AL I Sign[Im[θ]] ((θ-θc)Exp[-1/(B(θ-θc))]-(-θ-θc)Exp[-1/(B(-θ-θc))])

invDerivativeList[n_][f_][x_] := Module[
    {xp, dfs, fp, Pns},
  dfs = Rest[NestList[D[#, xp] &, f[xp], n]] /. xp -> x;
  Pns = FoldList[{Pm, m} \[Function]
        fp'[xp] D[Pm, xp] - (2 m - 1) fp''[xp] Pm, 1, Range[n - 1]] /. 
    Derivative[m_][fp][xp] :> dfs[[m]];
  MapIndexed[{Pn, i} \[Function] Pn/dfs[[1]]^(2 i[[1]] - 1), Pns]
  ]

dGdξLowList[n_, h_, params_:{}][m_] := Module[
  { ds, dF, df },
  ds = invDerivativeList[m+1][Function[θ, h[θ] / t[θ]^Δ[2] //. params]][θc //. params];
  dF = NestList[Function[f, D[f, θ]], RFReg[n][θ], m] + Table[dfLow[AL, B][k], {k, 0, m}] /. θ -> θc //. params;
  df = NestList[D[#, \[Theta]] &, 
      fp[\[Theta]]/t[\[Theta]]^2 - 1/(8 \[Pi]) Log[t[\[Theta]]^2], 
      m] /. 
   Map[Derivative[#][fp][\[Theta]] -> dF[[# + 1]] &, 
    Range[0, m]] /. θ -> θc //. params;
  Table[Sum[df[[k+1]] BellY[j, k, ds[[;; j - k + 1]]], {k, 0, j}]/(j!), {j, 0, m}]
]

dGdξList[n_, h_, params_:{}][m_, θp_] := Module[
  { ds, dF, df },
  ds = invDerivativeList[m+1][Function[θ, h[θ] / RealAbs[t[θ]]^Δ[2] //. params]][θp];
  dF = NestList[Function[f, D[f, θ]], RF[n][θ] //. params, m] /. θ -> θp //. params;
  df = NestList[D[#, \[Theta]] &, 
      fp[\[Theta]]/t[\[Theta]]^2 - 1/(8 \[Pi]) Log[t[\[Theta]]^2], 
      m] /. 
   Map[Derivative[#][fp][\[Theta]] -> dF[[# + 1]] &, 
    Range[0, m]] /. θ -> θp //. params;
  Table[Sum[df[[k+1]] BellY[j, k, ds[[;; j - k + 1]]], {k, 0, j}]/(j!), {j, 0, m}]
]

dΦdηList[n_, h_, params_:{}][m_, θp_] := Module[
  { ds, dF, df },
  ds = invDerivativeList[m+1][Function[θ, t[θ] / h[θ]^(1 / Δ[2]) //. params]][θp];
  dF = NestList[Function[f, D[f, θ]], RF[n][θ] //. params, m] /. θ -> θp;
  df = NestList[D[#, \[Theta]] &, 
      hh[θ]^(-2 / Δ[]) (fp[θ] - t[θ]^2 Log[hh[θ]^2] / (8 π Δ[])), 
      m] /. 
   Map[Derivative[#][fp][\[Theta]] -> dF[[# + 1]] &, 
    Range[0, m]] /. hh -> h /. θ -> θp //. params;
  Table[Sum[df[[k+1]] BellY[j, k, ds[[;; j - k + 1]]], {k, 0, j}]/(j!), {j, 0, m}]
]

ruleB[g_] := B -> (2 * os / π) * (- g'[θc] / t[θc]^Δ[2])
ξYL[g_] := g[I θ0]/(-t[I θ0])^Δ[2]/I
AYL[g_] := AH / ((g[I θ0]/ I)^(2 / Δ[2]) * (-η[g]'[I θ0] / (I))^(5/6))
ruleθ0[g_] := ξYL[g] - Around[0.18930, 0.00005]
ruleAL[g_] := AL -> Exp[Δ[2] t[θc]^(Δ[2] - 1) t'[θc] / (2 os / π g'[θc]) - t[θc]^Δ[2] g''[θc] / (4 os / π g'[θc]^2)] t[θc]^(1/8) os / (2 π) * g'[θc]
ruleAH[g_] := AYL[g] + Around[1.37, 0.02]

eqLowRHSReg[n_][m_] := dRFc[n][m]

eqLowLHS[h_][m_] :=D[
  t[θ]^2 (Gl[h[θ] t[θ]^-Δ[2]] + Log[t[θ]^2]/(8 π)),
  {θ, m} ] /. θ -> θc

eqLow[n_, h_][m_] := (eqLowRHSReg[n][m] - eqLowLHS[h][m]) / m!

eqHighRHS[F_][m_] := D[F[θ], {θ, m} ] /. θ -> 0

eqHighLHS[h_][m_] := D[(-t[θ])^2 (Gh[h[θ] (-t[θ])^-Δ[2]] + Log[(-t[θ])^2]/(8 π)), {θ, m} ] /. θ -> 0

eqHigh[n_, h_][m_] := (eqHighRHS[RF[n]][m] - eqHighLHS[h][m]) / m!

eqMid[F_, h_][m_] := D[
  F[θ] - t[θ]^2 Log[h[θ]^2]/(8 Δ[2]π) - h[θ]^((2-α[2])/Δ[2]) Φ[η]
    /. η -> t[θ] / h[θ]^(1 / Δ[2]),
  {θ, m} ] / m! /. θ -> 1

Φs = {
  -Gamma[1/3]Gamma[1/5]Gamma[7/15]/(2 π Gamma[2/3]Gamma[4/5]Gamma[8/15])(4 π^2 Gamma[13/16]^2 Gamma[3/4]/(Gamma[3/16]^2 Gamma[1/4]))^(8/15),
  -0.31881012489061,
(*
  Sin[π s(pp l - p k)/p]/Sin[π s/p(pp-p)](M π(ξ+1)Gamma[(2 ξ+2)/(3 ξ+6)]/(2^(2/3)Sqrt[3]Gamma[1/3]Gamma[ξ/(3 ξ+6)]))^(2 Δlk) Q[ξ, (ξ+1)l-ξ k] //. {
    l -> 1, k -> 3, p -> 3, pp -> 4, ξ -> p/(pp-p), Δlk -> ((pp l-p k)^2-(pp-p)^2)/(4 p pp),
    M -> 2^((ξ+5)/(3 ξ+6)) Sqrt[3] Gamma[1/3]Gamma[ξ/(3 ξ+6)]/π/Gamma[(2 ξ+2)/(3 ξ+6)](
      π^2 Gamma[(3 ξ+4)/(4 ξ+4)]^2 Gamma[1/2+1/(ξ+1)]/Gamma[ξ/(4 ξ+4)]^2/Gamma[1/2-1/(ξ+1)]
    )^((ξ+1)/(3 ξ + 6))
  } /. Q -> Function[{ξ, η}, Exp[
      NIntegrate[
      (Sinh[(ξ+2)t]Sinh[t(η-1)]Sinh[t(η+1)]/Sinh[3(ξ+2)t]/Sinh[2(ξ+1)t]/Sinh[ξ t](Cosh[3(ξ+2)t]+Cosh[(ξ+4)t]-Cosh[(3 ξ+4)t]+Cosh[ξ t]+1)-(η^2-1)/(2 ξ(ξ+1))Exp[-4 t]) /t
      ,{t, 0, Infinity}, WorkingPrecision->50]
    ]], 
*)
  Around[0.110886196683, 2.0 10^-12],
  Around[0.01642689465,  2.0 10^-11],
  Around[-2.639978 10^-4, 1.0 10^-10],
  Around[-5.140526 10^-4, 1.0 10^-10],
  Around[2.08865 10^-4,   1.0 10^-9],
  Around[-4.4819 10^-5,   1.0 10^-9],
  Around[3.16 10^-7,      1.0 10^-9],
  Around[4.31 10^-6,      0.01 10^-6],
  Around[-1.99 10^-6,     0.01 10^-6]
}

Gls = {
  0,
  -os,
  −1.000960328725262189480934955172097320572505951770117 Sqrt[2]/((2 )^(-7/8) (2^(3/16)/os)^2)/2/(12 \[Pi]),
  Around[ 0.038863932, 3.0 10^(-9)],
  Around[−0.068362119, 2.0 10^(-9)],
  Around[ 0.18388370,  1.0 10^(-8)],
  Around[-0.6591714,   1.0 10^(-7)],
  Around[ 2.937665,    3.0 10^(-6)],
  Around[-15.61,       1.0 10^(-2)],
  Around[ 96.76,       1.0 10^(-2)],
  Around[-6.79 10^2,   1.0],
  Around[ 5.34 10^3,    10.],
  Around[-4.66 10^4,   0.01 10^4],
  Around[ 4.46 10^5, 0.01 10^5],
  Around[-4.66 10^6, 0.01 10^6]
}

Ghs = {
  0,
  0,
  -1.000815260440212647119476363047210236937534925597789 Sqrt[2]/((2 )^(-7/8) (2^(3/16)/os)^2)/2,
  0,
  Around[  8.333711750, 5.0 10^(-9)],
  0,
  Around[-95.16896,     1.0 10^(-5)],
  0,
  Around[1457.62,       3.0 10^(-2)],
  0,
  Around[-2.5891 10^4,  2.0],
  0,
  Around[5.02 10^5, 0.01 10^5],
  0,
  Around[-1.04 10^7, 0.01 10^7]
}

dRule[sym_][f_, i_] := Derivative[i[[1]] - 1][sym][0] -> f (i[[1]] - 1)!

ΦRules = MapIndexed[dRule[Φ], Φs];
GlRules = MapIndexed[dRule[Gl], Gls];
GhRules = MapIndexed[dRule[Gh], Ghs];

rules[g_] := Join[ΦRules, GlRules, GhRules, {ruleAL[g], ruleB[g]}]

eqAround[n_, g_][m_, p_, q_] := Flatten[Join[{ruleAH[g], ruleθ0[g]}, eqLow[n, g][#] & /@ Range[0, m],eqMid[RF[n], g][#] & /@ Range[0, p], eqHigh[n, g] /@ Range[2, q, 2]]] //. rules[g]

formResiduals[data_, functions_, δ_:10^(-15)] := If[Head[#1]===Around,
  (#2-#1["Value"]) / Max[#1["Uncertainty"], δ],
  (#2-#1) / δ] & @@@ Thread[{data, functions}]

resLow[n_, g_, δ_][m_] := formResiduals[Gls[[;;m+1]], dGdξLowList[n, g][m], δ]

resHigh[n_, g_, δ_][m_] := Rest[formResiduals[Ghs[[;;m+1]], dGdξList[n, g][m, 0], δ][[;;;;2]]]

resMid[n_, g_, δ_][m_] := formResiduals[Φs[[;;m+1]], dΦdηList[n, g][m, 1], δ]

res[F_, g_, δ_][m_] := Join[resLow[F, g, δ][m], resHigh[F, g, δ][m]] //. rules[g]
chiSquared[F_, g_, δ_][m_] := Total[res[F, g, δ][m]^2]
resAll[F_, g_, δ_][m_] := Join[resLow[F, g, δ][m], resHigh[F, g, δ][m], (*resMid[F, g, δ][m],*) {ruleθ0[g] / 0.00005, ruleAH[g] / 0.02} /. Around[x_, _] :> x] //. rules[g]
resHighAll[F_, g_, δ_][m_] := Join[resLow[F, g, δ][m][[{1}]], resHigh[F, g, δ][m], (*resMid[F, g, δ][m],*) {ruleθ0[g] / 0.00005, ruleAH[g] / 0.02} /. Around[x_, _] :> x] //. rules[g]

newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[
  eqs,
  Join[
    {#1, #2 + γ * RandomVariate[NormalDistribution[]]} & @@@ (oldSol /. Rule -> List),
    Thread[{newVars, δ * RandomVariate[NormalDistribution[], Length[newVars]]}]
  ],
  MaxIterations -> 10000,
  opts
]

levenbergMarquardtStep[M_, λ_, g_] := Quiet[
  - LinearSolve[M + λ DiagonalMatrix[Max[10^-5, #] & /@Diagonal[M]], g]
  , LinearSolve::luc
]

levenbergMarquardtAcceleration[M_, λ_, v_, J_, rold_, rf_, h_, β_] := Quiet[
  - Re[LinearSolve[M + λ DiagonalMatrix[Max[10^-5, #] & /@Diagonal[M]], Transpose[J].(2/h((rf[β[[All,2]] + h v] - rold) / h - J . v))] / 2]
  , LinearSolve::luc
]

levenbergMarquardt[r_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] := 
 Module[
  {
    n = Length[β0],
    rf, Jf
   },
  PrintTemporary["Compiling the Jacobian..."];
  Quiet[
    rf = Compile[{{x, _Real, 1}}, 
      Evaluate[r /. Thread[Rule[First /@ β0, Part[x, #] & /@ Range[n]]]]
    ];
    Jf = Compile[{{x, _Real, 1}}, 
      Evaluate[D[r, {First /@ β0}] /. Thread[Rule[First /@ β0, Part[x, #] & /@ Range[n]]]]
    ];
  , Part::partd];
  levenbergMarquardtHelper[1, rf, Jf, β0, λ0, ν, ε]
 ]


levenbergMarquardtHelper[Δ_, rf_, Jf_, β0_, λ0_ : 1, ν_ : 1.5, ε_ : 10^-15] := 
 Module[
  {
    n = Length[β0],
    β = β0,
    λ = λ0,
    newβ = β0,
    x, oldJ, oldr, newr, M, g, δ, oldC, newC, vlast,
    h = 0.1, α = 0.1, a, v, μ = 5
   },
  oldr = rf[Append[β[[All, 2]], Δ]];
  oldJ = Jf[Append[β[[All, 2]], Δ]];
  oldC = Re[Total[oldr^2]];
  newC = oldC;
  g = Re[Transpose[oldJ] . oldr];
  M = Re[Transpose[oldJ] . oldJ];
  v = levenbergMarquardtStep[M, λ, g];
  a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β, Δ];
  δ = v + a;
  vlast = v;
  PrintTemporary["Current cost value: ", Dynamic[oldC]];
  While[Norm[δ] > ε,
    newβ[[All, 2]] += δ;
    newr = rf[Append[newβ[[All, 2]], Δ]];
    newC = Re[Total[newr^2]];
    While[((1 - v . vlast / Norm[v] / Norm[vlast])newC > oldC) || (2 Norm[a] / Norm[v] > α),
      λ *= ν;
      v = levenbergMarquardtStep[M, λ, g];
      a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β, Δ];
      δ = v + a;
      newβ = β;
      newβ[[All, 2]] += δ;
      newr = rf[Append[newβ[[All, 2]], Δ]];
      newC = Re[Total[newr^2]];
    ];
    oldC = newC;
    oldr = newr;
    β = newβ;
    oldJ = Jf[Append[β[[All, 2]], Δ]];
    g = Re[Transpose[oldJ] . oldr];
    M = Re[Transpose[oldJ] . oldJ];
    vlast = v;
    λ /= μ;
    v = levenbergMarquardtStep[M, λ, g];
    a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β, Δ];
    δ = v + a;
  ];
  {newC, β}
 ]

annealFit[r_, β0_, γ_ : 0, λ0_ : 1, ν_ : 2, μ_ : 3, ε_: 10^-15, h_ : 0.1, α_ : 0.1] :=
Module[ {
    rf, Jf,
    n = Length[β0],
    β = β0,
    λ = λ0,
    newβ = β0,
    x, oldJ, oldr, newr, M, g, δβ, oldC, newC, vlast,
    a, v
    },
  Quiet[
    rf = Compile[{{x, _Real, 1}}, 
      Evaluate[r /. Thread[Rule[First /@ β0, Part[x, #] & /@ Range[n]]]]
    ];
    Jf = Compile[{{x, _Real, 1}}, 
      Evaluate[D[r, {First /@ β0}] /. Thread[Rule[First /@ β0, Part[x, #] & /@ Range[n]]]]
    ];
  , Part::partd];
  If[γ > 0,
  β[[All, 2]] += RandomVariate[NormalDistribution[0, γ], n]];
  newβ = β;
  λ = λ0;
  oldr = rf[β[[All, 2]]];
  oldJ = Jf[β[[All, 2]]];
  oldC = Re[Total[oldr^2]];
  newC = oldC;
  g = Re[Transpose[oldJ] . oldr];
  M = Re[Transpose[oldJ] . oldJ];
  v = levenbergMarquardtStep[M, λ, g];
  a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β];
  δβ = v + a;
  vlast = v;
  While[Norm[δβ] > ε,
    newβ[[All, 2]] += δβ;
    newr = rf[newβ[[All, 2]]];
    newC = Re[Total[newr^2]];
    While[((1 - v . vlast / Norm[v] / Norm[vlast])newC > oldC) || (2 Norm[a] / Norm[v] > α),
      λ *= ν;
      v = levenbergMarquardtStep[M, λ, g];
      a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β];
      δβ = v + a;
      newβ = β;
      newβ[[All, 2]] += δβ;
      newr = rf[newβ[[All, 2]]];
      newC = Re[Total[newr^2]];
    ];
    oldC = newC;
    oldr = newr;
    β = newβ;
    oldJ = Jf[β[[All, 2]]];
    g = Re[Transpose[oldJ] . oldr];
    M = Re[Transpose[oldJ] . oldJ];
    vlast = v;
    λ /= μ;
    v = levenbergMarquardtStep[M, λ, g];
    a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β];
    δβ = v + a;
  ];
  {oldC, β}
]

fitCov[res_, s_] := 
 With[{J = D[res, {s[[All, 1]]}]}, 
  Symmetrize[Inverse[SetPrecision[Transpose[J] . J /. s, 60]]]]

EndPackage[]