Draw the warm-up start parameters on the scale of the targets

src/methods.jl

@@ -122,7 +122,8 @@ at `index` perturbed by `perturbation_size`. The warmup system is used to perfor
 using all other systems in the parameter sweep. It is very efficient for systems with minimal
 bifurcation in the parameter sweep. The Warm Up method should in theory guarantee to find all
 solutions, however, if the `start_parameters` is not proper (to close to the real line)
-it could miss some solutions.
+it could miss some solutions. Without explicit `start_parameters`, each one is drawn as a
+random complex number on the scale of the largest magnitude that parameter takes in the sweep.
 
 See[HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org/guides/many-systems/)
 for more information.

src/solve_homotopy.jl

@@ -189,17 +189,32 @@ function _prepare_input_params(
     return unique_fixed, tuple_to_vector.(input_array)
 end
 
+"""
+Random complex start parameters for the warm-up, on the scale of the targets: a standard
+complex normal draw per parameter times the largest magnitude it takes in
+`parameter_values`. The parameter homotopy to every target is then a relative change of
+order one. Draws of order one towards targets of order ten or more make the path tracking
+ill-conditioned and lose solutions.
+"""
+function random_start_parameters(parameter_values::AbstractArray{<:AbstractVector})
+    example = first(parameter_values)
+    scale = map(eachindex(example)) do i
+        magnitude = maximum(p -> abs(p[i]), parameter_values)
+        return iszero(magnitude) ? one(magnitude) : magnitude
+    end
+    return randn(complex(eltype(example)), length(example)) .* scale
+end
+
 "Uses HomotopyContinuation to solve `problem` at specified `parameter_values`."
 function _get_raw_solution(
     problem::HomotopyContinuationProblem, method::WarmUp, parameter_values; show_progress
 )
     warm_up_method = method.warm_up_method
 
-    example_p = parameter_values[1]
-    if isempty(method.start_parameters)
-        start_parameters = randn(complex(eltype(example_p)), length(example_p))
+    start_parameters = if isempty(method.start_parameters)
+        random_start_parameters(parameter_values)
     else
-        start_parameters = method.start_parameters
+        method.start_parameters
     end
 
     warmup_solution = HC.solve(

test/steady_states/jacobian.jl

@@ -108,6 +108,14 @@ end
 
     elapsed = @elapsed get_steady_states(harmonic_eq, swept, fixed; show_progress=false)
     @test elapsed < 10
+
+    # Issue #170: the random start point of the warm-up must not decide how many solutions
+    # are found. Seed 5 found 4 of the 9 when the start parameters ignored the target scale.
+    @testset "warm-up seed $seed" for seed in UInt32(1):UInt32(12)
+        method = WarmUp(; seed)
+        result = get_steady_states(harmonic_eq, method, swept, fixed; show_progress=false)
+        @test length(result.solutions[1]) == 9
+    end
 end
 
 @testset "autodiff backend gives the same Jacobian" begin