From bd59ba10bb4131cefcaca1ce8870d90acc2c7ac6 Mon Sep 17 00:00:00 2001
From: JoelTrent <79883375+JoelTrent@users.noreply.github.com>
Date: Wed, 28 Jun 2023 18:06:26 +1200
Subject: [PATCH] Improve performance of simultaneous method by allowing one of
 the internal or external points vectors to contain non-unique points for
 point pairs (using random numbers more efficiently)

---
 .../bivariate_profile_likelihood.jl           |  9 +++--
 src/bivariate_methods/vectorsearch.jl         | 35 ++++++++++++++++---
 src/types/bivariate_methods.jl                | 31 ++++++++++++----
 3 files changed, 59 insertions(+), 16 deletions(-)

diff --git a/src/bivariate_methods/bivariate_profile_likelihood.jl b/src/bivariate_methods/bivariate_profile_likelihood.jl
index b7d0d1f..953d38e 100644
--- a/src/bivariate_methods/bivariate_profile_likelihood.jl
+++ b/src/bivariate_methods/bivariate_profile_likelihood.jl
@@ -177,7 +177,8 @@ function bivariate_confidenceprofile(bivariate_optimiser::Function,
             boundary, internal = bivariate_confidenceprofile_vectorsearch(
                                     bivariate_optimiser, model, 
                                     num_points, consistent, ind1, ind2,
-                                    mle_targetll, save_internal_points, channel)
+                                    mle_targetll, save_internal_points, channel,
+                                    min_proportion_unique=method.min_proportion_unique)
 
         elseif method isa RadialRandomMethod
             boundary, internal = bivariate_confidenceprofile_vectorsearch(
@@ -370,10 +371,8 @@ function bivariate_confidenceprofiles!(model::LikelihoodModel,
     len_θcombinations = length(θcombinations)
     len_θcombinations > 0 || return nothing
 
-
-
-    tasks_per_profile = get_bivariate_method_tasknumbers(method, num_points)
-    totaltasks = tasks_per_profile * len_θcombinations
+    tasks_per_profile = get_bivariate_method_tasknumbers.(Ref(method), num_new_points)
+    totaltasks = sum(tasks_per_profile)
     # channel_buffer_size = min(ceil(Int, tasks_per_profile * 0.05), 30)
     channel_buffer_size = 2
     channel = RemoteChannel(() -> Channel{Bool}(channel_buffer_size))
diff --git a/src/bivariate_methods/vectorsearch.jl b/src/bivariate_methods/vectorsearch.jl
index 56718ec..8ac8b22 100644
--- a/src/bivariate_methods/vectorsearch.jl
+++ b/src/bivariate_methods/vectorsearch.jl
@@ -16,13 +16,16 @@ function findNpointpairs_simultaneous!(p::NamedTuple,
                                         ind2::Int,
                                         mle_targetll::Float64,
                                         save_internal_points::Bool,
-                                        biv_opt_is_ellipse_analytical::Bool)
+                                        biv_opt_is_ellipse_analytical::Bool,
+                                        min_proportion_unique::Real)
 
     internal  = zeros(2,num_points)
     internal_all = zeros(model.core.num_pars, save_internal_points ? num_points : 0)
     ll_values = zeros(save_internal_points ? num_points : 0)
     external = zeros(2,num_points)
 
+    min_num_unique = ceil(Int, min_proportion_unique*num_points)
+
     Ninside=0; Noutside=0
     iters=0
     while Noutside<num_points && Ninside<num_points
@@ -49,7 +52,7 @@ function findNpointpairs_simultaneous!(p::NamedTuple,
     end
 
     # while Ninside < N && iters < maxIters
-    while Ninside < num_points
+    while Ninside < min_num_unique
         x, y = generatepoint(model, ind1, ind2)
         p.pointa .= [x,y]
         g = bivariate_optimiser(0.0, p)
@@ -69,7 +72,7 @@ function findNpointpairs_simultaneous!(p::NamedTuple,
     end
 
     # while Noutside < N && iters < maxIters
-    while Noutside < num_points
+    while Noutside < min_num_unique
         x, y = generatepoint(model, ind1, ind2)
         p.pointa .= [x,y]
         if bivariate_optimiser(0.0, p) < 0
@@ -79,6 +82,28 @@ function findNpointpairs_simultaneous!(p::NamedTuple,
         iters+=1
     end
 
+    if Ninside < num_points
+        num_unique = Ninside*1
+        while Ninside < num_points
+            i, j = Ninside+1, min(num_points-Ninside, num_unique)
+            internal[:, i:(i+j-1)] .= internal[:,1:j]
+            Ninside += num_unique
+        end
+        
+        if save_internal_points
+            ll_values = ll_values[1:num_unique]
+            internal_all = internal_all[:, 1:num_unique]
+        end
+
+    elseif Noutside < num_points
+        num_unique = Ninside * 1
+        while Noutside < num_points
+            i, j = Noutside+1, min(num_points-Noutside, num_unique)
+            external[:, i:(i+j-1)] .= external[:, 1:j]
+            Noutside += num_unique
+        end
+    end
+
     if save_internal_points && biv_opt_is_ellipse_analytical
         get_ωs_bivariate_ellipse_analytical!(internal_all, num_points,
                                                     p.consistent, ind1, ind2, 
@@ -262,6 +287,7 @@ function bivariate_confidenceprofile_vectorsearch(bivariate_optimiser::Function,
                                                     save_internal_points::Bool,
                                                     channel::RemoteChannel;
                                                     num_radial_directions::Int=0,
+                                                    min_proportion_unique::Real=1.0,
                                                     ellipse_confidence_level::Float64=-1.0,
                                                     ellipse_start_point_shift::Float64=0.0,
                                                     ellipse_sqrt_distortion::Float64=0.0)
@@ -291,7 +317,8 @@ function bivariate_confidenceprofile_vectorsearch(bivariate_optimiser::Function,
         if num_radial_directions == 0
             internal, internal_all, ll_values, external = 
                 findNpointpairs_simultaneous!(p, bivariate_optimiser, model, num_points, ind1, ind2,
-                                                mle_targetll, save_internal_points, biv_opt_is_ellipse_analytical)
+                                                mle_targetll, save_internal_points, biv_opt_is_ellipse_analytical,
+                                                min_proportion_unique)
         else
             internal, internal_all, ll_values, external = 
                 findNpointpairs_radialrandom!(p, bivariate_optimiser, model, num_points, 
diff --git a/src/types/bivariate_methods.jl b/src/types/bivariate_methods.jl
index 721da76..a6af33b 100644
--- a/src/types/bivariate_methods.jl
+++ b/src/types/bivariate_methods.jl
@@ -233,15 +233,23 @@ struct RadialRandomMethod <: AbstractBivariateVectorMethod
 end
 
 """
-    SimultaneousMethod()
+    SimultaneousMethod(min_proportion_unique::Real=0.5)
 
-Method for finding the bivariate boundary of a confidence profile by finding internal and external boundary points using a uniform random distribution on provided bounds, pairing these points in the order they're found and bracketing between each pair (see [`PlaceholderLikelihood.findNpointpairs_simultaneous!`](@ref) and [`PlaceholderLikelihood.bivariate_confidenceprofile_vectorsearch`](@ref)).
+Method for finding the bivariate boundary of a confidence profile by finding internal and external boundary points using a uniform random distribution on provided bounds, pairing these points in the order they're found and bracketing between each pair (see [`PlaceholderLikelihood.findNpointpairs_simultaneous!`](@ref) and [`PlaceholderLikelihood.bivariate_confidenceprofile_vectorsearch`](@ref)). Values of `min_proportion_unique` lower than zero may improve performance if finding either internal points or external points is difficult given the specified bounds on interest parameters.
+    
+
+# Arguments
+- `min_proportion_unique`: a proportion ∈ (0.0, 1.0] for the minimum allowed proportion of unique points in one of the internal or external point vectors. One of these vectors will be fully unique. Whichever vector is not unique, will have at least `min_proportion_unique` unique points. Default is 0.5
 
 # Details
 
 Recommended for use with the [`LogLikelihood`](@ref) profile type. 
 
-The method uniformly samples the region specified by the bounds for the two interest parameters until as many internal and external boundary points as the desired number of boundary points are found. These points are paired in the order they are found. A bracketing method is then used to find a boundary point between the point pair (the external point and the internal point).
+When `min_proportion_unique = 1.0` the method uniformly samples the region specified by the bounds for the two interest parameters until as many internal and external boundary points as the desired number of boundary points are found. These points are paired in the order they are found. 
+
+If `min_proportion_unique < 1.0` the method uniformly samples the region specified by the bounds for the two interest parameters until either the internal or external number of boundary points found is as many as the number of desired boundary points. For whichever vector is less than the number of desired boundary points, we keep searching for that kind of point until at least `ceil(min_proportion_unique * num_points)` points have been found. Once `ceil(min_proportion_unique * num_points)` points have been found these points are used to fill the remainder of the vector, in order found. Internal and external points in each vector are then paired in the order they are found.
+
+A bracketing method is then used to find a boundary point between each point pair (the external point and the internal point).
 
 [`RadialRandomMethod`](@ref) and [`IterativeBoundaryMethod`](@ref) are preferred over this method from a computational performance standpoint as they require fewer log-likelihood evalutions (when [`RadialRandomMethod`](@ref) has parameter `num_radial_directions` > 1). 
 
@@ -253,18 +261,27 @@ This method can find multiple boundaries (if they exist).
 
 ## Impact of parameter bounds
 
-If a parameter bound is in the way of reaching the boundary, points will not be put on that bound. Additionally, if the true boundary is very close to a parameter bound, the method will struggle to find this region of the boundary. This is because finding the boundary in this location requires generating a random point between the boundary and the parameter bound, which becomes more difficult the closer they are. Interest parameter bounds that have ranges magnitudes larger than the range of the boundary will make finding internal points very difficult, requiring a lot of computational effort. Similarly, the inverse will be true if external points are hard to find. The method will fail if the interest parameter bounds are fully contained by the boundary.
+If a parameter bound is in the way of reaching the boundary, points will not be put on that bound. Additionally, if the true boundary is very close to a parameter bound, the method will struggle to find this region of the boundary. This is because finding the boundary in this location requires generating a random point between the boundary and the parameter bound, which becomes more difficult the closer they are. 
+
+Interest parameter bounds that have ranges magnitudes larger than the range of the boundary will make finding internal points difficult, requiring a lot of computational effort. Similarly, the inverse will be true if external points are hard to find. Smaller values of `min_proportion_unique` will improve performance in these cases. 
+
+The method will fail if the interest parameter bounds are fully contained by the boundary.
 
 # Internal Points
 
-Finds `num_points` internal points.
+Finds at least `ceil(min_proportion_unique * num_points)` internal points.
 
 # Supertype Hiearachy
 
 `SimultaneousMethod <: AbstractBivariateVectorMethod <: AbstractBivariateMethod <: Any`
 """
-struct SimultaneousMethod <: AbstractBivariateVectorMethod end
-
+struct SimultaneousMethod <: AbstractBivariateVectorMethod 
+    min_proportion_unique::Real
+    function SimultaneousMethod(min_proportion_unique::Real=0.5)
+        (0.0 < min_proportion_unique && min_proportion_unique <= 1.0) || throw(DomainError("min_proportion_unique must be in the interval (0.0,1.0]"))
+        return new(min_proportion_unique)
+    end
+end
 
 """
     Fix1AxisMethod()