diff options
Diffstat (limited to 'schofield.wl')
| -rw-r--r-- | schofield.wl | 60 |
1 files changed, 30 insertions, 30 deletions
diff --git a/schofield.wl b/schofield.wl index f32d855..3f45f5d 100644 --- a/schofield.wl +++ b/schofield.wl @@ -96,9 +96,11 @@ dΦdηList[n_, h_][m_, θp_] := Module[ ] ruleB[g_] := B -> (2 * OverBar[s] / π) * (- g'[θc] / t[θc]^Δ[2]) -ruleθ0[g_] := Simplify[g[I θ0]/(-t[I θ0])^Δ[2]/I] - Around[0.18930, 0.00005] +ξ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] -ruleAH[g_] := AH / ((g[I θ0]/ I)^(2 / Δ[2]) * (-η[g]'[I θ0] / (2 θ0 I))^(5/6)) + Around[1.37, 0.02] +ruleAH[g_] := AYL[g] + Around[1.37, 0.02] eqLowRHSReg[n_][m_] := dRFc[n][m] @@ -144,9 +146,9 @@ Gls = { Around[ 2.937665, 3.0 10^(-6)], Around[-15.61, 1.0 10^(-2)], Around[ 96.76, 1.0 10^(-2)], - -6.79 10^2, - 5.34 10^3, - -4.66 10^4, + Around[-6.79 10^2, 1.0], + Around[5.34 10^3, 10.], + Around[-4.66 10^4, 0.01 10^4], 4.46 10^5, -4.66 10^6 } @@ -187,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], {ruleθ0[g] / 0.00005(*, ruleAH[g] / 0.02*)} /. Around[x_, _] :> x] +res[F_, g_, δ_][m_] := Join[resLow[F, g, δ][m], resHigh[F, g, δ][m]] chiSquared[F_, g_, δ_][m_] := Total[res[F, g, δ][m]^2] newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[ @@ -205,8 +207,9 @@ levenbergMarquardtStep[M_, λ_, g_] := Quiet[ , LinearSolve::luc ] -levenbergMarquardtAcceleration[M_, λ_, v_, J_, rold_, rf_, h_, β_, δ_] := Quiet[ - - Re[LinearSolve[M + λ DiagonalMatrix[Diagonal[M]], Transpose[J].(2/h((rf[Append[β[[All,2]] + h v, δ]] - rold) / h - J . v))] / 2] +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] , LinearSolve::luc ] @@ -229,7 +232,7 @@ levenbergMarquardt[r_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] := ] -levenbergMarquardtHelper[Δ_, rf_, Jf_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15, accOff_ : False] := +levenbergMarquardtHelper[r_, J_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] := Module[ { n = Length[β0], @@ -240,66 +243,63 @@ levenbergMarquardtHelper[Δ_, rf_, Jf_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15, a h = 0.1, fvv, α = 0.1, a, v, U, nSteps = 0, reject, μ =3 }, PrintTemporary["Beginning the algorithm."]; - oldr = rf[Append[β[[All, 2]], Δ]]; - oldJ = Jf[Append[β[[All, 2]], Δ]]; + oldr = r /. β; + oldJ = J /. β; U = First[SingularValueDecomposition[oldJ]]; 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 + If[accOff, 0, a]; + a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, r, h, β]; + δ = v + a; vlast = v; PrintTemporary["Current cost value: ", Dynamic[oldC]]; While[Norm[δ] > ε, newβ[[All, 2]] += δ; - newr = rf[Append[newβ[[All, 2]], Δ]]; + newr = r /. newβ; 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 + If[accOff, 0, a]; + a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, r, h, β]; + δ = v + a; newβ = β; newβ[[All, 2]] += δ; - newr = rf[Append[newβ[[All, 2]], Δ]]; + newr = r /. newβ; newC = Re[Total[newr^2]]; ]; oldC = newC; oldr = newr; β = newβ; - oldJ = Jf[Append[β[[All, 2]], Δ]]; + oldJ = J /. β; U = First[SingularValueDecomposition[oldJ]]; 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 + If[accOff, 0, a]; + a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, r, 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}, rf, Jf, δ = δ0, n = Length[β0]}, - Quiet[ - rf = Compile[{{x, _Real, 1}}, - Evaluate[r /. Thread[Rule[Append[First /@ β0, Δ], Part[x, #] & /@ Range[n + 1]]]] - ]; - Jf = Compile[{{x, _Real, 1}}, - Evaluate[D[r, {First /@ β0}] /. Thread[Rule[Append[First /@ β0, Δ], Part[x, #] & /@ Range[n + 1]]]] - ]; - , Part::partd]; +Module[ {lastfit = {1, β0}, dr, δ = δ0}, + dr = D[r, {First /@ β0}]; While[δ >= δ1, - lastfit = levenbergMarquardtHelper[δ, rf, Jf, lastfit[[2]]]; + lastfit = levenbergMarquardtHelper[r /. Δ -> δ, dr /. Δ -> δ, lastfit[[2]], λ0, ν, ε]; Print[lastfit[[1]]]; δ /= δδ; ]; lastfit ] +fitCov[res_, s_] := + With[{J = D[res, {s[[All, 1]]}]}, + Symmetrize[Inverse[SetPrecision[Transpose[J] . J /. s, 60]]]] + EndPackage[] |