Merge branch 'main' into mooncake-ext

Author: thorek1

src/algorithms/nonlinear_solver.jl

@@ -482,10 +482,16 @@ function newton(
                 new_residuals .= guess_update
             else
                 fact∇ = ℒ.lu!(∇, check = false)
-                if !ℒ.issuccess(fact∇)
-                    fact∇ = ℒ.qr(∇, ℒ.ColumnNorm())
+                try
+                    if !ℒ.issuccess(fact∇)
+                        fact∇ = ℒ.qr(∇, ℒ.ColumnNorm())
+                    end
+                    ℒ.ldiv!(fact∇, new_residuals)
+                catch
+                    rel_xtol_reached = typemax(T)
+                    new_residuals_norm = typemax(T)
+                    break
                 end
-                ℒ.ldiv!(fact∇, new_residuals)
             end
 
             guess_update_norm = ℒ.norm(new_residuals)
@@ -530,10 +536,16 @@ function newton(
             new_residuals .= guess_update
         else
             fact∇ = ℒ.lu!(∇, check = false)
-            if !ℒ.issuccess(fact∇)
-                fact∇ = ℒ.qr(∇, ℒ.ColumnNorm())
+            try
+                if !ℒ.issuccess(fact∇)
+                    fact∇ = ℒ.qr(∇, ℒ.ColumnNorm())
+                end
+                ℒ.ldiv!(fact∇, new_residuals)
+            catch
+                rel_xtol_reached = typemax(T)
+                new_residuals_norm = typemax(T)
+                break
             end
-            ℒ.ldiv!(fact∇, new_residuals)
         end
 
         guess_update_norm = ℒ.norm(new_residuals)

src/filter/find_shocks.jl

@@ -77,16 +77,14 @@ function find_shocks(::Val{:LagrangeNewton},
         #     return x, false
         # end
 
-        f̂xλp = try 
-            ℒ.factorize(fxλp)
+        try
+            f̂xλp = ℒ.factorize(fxλp)
+            ℒ.ldiv!(Δxλ, f̂xλp, fxλ)
         catch
             # ℒ.svd(fxλp)
             # println("factorization fails")
             return x, false
         end
-
-        # Δxλ = fxλp \ fxλ
-        ℒ.ldiv!(Δxλ, f̂xλp, fxλ)
 
         if !all(isfinite,Δxλ) break end
 
@@ -283,16 +281,14 @@ function find_shocks(::Val{:LagrangeNewton},
         # fXλp = [reshape((2 * 𝐒ⁱ²ᵉ + 6 * 𝐒ⁱ³ᵉ * ℒ.kron(ℒ.I(length(x)), ℒ.kron(ℒ.I(length(x)),x)))' * λ, size(𝐒ⁱ, 2), size(𝐒ⁱ, 2)) - 2*ℒ.I(size(𝐒ⁱ, 2))  (𝐒ⁱ + 2 * 𝐒ⁱ²ᵉ * ℒ.kron(ℒ.I(length(x)), x) + 3 * 𝐒ⁱ³ᵉ * ℒ.kron(ℒ.I(length(x)), ℒ.kron(x, x)))'
         #         -(𝐒ⁱ + 2 * 𝐒ⁱ²ᵉ * ℒ.kron(ℒ.I(length(x)), x) + 3 * 𝐒ⁱ³ᵉ * ℒ.kron(ℒ.I(length(x)), ℒ.kron(x, x)))  zeros(size(𝐒ⁱ, 1),size(𝐒ⁱ, 1))]
 
-        f̂xλp = try 
-            ℒ.factorize(fxλp)
+        try
+            f̂xλp = ℒ.factorize(fxλp)
+            ℒ.ldiv!(Δxλ, f̂xλp, fxλ)
         catch
             # ℒ.svd(fxλp)
             # println("factorization fails")
             return x, false
         end
-
-        # Δxλ = fxλp \ fxλ
-        ℒ.ldiv!(Δxλ, f̂xλp, fxλ)
 
         if !all(isfinite,Δxλ) break end
 

src/filter/inversion.jl

@@ -729,10 +729,14 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood),
 
         copy!(jacct, jacc[i]')
 
-        jacc_fact = ℒ.factorize(jacct) # otherwise this fails for nshocks > nexo
-
-        # λ[i] = jacc[i]' \ x[i] * 2
-        ℒ.ldiv!(λ[i], jacc_fact, x[i])
+        try
+            jacc_fact = ℒ.factorize(jacct) # otherwise this fails for nshocks > nexo
+            # λ[i] = jacc[i]' \ x[i] * 2
+            ℒ.ldiv!(λ[i], jacc_fact, x[i])
+        catch
+            if opts.verbose println("Inversion filter failed at step $i") end
+            return -Inf, x -> NoTangent(), NoTangent(),  NoTangent(), NoTangent(), NoTangent(), NoTangent(),  NoTangent(),  NoTangent(),  NoTangent(), NoTangent()
+        end
 
         ℒ.rmul!(λ[i], 2)
 
@@ -1410,10 +1414,14 @@ function rrule(::typeof(calculate_inversion_filter_loglikelihood),
 
         copy!(jacct, jacc[i]')
 
-        jacc_fact = ℒ.factorize(jacct)
-
-        # λ[i] = jacc[i]' \ x[i] * 2
-        ℒ.ldiv!(λ[i], jacc_fact, x[i])
+        try
+            jacc_fact = ℒ.factorize(jacct)
+            # λ[i] = jacc[i]' \ x[i] * 2
+            ℒ.ldiv!(λ[i], jacc_fact, x[i])
+        catch
+            if opts.verbose println("Inversion filter failed at step $i") end
+            return -Inf, x -> NoTangent(), NoTangent(),  NoTangent(), NoTangent(), NoTangent(), NoTangent(),  NoTangent(),  NoTangent(),  NoTangent(), NoTangent()
+        end
 
         # ℒ.ldiv!(λ[i], jacc_fact', x[i])
         ℒ.rmul!(λ[i], 2)