diff options
Diffstat (limited to 'schofield.wl')
-rw-r--r-- | schofield.wl | 109 |
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_] := |