summaryrefslogtreecommitdiff
path: root/schofield.wl
diff options
context:
space:
mode:
Diffstat (limited to 'schofield.wl')
-rw-r--r--schofield.wl109
1 files changed, 77 insertions, 32 deletions
diff --git a/schofield.wl b/schofield.wl
index 3f45f5d..98f9721 100644
--- a/schofield.wl
+++ b/schofield.wl
@@ -32,7 +32,7 @@ $Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0}
Δ[D_:2] := β[D] δ[D]
-OverBar[s] := 2^(1/12) Exp[-1/8] Glaisher^(3/2)
+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}]
@@ -95,11 +95,11 @@ dΦdηList[n_, h_][m_, θp_] := Module[
Table[Sum[df[[k+1]] BellY[j, k, ds[[;; j - k + 1]]], {k, 0, j}]/(j!), {j, 0, m}]
]
-ruleB[g_] := B -> (2 * OverBar[s] / π) * (- g'[θc] / t[θc]^Δ[2])
+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] / (2 θ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 OverBar[s] / π g'[θc]) - t[θc]^Δ[2] g''[θc] / (4 OverBar[s] / π g'[θc]^2)] t[θc]^(1/8) OverBar[s] / (2 π) * g'[θc]
+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]
@@ -137,8 +137,8 @@ eqMid[F_, h_][m_] := D[
Gls = {
0,
- -OverBar[s],
- −1.000960328725262189480934955172097320572505951770117 Sqrt[2]/((2 )^(-7/8) (2^(3/16)/OverBar[s])^2)/2/(12 \[Pi]),
+ -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.18388371, 1.0 10^(-8)],
@@ -156,7 +156,7 @@ Gls = {
Ghs = {
0,
0,
- -1.000815260440212647119476363047210236937534925597789 Sqrt[2]/((2 )^(-7/8) (2^(3/16)/OverBar[s])^2)/2,
+ -1.000815260440212647119476363047210236937534925597789 Sqrt[2]/((2 )^(-7/8) (2^(3/16)/os)^2)/2,
0,
Around[ 8.333711750, 5.0 10^(-9)],
0,
@@ -189,7 +189,7 @@ 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]]]
-res[F_, g_, δ_][m_] := Join[resLow[F, g, δ][m], resHigh[F, g, δ][m]]
+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]
newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[
@@ -203,13 +203,12 @@ newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[
]
levenbergMarquardtStep[M_, λ_, g_] := Quiet[
- - LinearSolve[M + λ DiagonalMatrix[Diagonal[M]], g]
+ - LinearSolve[M + λ DiagonalMatrix[Max[10^-5, #] & /@Diagonal[M]], g]
, LinearSolve::luc
]
levenbergMarquardtAcceleration[M_, λ_, v_, J_, rold_, rf_, h_, β_] := Quiet[
- δβ = β; δβ[[All, 2]] += h v;
- - Re[LinearSolve[M + λ DiagonalMatrix[Diagonal[M]], Transpose[J].(2/h(((rf /. δβ) - rold) / h - J . v))] / 2]
+ - 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
]
@@ -232,7 +231,7 @@ levenbergMarquardt[r_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] :=
]
-levenbergMarquardtHelper[r_, J_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] :=
+levenbergMarquardtHelper[Δ_, rf_, Jf_, β0_, λ0_ : 1, ν_ : 1.5, ε_ : 10^-15] :=
Module[
{
n = Length[β0],
@@ -240,61 +239,107 @@ levenbergMarquardtHelper[r_, J_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] :=
λ = λ0,
newβ = β0,
x, oldJ, oldr, newr, M, g, δ, oldC, newC, vlast,
- h = 0.1, fvv, α = 0.1, a, v, U, nSteps = 0, reject, μ =3
+ h = 0.1, α = 0.1, a, v, μ = 5
},
- PrintTemporary["Beginning the algorithm."];
- oldr = r /. β;
- oldJ = J /. β;
- U = First[SingularValueDecomposition[oldJ]];
+ 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, r, h, β];
+ a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β, Δ];
δ = v + a;
vlast = v;
PrintTemporary["Current cost value: ", Dynamic[oldC]];
While[Norm[δ] > ε,
newβ[[All, 2]] += δ;
- newr = r /. newβ;
+ 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, r, h, β];
+ a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β, Δ];
δ = v + a;
newβ = β;
newβ[[All, 2]] += δ;
- newr = r /. newβ;
+ newr = rf[Append[newβ[[All, 2]], Δ]];
newC = Re[Total[newr^2]];
];
oldC = newC;
oldr = newr;
β = newβ;
- oldJ = J /. β;
- U = First[SingularValueDecomposition[oldJ]];
+ 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, r, h, β];
+ a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β, Δ];
δ = v + a;
- nSteps += 1;
];
{newC, β}
]
-annealFit[r_, Δ_, β0_, δ0_ : 10^-5, δ1_ : 10^-13, δδ_ : 100, λ0_ : 1, ν_ : 2, ε_: 10^-15] :=
-Module[ {lastfit = {1, β0}, dr, δ = δ0},
- dr = D[r, {First /@ β0}];
- While[δ >= δ1,
- lastfit = levenbergMarquardtHelper[r /. Δ -> δ, dr /. Δ -> δ, lastfit[[2]], λ0, ν, ε];
- Print[lastfit[[1]]];
- δ /= δδ;
+annealFit[r_, β0_, γ_ : 0, λ0_ : 1, ν_ : 2, μ_ : 3, ε_: 10^-14, 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;
];
- lastfit
+ {oldC, β}
]
fitCov[res_, s_] :=