Work.

1 files changed, 59 insertions, 31 deletions

diff --git a/schofield.wl b/schofield.wl
index 77dce96..40253d4 100644
--- a/schofield.wl
+++ b/schofield.wl
@@ -213,73 +213,101 @@ 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[β[[All,2]] + h v] - rold) / h - J . v))] / 2]
+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]
   , 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, ν, ε]
+ ]
+
 
-levenburgMarquardt[r_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] := 
+levenbergMarquardtHelper[Δ_, rf_, Jf_, β0_, λ0_ : 1, ν_ : 1.5, ε_ : 10^-15, accOff_ : False] := 
  Module[
   {
     n = Length[β0],
     β = β0,
     λ = λ0,
     newβ = β0,
-    x, rf, Jf, oldJ, oldr, newr, M, g, δ, oldC, newC, δlast,
-    h = 0.1, fvv, α = 0.75, a, v, U, nSteps = 0
+    x, oldJ, oldr, newr, M, g, δ, oldC, newC, vlast,
+    h = 0.1, fvv, α = 0.75, a, v, U, nSteps = 0, reject, μ =3
    },
-  PrintTemporary["Compiling the Jacobian..."];
-  rf = Quiet[Compile[{{x, _Real, 1}}, 
-    Evaluate[r /. Thread[Rule[First /@ β, Part[x, #] & /@ Range[n]]]]
-  ], Part::partd];
-  Jf = Quiet[Compile[{{x, _Real, 1}}, 
-    Evaluate[D[r, {First /@ β}] /. Thread[Rule[First /@ β, Part[x, #] & /@ Range[n]]]]
-  ], Part::partd];
   PrintTemporary["Beginning the algorithm."];
-  oldr = rf[β[[All, 2]]];
-  oldJ = Jf[β[[All, 2]]];
+  oldr = rf[Append[β[[All, 2]], Δ]];
+  oldJ = Jf[Append[β[[All, 2]], Δ]];
   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 + a;
-  δlast = δ;
+  v = levenbergMarquardtStep[M, λ, g];
+  a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β, Δ];
+  δ = v + If[accOff, 0, a];
+  vlast = v;
   PrintTemporary["Current cost value: ", Dynamic[oldC]];
   While[Norm[δ] > ε,
     newβ[[All, 2]] += δ;
-    newr = rf[newβ[[All, 2]]];
+    newr = rf[Append[newβ[[All, 2]], Δ]];
     newC = Re[Total[newr^2]];
-    While[((1 - δ . δlast / Norm[δ] / Norm[δlast])newC > oldC) || (2 Norm[a] / Norm[v] > α),
+    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;
+      a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β, Δ];
+      δ = v + If[accOff, 0, a];
       newβ = β;
       newβ[[All, 2]] += δ;
-      newr = rf[newβ[[All, 2]]];
+      newr = rf[Append[newβ[[All, 2]], Δ]];
       newC = Re[Total[newr^2]];
-      λ *= ν;
     ];
-    λ /= ν;
     oldC = newC;
     oldr = newr;
     β = newβ;
-    oldJ = Jf[β[[All, 2]]];
+    oldJ = Jf[Append[β[[All, 2]], Δ]];
     U = First[SingularValueDecomposition[oldJ]];
     g = Re[Transpose[oldJ] . oldr];
     M = Re[Transpose[oldJ] . oldJ];
-    v = levenbergMarquardtStep[M, λ / ν, g];
-    a = levenbergMarquardtAcceleration[M, λ / ν, v, oldJ, oldr, rf, h, β];
-    δlast = δ;
-    δ = v + a;
+    vlast = v;
+    λ /= μ;
+    v = levenbergMarquardtStep[M, λ, g];
+    a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β, Δ];
+    δ = v + If[accOff, 0, a];
     nSteps += 1;
   ];
-  Print[nSteps];
   {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];
+  While[δ >= δ1,
+    lastfit = levenbergMarquardtHelper[δ, rf, Jf, lastfit[[2]]];
+    Print[lastfit[[1]]];
+    δ /= δδ;
+  ];
+  lastfit
+]
+
 EndPackage[]