More changes.

1 files changed, 26 insertions, 15 deletions

diff --git a/schofield.wl b/schofield.wl
index eb0525c..908ca94 100644
--- a/schofield.wl
+++ b/schofield.wl
@@ -208,40 +208,48 @@ newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[
   opts
 ]
 
-levenbergMarquardtStep[M_, λ_, g_] := LinearSolve[M + λ DiagonalMatrix[Diagonal[M]], g]
+levenbergMarquardtStep[M_, λ_, g_] := Quiet[LinearSolve[M + λ DiagonalMatrix[Diagonal[M]], g], LinearSolve::luc]
 
-levenburgMarquardt[r_, β0_, λ0_ : 1, ν_ : 2] := 
+levenburgMarquardt[r_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] := 
  Module[
   {
     n = Length[β0],
     β = β0,
     λ = λ0,
     newβ = β0,
-    x, rf, Jf, oldJ, oldr, newr, M, g, δ, oldC, newC
+    x, rf, Jf, oldJ, oldr, newr, M, g, δ, oldC, newC,
+    h = 0.1, fvv, α, a, v
    },
   PrintTemporary["Compiling the Jacobian..."];
-  rf = Compile[{{x, _Real, 1}}, 
+  rf = Quiet[Compile[{{x, _Real, 1}}, 
     Evaluate[r /. Thread[Rule[First /@ β, Part[x, #] & /@ Range[n]]]]
-  ];
-  Jf = Compile[{{x, _Real, 1}}, 
-    Evaluate[D[r, {First /@ β}] /. 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]]];
   oldC = Re[Total[oldr^2]];
+  newC = oldC;
   g = Re[Transpose[oldJ] . oldr];
   M = Re[Transpose[oldJ] . oldJ];
-  δ = levenbergMarquardtStep[M, λ / ν, g];
-  PrintTemporary[Dynamic[oldC]]
-  While[Norm[δ] > 10^-15,
-    newβ[[All, 2]] -= δ;
+  v = levenbergMarquardtStep[M, λ / ν, g];
+  a = LinearSolve[oldJ, -2/h((rf[β[[All,2]] + h v] - oldr) / h - oldJ . v)];
+  If[2 Norm[a] / Norm[v] > α, a /= (2 / Norm[v])];
+  δ = v + a / 2;
+  PrintTemporary["Current cost value: ", Dynamic[oldC]]
+  While[Norm[δ] > ε,
+    newβ[[All, 2]] += δ;
     newr = rf[newβ[[All, 2]]];
     newC = Re[Total[newr^2]];
     While[newC > oldC,
-      δ = levenbergMarquardtStep[M, λ, g];
+      v = levenbergMarquardtStep[M, λ, g];
+      a = LinearSolve[oldJ, -2/h((rf[β[[All,2]] + h v] - oldr) / h - oldJ . v)];
+      If[2 Norm[a] / Norm[v] > α, a /= (2 / Norm[v])];
+      δ = v + a / 2;
       newβ = β;
-      newβ[[All, 2]] -= δ;
+      newβ[[All, 2]] += δ;
       newr = rf[newβ[[All, 2]]];
       newC = Re[Total[newr^2]];
       λ *= ν;
@@ -253,7 +261,10 @@ levenburgMarquardt[r_, β0_, λ0_ : 1, ν_ : 2] :=
     oldJ = Jf[β[[All, 2]]];
     g = Re[Transpose[oldJ] . oldr];
     M = Re[Transpose[oldJ] . oldJ];
-    δ = levenbergMarquardtStep[M, λ / ν, g];
+    v = levenbergMarquardtStep[M, λ / ν, g];
+    a = LinearSolve[oldJ, -2/h((rf[β[[All,2]] + h v] - oldr) / h - oldJ . v)];
+    If[2 Norm[a] / Norm[v] > α, a /= (2 / Norm[v])];
+    δ = v + a / 2;
   ];
   {newC, β}
  ]