fix: per-interpolator error handling to prevent HC crashes from losing results

Wrap per-interpolator body in try/catch so an InexactError(NaN) in one
interpolator (e.g. S3SE) doesn't crash the entire estimation, preserving
results from other interpolators. Also add per-point try/catch in
solve_with_hc_parameterized so a single bad shooting point doesn't kill
the whole interpolator pass. Includes resolve_states try/catch fallback
and fhdn(5) fix.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: Oren Bassik

src/core/homotopy_continuation.jl

@@ -917,6 +917,7 @@ function solve_with_hc_parameterized(poly_system, solve_vars, data_vars, param_v
 	initial_solution_count = 0
 
 	for (i, current_params) in enumerate(param_values_list)
+		try
 		if i == 1 || isnothing(prev_all_solutions) || isempty(prev_all_solutions)
 			# Fresh solve at first point - get ALL solutions
 			if debug
@@ -989,6 +990,20 @@ function solve_with_hc_parameterized(poly_system, solve_vars, data_vars, param_v
 		# Track ALL solutions to next point (real + complex)
 		prev_all_solutions = all_solutions
 		prev_params = current_params
+
+		catch e
+			@error "[HC-PARAM] Point $i failed" exception=(e, catch_backtrace())
+			println(stderr, "[HC-CRASH] Point $i threw $(typeof(e)): $e")
+			println(stderr, "[HC-CRASH] current_params = $current_params")
+			if any(isnan, current_params) || any(isinf, current_params)
+				println(stderr, "[HC-CRASH] WARNING: current_params contains NaN/Inf!")
+			end
+			# Push empty results for this point
+			push!(all_real_results, Vector{Vector{Float64}}())
+			# Reset tracking so next point does a fresh solve
+			prev_all_solutions = nothing
+			prev_params = nothing
+		end
 	end
 
 	return all_real_results

src/core/optimized_multishot_estimation.jl

@@ -1255,6 +1255,8 @@ function optimized_multishot_parameter_estimation(PEP::ParameterEstimationProble
 			interp_solutions = []
 			interp_time_indices = Int[]
 
+		try  # Catch errors in individual interpolators so one crash doesn't lose all results
+
 		if use_param_homotopy
 			# ============================================================================
 			# PARAMETER HOMOTOPY PATH
@@ -1541,6 +1543,20 @@ function optimized_multishot_parameter_estimation(PEP::ParameterEstimationProble
 			end
 		end  # end if use_param_homotopy
 
+		catch e
+			@error "Interpolator $interp_sym failed" exception=(e, catch_backtrace())
+			println(stderr, "[HC-CRASH] Interpolator $interp_sym threw $(typeof(e)): $e")
+			if @isdefined(param_values_list) && isa(param_values_list, AbstractVector)
+				println(stderr, "[HC-CRASH] param_values_list had $(length(param_values_list)) points")
+				for (pi, pv) in enumerate(param_values_list)
+					if any(isnan, pv) || any(isinf, pv)
+						println(stderr, "[HC-CRASH]   Point $pi: contains NaN/Inf: $pv")
+					end
+				end
+			end
+			# interp_solutions is empty, so accumulation below appends nothing
+		end  # end try/catch for interpolator
+
 			# Accumulate this pass's solutions into the global pool
 			append!(all_solutions, interp_solutions)
 			append!(solution_time_indices, interp_time_indices)
@@ -1650,20 +1666,29 @@ function optimized_multishot_parameter_estimation(PEP::ParameterEstimationProble
 					known_param_dict = OrderedDict{Any, Float64}(k => Float64(v) for (k, v) in params_dict)
 
 					# Algebraic re-solve at t=0 with fixed parameters
-					state_solutions, state_vars = resolve_states_with_fixed_params(
-						PEP.model.system,
-						PEP.measured_quantities,
-						PEP.data_sample,
-						good_deriv_level,
-						good_udict,
-						good_varlist,
-						good_DD,
-						known_param_dict,
-						interpolants;
-						si_template = si_template,
-						time_index = 1,
-						diagnostics = opts.diagnostics,
-					)
+					state_solutions = Vector{Vector{Float64}}()
+					state_vars = Any[]
+					try
+						state_solutions, state_vars = resolve_states_with_fixed_params(
+							PEP.model.system,
+							PEP.measured_quantities,
+							PEP.data_sample,
+							good_deriv_level,
+							good_udict,
+							good_varlist,
+							good_DD,
+							known_param_dict,
+							interpolants;
+							si_template = si_template,
+							time_index = 1,
+							diagnostics = opts.diagnostics,
+						)
+					catch e
+						@warn "[RESOLVE] Algebraic re-solve failed for parameter set" exception = (e, catch_backtrace())
+						if !opts.nooutput
+							println("  Algebraic re-solve failed ($(typeof(e))), falling back to data-derived ICs")
+						end
+					end
 
 					if !isempty(state_solutions)
 						# Build ParameterEstimationResult for each algebraic state solution

src/types/estimation_options.jl

@@ -344,7 +344,7 @@ function get_interpolator_function(method::InterpolatorMethod, custom::Union{Not
 	elseif method == InterpolatorAAADOld
 		return aaad_old_reliable
 	elseif method == InterpolatorFHD
-		return fhd5  # Default to degree 5 FHD
+		return fhdn(5)  # Default to degree 5 FHD
 	elseif method == InterpolatorAGP
 		return agp_gpr
 	elseif method == InterpolatorAGPRobust