From 60d605279a658bc5de8d712417676c8fb93927e8 Mon Sep 17 00:00:00 2001
From: John Cobb <jdcobb3@gmail.com>
Date: Thu, 7 May 2026 10:07:33 -0400
Subject: [PATCH 1/3] Add start-solution expansion controls & interrupts

Introduce configurable start-solution expansion and interrupt handling. Added StartSolutionExpansionResult type and iterator helpers; added flags expand_start_solutions, expand_start_solutions_newton, expand_start_solutions_gradient_flow and catch_interrupt to critical_points/_expand_start_solutions/_solve_and_trace to allow selective expansion and safe interruption. Improve robustness by catching interrupt exceptions during Newton and gradient-flow phases and returning partial results when appropriate. Wrap graph connection loops with interrupt-aware try/catch and add catch_interrupt arg. Update tests to exercise new expansion options and adjust Project.toml extras/targets for testing.
---
 Project.toml           |   9 ++
 src/critical_points.jl | 222 +++++++++++++++++++++++++++++------------
 src/graph.jl           | 154 ++++++++++++++--------------
 test/runtests.jl       |  35 ++++++-
 4 files changed, 280 insertions(+), 140 deletions(-)

diff --git a/Project.toml b/Project.toml
index 2ac792c..06e10d7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,3 +23,12 @@ Plots = "1.41.6"
 ProgressMeter = "1.11"
 Reexport = "1"
 SciMLBase = "2"
+
+[extras]
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Logging", "Pkg", "Random", "Test"]
diff --git a/src/critical_points.jl b/src/critical_points.jl
index c6d702d..f3f3a14 100644
--- a/src/critical_points.jl
+++ b/src/critical_points.jl
@@ -3,10 +3,38 @@ export critical_points
 
 import HomotopyContinuation: MonodromyOptions, UniquePoints, EndgameTracker
 
+struct StartSolutionExpansionResult{TS,TN}
+    start_solutions::TS
+    routing_points::TN
+    interrupted::Bool
+end
+
+function Base.iterate(result::StartSolutionExpansionResult, state = 1)
+    if state == 1
+        return result.start_solutions, 2
+    elseif state == 2
+        return result.routing_points, 3
+    else
+        return nothing
+    end
+end
+
+Base.length(::StartSolutionExpansionResult) = 2
+
+_is_interrupt_exception(e) =
+    e isa InterruptException ||
+    (e isa TaskFailedException && _is_interrupt_exception(e.task.exception))
+
 """
     critical_points(r, S0, rhs0; kwargs...)
 
 Find critical points of the routing function using monodromy and gradient flow.
+
+The start-solution expansion can be disabled with `expand_start_solutions = false`.
+When expansion is enabled, the Newton and gradient-flow substeps can be controlled
+independently with `expand_start_solutions_newton` and
+`expand_start_solutions_gradient_flow`. If `catch_interrupt = true`, interrupting a
+long-running phase returns the solutions found so far when possible.
 """
 function critical_points(
     r::RoutingFunction,
@@ -16,6 +44,10 @@ function critical_points(
     start_grid_width = 5,
     start_grid_stepsize = 0.2,
     start_grid_center = nothing,
+    expand_start_solutions = true,
+    expand_start_solutions_newton = true,
+    expand_start_solutions_gradient_flow = true,
+    catch_interrupt = true,
     monodromy_at_zero = false,
     options = MonodromyOptions(
         parameter_sampler = p -> 10 .* randn(ComplexF64, length(p)),
@@ -35,20 +67,31 @@ function critical_points(
     )
 
     # Step 2: Expand start solutions via Newton's method and gradient flow
-    S0, new_pts = _expand_start_solutions(
+    expansion_result = _expand_start_solutions(
         ∇r, H, S0, rhs0, k;
         verbose = verbose,
         start_grid_width = start_grid_width,
         start_grid_stepsize = start_grid_stepsize,
         start_grid_center = start_grid_center,
+        expand_start_solutions = expand_start_solutions,
+        expand_start_solutions_newton = expand_start_solutions_newton,
+        expand_start_solutions_gradient_flow = expand_start_solutions_gradient_flow,
+        catch_interrupt = catch_interrupt,
         monodromy_at_zero = monodromy_at_zero,
     )
+    S0, new_pts = expansion_result
+    if expansion_result.interrupted
+        verbose && @warn "Interrupted while expanding start solutions. Returning routing points found so far."
+        return real.(new_pts), nothing, nothing
+    end
 
     # Step 3: Solve and trace to critical points
     return _solve_and_trace(
         MS, H, S0, rhs0, new_pts;
         monodromy_at_zero = monodromy_at_zero,
-        start_grid_width = start_grid_width,
+        expand_start_solutions = expand_start_solutions && start_grid_width > 0,
+        catch_interrupt = catch_interrupt,
+        seed = seed,
     )
 end
 
@@ -128,12 +171,17 @@ function _expand_start_solutions(
     start_grid_width = 5,
     start_grid_stepsize = 0.2,
     start_grid_center = nothing,
+    expand_start_solutions = true,
+    expand_start_solutions_newton = true,
+    expand_start_solutions_gradient_flow = true,
+    catch_interrupt = true,
     monodromy_at_zero = false,
 )
     new_pts = Vector{ComplexF64}[]
     # Setting up grid
-    if start_grid_width <= 0
-        return S0, new_pts
+    if !expand_start_solutions || start_grid_width <= 0 ||
+       (!expand_start_solutions_newton && !expand_start_solutions_gradient_flow)
+        return StartSolutionExpansionResult(S0, new_pts, false)
     end
 
     if isnothing(start_grid_center)
@@ -141,96 +189,127 @@ function _expand_start_solutions(
     end
 
     w = (start_grid_width / 2)
-    grid = [
-        (start_grid_center[i]-w):start_grid_stepsize:(start_grid_center[i]+w) for
-        i = 1:k
-    ]
-    newton_w = 10*w
-    newton_grid = [
-        (start_grid_center[i]-newton_w):start_grid_stepsize:(start_grid_center[i]+newton_w) for
-        i = 1:k
-    ]
+    interrupted = false
 
     # First we try to find start solutions via blindly applying Newton's method to ∇r=0.
-    verbose && println("Expanding start solutions via Newton's method...")
+    newton_pts = Vector{ComplexF64}[]
     newton_success_count = 0
-    start_pt = zeros(ComplexF64, k)
-    ProgressMeter.@showprogress for start_point in Iterators.product(newton_grid...)
+    newton_total_count = 0
+    if expand_start_solutions_newton
+        newton_w = 10*w
+        newton_grid = [
+            (start_grid_center[i]-newton_w):start_grid_stepsize:(start_grid_center[i]+newton_w) for
+            i = 1:k
+        ]
+        newton_total_count = prod(length.(newton_grid))
+        verbose && println("Expanding start solutions via Newton's method...")
+        start_pt = zeros(ComplexF64, k)
+        ProgressMeter.@showprogress for start_point in Iterators.product(newton_grid...)
             start_pt .= start_point # this avoids allocations from splatting the tuple into the newton function
             # Newton's method on each initial guess
             try
                 pt = newton(∇r, start_pt, rhs0; max_iters = 200) |> solution
                 if norm(evaluate(∇r, pt, rhs0)) < 1e-10
                     newton_success_count += 1
-                    push!(new_pts, pt)
+                    push!(newton_pts, pt)
                 end
             catch e
+                if _is_interrupt_exception(e)
+                    catch_interrupt || rethrow(e)
+                    interrupted = true
+                    break
+                end
                 continue
             end
-
+        end
     end
 
     num_newton_pts = 0
-    if length(new_pts) > 0
-        new_pts = HC.unique_points(new_pts)
-        num_newton_pts += length(new_pts)
+    if length(newton_pts) > 0
+        newton_pts = HC.unique_points(newton_pts)
+        num_newton_pts = length(newton_pts)
     end
 
     # since the points obtained via Newton's method are already solutions to ∇r=rhs0, we can directly add them to S0 if monodromy_at_zero is false
     if !monodromy_at_zero
-        S0 = HC.unique_points([S0; new_pts])
-        # to avoid redundant work later, we remove the points found via Newton's method from new_pts so that we aren't tracing them using monodromy
-        empty!(new_pts)
+        S0 = HC.unique_points([S0; newton_pts])
+    elseif !isempty(newton_pts)
+        new_pts = HC.unique_points([new_pts; newton_pts])
     end
 
-    verbose && println("Successful Newton's method attempts: $(newton_success_count) out of $(length(newton_grid[1])^k) ($(round(newton_success_count / (length(newton_grid[1])^k) * 100, digits=2))%)")
-    verbose && println("Found $num_newton_pts solutions to ∇r(z)=rhs0.")
+    if expand_start_solutions_newton
+        verbose && println("Successful Newton's method attempts: $(newton_success_count) out of $(newton_total_count) ($(round(newton_success_count / newton_total_count * 100, digits=2))%)")
+        verbose && println("Found $num_newton_pts solutions to ∇r(z)=rhs0.")
+    end
+
+    if interrupted
+        return StartSolutionExpansionResult(S0, new_pts, true)
+    end
     
     # Now we try gradient flow
-    g(x, param, t) = real(evaluate(∇r, x))
-    tspan = (0.0, 1e4)
-    
-    verbose && println("Expanding the set of start solutions via gradient flow...")
-
     gradient_success_count = 0
-    start_pt = zeros(k)
-    ProgressMeter.@showprogress for start_point in Iterators.product(grid...)
-        try
-            start_pt .= start_point # this avoids allocations from splatting the tuple into the ODEProblem
-            prob = SciMLBase.ODEProblem(g, start_pt, tspan)
-            sol = DE.solve(prob, reltol = 1e-6, abstol = 1e-6)
-            convergence_point = last(sol.u)
-            improved_point = newton(∇r, convergence_point) |> solution
-            push!(new_pts, improved_point)
-            gradient_success_count += 1                
-        catch e
-            continue
+    gradient_total_count = 0
+    if expand_start_solutions_gradient_flow
+        grid = [
+            (start_grid_center[i]-w):start_grid_stepsize:(start_grid_center[i]+w) for
+            i = 1:k
+        ]
+        gradient_total_count = prod(length.(grid))
+        g(x, param, t) = real(evaluate(∇r, x))
+        tspan = (0.0, 1e4)
+
+        verbose && println("Expanding the set of start solutions via gradient flow...")
+
+        start_pt = zeros(k)
+        ProgressMeter.@showprogress for start_point in Iterators.product(grid...)
+            try
+                start_pt .= start_point # this avoids allocations from splatting the tuple into the ODEProblem
+                prob = SciMLBase.ODEProblem(g, start_pt, tspan)
+                sol = DE.solve(prob, reltol = 1e-6, abstol = 1e-6)
+                convergence_point = last(sol.u)
+                improved_point = newton(∇r, convergence_point) |> solution
+                push!(new_pts, improved_point)
+                gradient_success_count += 1
+            catch e
+                if _is_interrupt_exception(e)
+                    catch_interrupt || rethrow(e)
+                    interrupted = true
+                    break
+                end
+                continue
+            end
+        end
+        if length(new_pts) > 0
+            new_pts = HC.unique_points(new_pts)
         end
+        verbose && println("Successful gradient flow attempts: $(gradient_success_count) out of $(gradient_total_count) ($(round(gradient_success_count / gradient_total_count * 100, digits=2))%)")
+        verbose && println("Found $(length(new_pts)) routing points via gradient flow.")
     end
-    if length(new_pts) > 0
-        new_pts = HC.unique_points(new_pts)
+
+    if interrupted
+        return StartSolutionExpansionResult(S0, new_pts, true)
     end
-    verbose && println("Successful gradient flow attempts: $(gradient_success_count) out of $(length(grid[1])^k) ($(round(gradient_success_count / (length(grid[1])^k) * 100, digits=2))%)")
-    verbose && println("Found $(length(new_pts)) routing points via gradient flow.")
 
     if !monodromy_at_zero
-        start_parameters!(H, zeros(ComplexF64, length(rhs0)))
-        target_parameters!(H, rhs0)
-        S0_new_sols = HC.solve(H, new_pts) |> solutions
-        number_of_old_sols = length(S0)
-        S0 = HC.unique_points([S0; S0_new_sols])
-        verbose && println(
-            "Traced to $(length(S0)-number_of_old_sols) additional start solutions for the monodromy.",
-        )
+        if !isempty(new_pts)
+            start_parameters!(H, zeros(ComplexF64, length(rhs0)))
+            target_parameters!(H, rhs0)
+            S0_new_sols = HC.solve(H, new_pts; catch_interrupt = catch_interrupt) |> solutions
+            number_of_old_sols = length(S0)
+            S0 = HC.unique_points([S0; S0_new_sols])
+            verbose && println(
+                "Traced to $(length(S0)-number_of_old_sols) additional start solutions for the monodromy.",
+            )
+        end
     else
         S0 = [S0; new_pts]
     end
 
-    return S0, new_pts
+    return StartSolutionExpansionResult(S0, new_pts, false)
 end
 
 """
-    _solve_and_trace(MS, H, S0, rhs0, new_pts; monodromy_at_zero, start_grid_width)
+    _solve_and_trace(MS, H, S0, rhs0, new_pts; monodromy_at_zero, expand_start_solutions)
 
 Perform monodromy solving and trace solutions to ∇r=0.
 Returns (routing_points, result, mon_result).
@@ -242,24 +321,43 @@ function _solve_and_trace(
     rhs0::AbstractVector{<:Number},
     new_pts::AbstractVector{<:AbstractVector{<:Number}};
     monodromy_at_zero = false,
-    start_grid_width = 5,
+    expand_start_solutions = true,
+    catch_interrupt = true,
+    seed = rand(UInt32),
 )
     ### Monodromy
-    mon_result = monodromy_solve(MS, S0, rhs0, rand(UInt32);)
+    mon_result = monodromy_solve(MS, S0, rhs0, seed; catch_interrupt = catch_interrupt)
 
     ### Trace to ∇r=0
     if !monodromy_at_zero
+        if isempty(solutions(mon_result))
+            routing_points = expand_start_solutions ? real.(new_pts) : Vector{Float64}[]
+            return routing_points, nothing, mon_result
+        end
+
         intermediate_rhs = randn(ComplexF64, length(rhs0))
         start_parameters!(H, rhs0)
         target_parameters!(H, intermediate_rhs)
-        result_intermediate = HomotopyContinuation.solve(H, solutions(mon_result))
+        result_intermediate = HomotopyContinuation.solve(
+            H,
+            solutions(mon_result);
+            catch_interrupt = catch_interrupt,
+        )
+        if isempty(solutions(result_intermediate))
+            routing_points = expand_start_solutions ? real.(new_pts) : Vector{Float64}[]
+            return routing_points, result_intermediate, mon_result
+        end
         start_parameters!(H, intermediate_rhs)
         target_parameters!(H, zeros(ComplexF64, length(rhs0)))
-        result = HomotopyContinuation.solve(H, result_intermediate)
+        result = HomotopyContinuation.solve(
+            H,
+            result_intermediate;
+            catch_interrupt = catch_interrupt,
+        )
         routing_points = real_solutions(result)
 
         # Make sure none of the routing points found via gradient flow are lost
-        if start_grid_width > 0
+        if expand_start_solutions
             routing_points = HC.unique_points([routing_points; real.(new_pts)])
         end
         return routing_points, result, mon_result
diff --git a/src/graph.jl b/src/graph.jl
index 3f137a1..bec5f60 100644
--- a/src/graph.jl
+++ b/src/graph.jl
@@ -67,7 +67,8 @@ function partition_of_critical_points(
     crit_pts::Vector{Vector{Float64}},
     epsilon::Float64 = 1e-6,
     reltol::Float64 = 1e-6,
-    abstol::Float64 = 1e-9,
+    abstol::Float64 = 1e-9;
+    catch_interrupt = true,
 )
 
     ∇r = RoutingGradient(r)
@@ -95,108 +96,113 @@ function partition_of_critical_points(
 
 
     failed_info_list = []
-    ProgressMeter.@showprogress for i = 1:length(index_list)
-        if connectivity_status[i] == 0 && index_list[i] == 1
-            # need to do path tracking in two directions
-            critical_point_index = critical_points_indices[i]
-            unstable_eigenvector = unstable_eigenvector_list[critical_point_index]
-            pair_pos, failed_info_pos = limit_critical_point_from_critical_point(
-                ode_log!,
-                crit_pts,
-                i,
-                index_list,
-                unstable_eigenvector,
-                epsilon,
-                reltol,
-                abstol,
-            )
-            if isempty(failed_info_pos)
-                LightGraphs.add_edge!(graph, pair_pos[1], pair_pos[2])
-            else
-                push!(
-                    failed_info_list,
-                    [critical_point_index, epsilon, unstable_eigenvector, failed_info_pos],
-                )
-            end
-
-            pair_neg, failed_info_neg = limit_critical_point_from_critical_point(
-                ode_log!,
-                crit_pts,
-                i,
-                index_list,
-                unstable_eigenvector,
-                -epsilon,
-                reltol,
-                abstol,
-            )
-
-            if isempty(failed_info_neg)
-                LightGraphs.add_edge!(graph, pair_neg[1], pair_neg[2])
-            else
-                push!(
-                    failed_info_list,
-                    [critical_point_index, -epsilon, unstable_eigenvector, failed_info_neg],
-                )
-            end
-            connectivity_status[i] = 1
-        end
-    end
-
-    ProgressMeter.@showprogress for i = 1:length(index_list)
-        if connectivity_status[i] == 0 && index_list[i] > 1
-            # track paths and stop whenever one path converges
-            critical_point_index = critical_points_indices[i]
-            sub_failed_info_list = []
-            for v in eachcol(unstable_eigenvector_list[critical_point_index])
-                pair, failed_info = limit_critical_point_from_critical_point(
+    try
+        ProgressMeter.@showprogress for i = 1:length(index_list)
+            if connectivity_status[i] == 0 && index_list[i] == 1
+                # need to do path tracking in two directions
+                critical_point_index = critical_points_indices[i]
+                unstable_eigenvector = unstable_eigenvector_list[critical_point_index]
+                pair_pos, failed_info_pos = limit_critical_point_from_critical_point(
                     ode_log!,
                     crit_pts,
                     i,
                     index_list,
-                    v,
+                    unstable_eigenvector,
                     epsilon,
                     reltol,
                     abstol,
                 )
-                if isempty(failed_info)
-                    LightGraphs.add_edge!(graph, pair[1], pair[2])
-                    connectivity_status[i] = 1
-                    sub_failed_info_list = []
-                    break
+                if isempty(failed_info_pos)
+                    LightGraphs.add_edge!(graph, pair_pos[1], pair_pos[2])
                 else
                     push!(
-                        sub_failed_info_list,
-                        [critical_point_index, epsilon, v, failed_info],
+                        failed_info_list,
+                        [critical_point_index, epsilon, unstable_eigenvector, failed_info_pos],
                     )
                 end
 
-                pair, failed_info = limit_critical_point_from_critical_point(
+                pair_neg, failed_info_neg = limit_critical_point_from_critical_point(
                     ode_log!,
                     crit_pts,
                     i,
                     index_list,
-                    v,
+                    unstable_eigenvector,
                     -epsilon,
                     reltol,
                     abstol,
                 )
-                if isempty(failed_info)
-                    LightGraphs.add_edge!(graph, pair[1], pair[2])
-                    connectivity_status[i] = 1
-                    sub_failed_info_list = []
-                    break
+
+                if isempty(failed_info_neg)
+                    LightGraphs.add_edge!(graph, pair_neg[1], pair_neg[2])
                 else
                     push!(
-                        sub_failed_info_list,
-                        [critical_point_index, epsilon, v, failed_info],
+                        failed_info_list,
+                        [critical_point_index, -epsilon, unstable_eigenvector, failed_info_neg],
                     )
                 end
+                connectivity_status[i] = 1
             end
+        end
 
-            if !isempty(sub_failed_info_list)
-                push!(failed_info_list, sub_failed_info_list)
+        ProgressMeter.@showprogress for i = 1:length(index_list)
+            if connectivity_status[i] == 0 && index_list[i] > 1
+                # track paths and stop whenever one path converges
+                critical_point_index = critical_points_indices[i]
+                sub_failed_info_list = []
+                for v in eachcol(unstable_eigenvector_list[critical_point_index])
+                    pair, failed_info = limit_critical_point_from_critical_point(
+                        ode_log!,
+                        crit_pts,
+                        i,
+                        index_list,
+                        v,
+                        epsilon,
+                        reltol,
+                        abstol,
+                    )
+                    if isempty(failed_info)
+                        LightGraphs.add_edge!(graph, pair[1], pair[2])
+                        connectivity_status[i] = 1
+                        sub_failed_info_list = []
+                        break
+                    else
+                        push!(
+                            sub_failed_info_list,
+                            [critical_point_index, epsilon, v, failed_info],
+                        )
+                    end
+
+                    pair, failed_info = limit_critical_point_from_critical_point(
+                        ode_log!,
+                        crit_pts,
+                        i,
+                        index_list,
+                        v,
+                        -epsilon,
+                        reltol,
+                        abstol,
+                    )
+                    if isempty(failed_info)
+                        LightGraphs.add_edge!(graph, pair[1], pair[2])
+                        connectivity_status[i] = 1
+                        sub_failed_info_list = []
+                        break
+                    else
+                        push!(
+                            sub_failed_info_list,
+                            [critical_point_index, epsilon, v, failed_info],
+                        )
+                    end
+                end
+
+                if !isempty(sub_failed_info_list)
+                    push!(failed_info_list, sub_failed_info_list)
+                end
             end
         end
+    catch e
+        catch_interrupt && _is_interrupt_exception(e) || rethrow(e)
+        @warn "Interrupted while connecting critical points. Returning the partial partition."
     end
 
 
diff --git a/test/runtests.jl b/test/runtests.jl
index cd6e54e..26f8870 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -89,9 +89,36 @@ end
 
     # Test that the expansion of start solutions works
     ∇r = RoutingGradient(r)
-    MS, H, S0, rhs0, k = ProjectedHypersurfaceRegions._setup_monodromy_solver(∇r)
+    MS, H, S0_initial, rhs0, k = ProjectedHypersurfaceRegions._setup_monodromy_solver(∇r)
+    S0_skip, new_pts_skip = ProjectedHypersurfaceRegions._expand_start_solutions(
+        ∇r, H, S0_initial, rhs0, k;
+        start_grid_width = 10,
+        start_grid_stepsize = 1,
+        expand_start_solutions = false,
+    )
+    @test S0_skip == S0_initial
+    @test isempty(new_pts_skip)
+
+    S0_no_gradient, new_pts_no_gradient = ProjectedHypersurfaceRegions._expand_start_solutions(
+        ∇r, H, S0_initial, rhs0, k;
+        start_grid_width = 1,
+        start_grid_stepsize = 1,
+        expand_start_solutions_gradient_flow = false,
+    )
+    @test all(norm(∇r(z)-rhs0) < 1e-10 for z in S0_no_gradient)
+    @test isempty(new_pts_no_gradient)
+
+    S0_no_newton, new_pts_no_newton = ProjectedHypersurfaceRegions._expand_start_solutions(
+        ∇r, H, S0_initial, rhs0, k;
+        start_grid_width = 1,
+        start_grid_stepsize = 1,
+        expand_start_solutions_newton = false,
+    )
+    @test all(norm(∇r(z)-rhs0) < 1e-10 for z in S0_no_newton)
+    @test all(norm(∇r(z)) < 1e-10 for z in new_pts_no_newton)
+
     S0, new_pts = ProjectedHypersurfaceRegions._expand_start_solutions(
-        ∇r, H, S0, rhs0, k;
+        ∇r, H, S0_initial, rhs0, k;
         start_grid_width = 10,
         start_grid_stepsize = 1,
     )
@@ -101,7 +128,7 @@ end
 
     # Check critical points
     options = MonodromyOptions(target_solutions_count = 2)
-    pts, res0, mon_res = critical_points(r, start_grid_width=0, options=options)    
+    pts, res0, mon_res = critical_points(r, expand_start_solutions=false, options=options)
     @test all(norm.(∇r_symbolic.(solutions(res0))) .< 1e-12)
 
     pl = generate_plot(
@@ -299,4 +326,4 @@ end
     [8*p[1]/(p[1]^2 - 4*p[2])^2  -16/(p[1]^2 - 4*p[2])^2]]
     @test Hess_log_abs_h(pt) - ProjectedHypersurfaceRegions.gradient_and_hessian(h, pt)[2] |> norm < 1e-6
 
-end
\ No newline at end of file
+end

From 1b87c438d6380a63a3839e2f5ff8a29c9ec0f7a3 Mon Sep 17 00:00:00 2001
From: John Cobb <jdcobb3@gmail.com>
Date: Thu, 7 May 2026 12:55:15 -0400
Subject: [PATCH 2/3] Fix interrupt handling and Newton expansion controls

---
 src/critical_points.jl | 44 ++++++++++++++++++++++++++++--------------
 test/runtests.jl       | 15 +++++++++++++-
 2 files changed, 44 insertions(+), 15 deletions(-)

diff --git a/src/critical_points.jl b/src/critical_points.jl
index f3f3a14..7c4b521 100644
--- a/src/critical_points.jl
+++ b/src/critical_points.jl
@@ -25,6 +25,13 @@ _is_interrupt_exception(e) =
     e isa InterruptException ||
     (e isa TaskFailedException && _is_interrupt_exception(e.task.exception))
 
+function _solve_homotopy(H, starts; catch_interrupt)
+    starts = collect(starts)
+    result = HomotopyContinuation.solve(H, starts; catch_interrupt = catch_interrupt)
+    interrupted = catch_interrupt && length(result) < length(starts)
+    return result, interrupted
+end
+
 """
     critical_points(r, S0, rhs0; kwargs...)
 
@@ -267,9 +274,17 @@ function _expand_start_solutions(
                 prob = SciMLBase.ODEProblem(g, start_pt, tspan)
                 sol = DE.solve(prob, reltol = 1e-6, abstol = 1e-6)
                 convergence_point = last(sol.u)
-                improved_point = newton(∇r, convergence_point) |> solution
-                push!(new_pts, improved_point)
-                gradient_success_count += 1
+                if expand_start_solutions_newton
+                    candidate_point = newton(∇r, convergence_point) |> solution
+                    residual_tolerance = 1e-10
+                else
+                    candidate_point = ComplexF64.(convergence_point)
+                    residual_tolerance = 1e-6
+                end
+                if norm(evaluate(∇r, candidate_point)) < residual_tolerance
+                    push!(new_pts, candidate_point)
+                    gradient_success_count += 1
+                end
             catch e
                 if _is_interrupt_exception(e)
                     catch_interrupt || rethrow(e)
@@ -294,12 +309,16 @@ function _expand_start_solutions(
         if !isempty(new_pts)
             start_parameters!(H, zeros(ComplexF64, length(rhs0)))
             target_parameters!(H, rhs0)
-            S0_new_sols = HC.solve(H, new_pts; catch_interrupt = catch_interrupt) |> solutions
+            S0_result, interrupted = _solve_homotopy(H, new_pts; catch_interrupt = catch_interrupt)
+            S0_new_sols = solutions(S0_result)
             number_of_old_sols = length(S0)
             S0 = HC.unique_points([S0; S0_new_sols])
             verbose && println(
                 "Traced to $(length(S0)-number_of_old_sols) additional start solutions for the monodromy.",
             )
+            if interrupted
+                return StartSolutionExpansionResult(S0, new_pts, true)
+            end
         end
     else
         S0 = [S0; new_pts]
@@ -330,28 +349,25 @@ function _solve_and_trace(
 
     ### Trace to ∇r=0
     if !monodromy_at_zero
-        if isempty(solutions(mon_result))
-            routing_points = expand_start_solutions ? real.(new_pts) : Vector{Float64}[]
-            return routing_points, nothing, mon_result
-        end
-
         intermediate_rhs = randn(ComplexF64, length(rhs0))
         start_parameters!(H, rhs0)
         target_parameters!(H, intermediate_rhs)
-        result_intermediate = HomotopyContinuation.solve(
+        monodromy_solutions = solutions(mon_result)
+        result_intermediate, interrupted = _solve_homotopy(
             H,
-            solutions(mon_result);
+            monodromy_solutions;
             catch_interrupt = catch_interrupt,
         )
-        if isempty(solutions(result_intermediate))
+        if interrupted
             routing_points = expand_start_solutions ? real.(new_pts) : Vector{Float64}[]
             return routing_points, result_intermediate, mon_result
         end
         start_parameters!(H, intermediate_rhs)
         target_parameters!(H, zeros(ComplexF64, length(rhs0)))
-        result = HomotopyContinuation.solve(
+        intermediate_solutions = solutions(result_intermediate)
+        result, _ = _solve_homotopy(
             H,
-            result_intermediate;
+            intermediate_solutions;
             catch_interrupt = catch_interrupt,
         )
         routing_points = real_solutions(result)
diff --git a/test/runtests.jl b/test/runtests.jl
index 26f8870..819c291 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -115,7 +115,7 @@ end
         expand_start_solutions_newton = false,
     )
     @test all(norm(∇r(z)-rhs0) < 1e-10 for z in S0_no_newton)
-    @test all(norm(∇r(z)) < 1e-10 for z in new_pts_no_newton)
+    @test all(norm(∇r(z)) < 1e-6 for z in new_pts_no_newton)
 
     S0, new_pts = ProjectedHypersurfaceRegions._expand_start_solutions(
         ∇r, H, S0_initial, rhs0, k;
@@ -131,6 +131,19 @@ end
     pts, res0, mon_res = critical_points(r, expand_start_solutions=false, options=options)
     @test all(norm.(∇r_symbolic.(solutions(res0))) .< 1e-12)
 
+    empty_start_solutions = Vector{Vector{ComplexF64}}()
+    empty_routing_points, empty_result, empty_mon_res = ProjectedHypersurfaceRegions._solve_and_trace(
+        MS,
+        H,
+        empty_start_solutions,
+        rhs0,
+        Vector{Vector{ComplexF64}}();
+        expand_start_solutions = false,
+    )
+    @test isempty(empty_routing_points)
+    @test !isnothing(empty_result)
+    @test isempty(solutions(empty_result))
+
     pl = generate_plot(
         r,
         [[0.0, 0.0]],

From 1a98ab619c76a1646bac55773e9d46da572128dd Mon Sep 17 00:00:00 2001
From: John Cobb <jdcobb3@gmail.com>
Date: Thu, 7 May 2026 13:34:01 -0400
Subject: [PATCH 3/3] Do not return intermediate trace as critical-point result

---
 src/critical_points.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/critical_points.jl b/src/critical_points.jl
index 7c4b521..44cce3d 100644
--- a/src/critical_points.jl
+++ b/src/critical_points.jl
@@ -360,7 +360,7 @@ function _solve_and_trace(
         )
         if interrupted
             routing_points = expand_start_solutions ? real.(new_pts) : Vector{Float64}[]
-            return routing_points, result_intermediate, mon_result
+            return routing_points, nothing, mon_result
         end
         start_parameters!(H, intermediate_rhs)
         target_parameters!(H, zeros(ComplexF64, length(rhs0)))