From ffa797424814f70e079a9782b07924cb24583542 Mon Sep 17 00:00:00 2001 From: itsdfish Date: Fri, 8 Jul 2022 07:46:33 -0400 Subject: [PATCH] update sampler and increment version --- Project.toml | 2 +- src/ParameterSpacePartitions.jl | 4 +- src/intersection_test.jl | 251 -------------------------------- test/intersection_tests.jl | 115 --------------- test/runtests.jl | 3 +- 5 files changed, 3 insertions(+), 372 deletions(-) delete mode 100644 src/intersection_test.jl delete mode 100644 test/intersection_tests.jl diff --git a/Project.toml b/Project.toml index 95fe3fc..0bb2a2f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ParameterSpacePartitions" uuid = "64931d76-8770-46b0-8423-0a5d3a7d2d72" authors = ["itsdfish"] -version = "0.3.10" +version = "0.4.0" [deps] ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" diff --git a/src/ParameterSpacePartitions.jl b/src/ParameterSpacePartitions.jl index ae4ad53..0a96d86 100644 --- a/src/ParameterSpacePartitions.jl +++ b/src/ParameterSpacePartitions.jl @@ -7,8 +7,7 @@ module ParameterSpacePartitions no_adaption!, estimate_volume, bias_correction, - psp_slices, - intersects + psp_slices export Options, Results @@ -17,7 +16,6 @@ module ParameterSpacePartitions include("structs.jl") include("sampler.jl") include("volume.jl") - include("intersection_test.jl") include("TestModels.jl") end diff --git a/src/intersection_test.jl b/src/intersection_test.jl deleted file mode 100644 index 22fbd57..0000000 --- a/src/intersection_test.jl +++ /dev/null @@ -1,251 +0,0 @@ -""" - intersects(μ1, μ2, cov1, cov2, c=2) - -Tests whether two hyperellipsoids intersect. The test returns true if -the hyperellipsoids intersection and false otherwise. - -# Arguments - -- `μ1`: centroid of ellipsoid 1 -- `μ2`: centroid of ellipsoid 2 -- `cov1`: covariance matrix of ellipsoid 1 -- `cov2`: covariance matrix of ellipsoid 2 -- `c=2`: ellipse sclar -""" -function intersects(μ1, μ2, cov1, cov2, c=2) - cov1 .*= c^2 - cov2 .*= c^2 - cho1 = cholesky(cov1) - inv1 = inv(cho1.U) - - _Q2b = inv1' * cov2 * inv1 - Q2b = Symmetric(_Q2b) - - if !isposdef(Q2b) - println("cov1") - cov1 ./= c^2 - println(cov1) - println("cov2") - cov2 ./= c^2 - println(cov2) - end - - if !issymmetric(_Q2b) - ϵ = sum(abs.(_Q2b .- _Q2b')) / length(_Q2b) - if ϵ > 1e-10 - @warn "Q2b is not symmetric: ϵ = $(ϵ). Q2b = $(Q2b)" - end - end - c2b = inv1' * (μ2 - μ1) - c2c = (c2b' * inv(cholesky(Q2b).U))' - v2c = -c2c / sqrt(c2c' * c2c) - - test_point = (v2c' * cholesky(Q2b).U)' .+ c2b - if test_point' * test_point < 1 - return true - elseif all(sign.(test_point) ≠ sign.(c2b)) - return true - else - return false - end -end - -""" - intersects(chain1, chain2, c=2) - -Tests whether two hyperellipsoids intersect. The test returns true if -the hyperellipsoids intersection and false otherwise. - -# Arguments - -- `chain1`: chain object -- `chain2`: chain object -- `c=2`: ellipse sclar -""" -function intersects(chain1, chain2, c=2) - mat1 = to_matrix(chain1) - mat2 = to_matrix(chain2) - μ1 = mean(mat1, dims=1)[:] - μ2 = mean(mat2, dims=1)[:] - cov1 = cov(mat1) - add_variance!(cov1) - cov2 = cov(mat2) - add_variance!(cov2) - if !isposdef(cov1) - println("cov1") - println(cov1) - end - if !isposdef(cov2) - println("cov2") - println(cov2) - end - return intersects(μ1, μ2, cov1, cov2) -end - -function add_variance!(x) - if any(x -> isapprox(x, 0; atol = 1e-10), diag(x)) - x[diagind(x)] .= eps() - end - return nothing -end - -to_matrix(x) = reduce(vcat, transpose.(x.all_parms)) - -""" - get_group_indices(chains, chain_indices) - -Sorts chains of the same pattern into non-overlapping groups. The vector -[[1,2],[3,4]] indices chains 1 and 2 are located in region R₁ and chains -3 and 4 are located in region R₂. - -# Arguments - -- `chains`: a vector of chains -- `group`: a vector of indices corresponding to chains with the same pattern -""" -function get_group_indices(chains, chain_indices) - # group chain indices according to region - indices = Vector{Vector{Int}}() - # first index group will have c = 1 - push!(indices, [chain_indices[1]]) - n_groups = length(indices) - n_chains = length(chain_indices) - # group index - g = 1 - # loop through each chain - for i ∈ 2:n_chains - # loop through each index group - c = chain_indices[i] - while g ≤ n_groups - # if chain c matches region index in group g, - # push c into index group g - if intersects(chains[indices[g][1]], chains[c]) - push!(indices[g], c) - merge_chains!(chains[indices[g][1]], chains[c]) - break - end - g += 1 - end - # add new index group - if g > n_groups - # add new group - push!(indices, [c]) - # increment number of groups - n_groups += 1 - end - # reset group counter - g = 1 - end - return indices -end - -""" - remove_redundant_chains!(chains, indices) - -Removes chains that have the same pattern and location in the parameter space. -For example, in the vector `indices` [[1,3],[34,50]], chains indexed in positions 1 and 3 have the same -pattern and location, as do chains indexed at positions 34 and 50. Only the first element of each sub-vector -is retained (i.e., [1,50]) - -# Arguments - -- `chains`: a vector containing all chain objects -- `indices`: a nested vector that maps chains of the same pattern and location to `chains` -""" -function remove_redundant_chains!(chains, indices) - k_indices = Int[] - for i in indices - push!(k_indices, i[1]) - end - r_indices = setdiff(vcat(indices...), k_indices) - sort!(r_indices) - deleteat!(chains, r_indices) - return nothing -end - -""" - group_by_pattern(chains) - -Groups chains according to pattern and returns a nested vector of chain indices. The vector -[[1,2],[3,4]] indices chains 1 and 2 are located in region R₁ and chains -3 and 4 are located in region R₂. - -# Arguments - -- `chains`: a vector of all chains -""" -function group_by_pattern(chains) - patterns = map(c -> c.pattern, chains) - u_patterns = unique(patterns) - n_patterns = length(u_patterns) - chain_indices = [Vector{Int}() for _ in 1:n_patterns] - for c in 1:length(chains) - g = findfirst(p -> chains[c].pattern == p, u_patterns) - push!(chain_indices[g], c) - end - return chain_indices -end - -function remove_nonposdef!(chains) - n_chains = length(chains) - test_vals = fill(false, n_chains) - for c in 1:n_chains - mat = to_matrix(chains[c]) - covar = cov(mat) - test_vals[c] = !isposdef(covar) - end - deleteat!(chains, test_vals) - return nothing -end - -""" - make_unique!(chains, options; timer=nothing, show_time=false) - -This function sorts chains by pattern and merges chains that are in the same region - -# Arguments - -- `chains`: a vector of all chains -- `options`: an Options configuration object for the PSP algorithm - -# Keywords - -- `timer`: monitors time -- `show_time=false`: shows time if true -""" -function make_unique!(chains, options; timer=nothing, show_timer=false) - for iter in 1:2 - _make_unique(chains, options; iter, timer, show_timer) - end - return nothing -end - -function _make_unique(chains, options; iter=0, timer=nothing, show_timer = false) - remove_nonposdef!(chains) - chain_indices = group_by_pattern(chains) - all_indices = Vector{Vector{Int}}() - for c in chain_indices - temp = get_group_indices(chains, c) - push!(all_indices, temp...) - show_timer ? next!(timer; showvalues=show_values(iter)) : nothing - end - remove_redundant_chains!(chains, all_indices) - return nothing -end - -""" - merge_chains!(chain1, chain2) - -Merges chain2 into chain1 on fields `all_parms`, `acceptance`, and `radii`. - -# Arguments - -- `chain1`: a chain object -- `chain2`: a chain object -""" -function merge_chains!(chain1, chain2) - push!(chain1.all_parms, chain2.all_parms...) - push!(chain1.acceptance, chain2.acceptance...) - push!(chain1.radii, chain2.radii...) - return nothing -end \ No newline at end of file diff --git a/test/intersection_tests.jl b/test/intersection_tests.jl deleted file mode 100644 index f5abbee..0000000 --- a/test/intersection_tests.jl +++ /dev/null @@ -1,115 +0,0 @@ -@testset verbose = true "Intersection Test" begin - @safetestset "Non-intersecting hyperellipsoids" begin - # Let d be the distance between centroids of ellispoids e1 and e2. The test should always return false - # if d > m_1 / 2 + m_2 /2, where m_1 is the major axis of e1 and m_2 is the major axis of e2. - # returns axes of ellipsoid with scaling factor of 2 - using ParameterSpacePartitions, Test, LinearAlgebra - using Random - include("intersection_utilities.jl") - Random.seed!(8744) - - for _ in 1:100 - n_dims = rand(2:10) - μ₁ = fill(0, n_dims) - Σ₁ = rand_cov_mat(n_dims) - Σ₂ = rand_cov_mat(n_dims) - r₁ = get_axes(Σ₁) / 2 - r₂ = get_axes(Σ₂) / 2 - r_max = maximum(r₁) + maximum(r₂) - x = randn(n_dims) - μ₂ = (x / sqrt(x' * x)) * r_max * 1.01 - @test intersects(μ₁, μ₂, Σ₁, Σ₂) == false - end - end - - @safetestset "Intersecting hyperellipsoids" begin - using ParameterSpacePartitions, Test, LinearAlgebra - using Random - include("intersection_utilities.jl") - Random.seed!(354) - - for _ in 1:100 - n_dims = rand(2:10) - μ₁ = fill(0, n_dims) - Σ₁ = rand_cov_mat(n_dims) - Σ₂ = rand_cov_mat(n_dims) - r₁ = get_axes(Σ₁) / 2 - r₂ = get_axes(Σ₂) / 2 - r_min = min(minimum(r₁), minimum(r₂)) - x = randn(n_dims) - # should always intersect if less than twice minimum axis - μ₂ = (x / sqrt(x' * x)) * 2 * r_min * .99 - @test intersects(μ₁, μ₂, Σ₁, Σ₂) - end - end - - @safetestset "hyperspheres" begin - using ParameterSpacePartitions, Test, LinearAlgebra - using Random - include("intersection_utilities.jl") - Random.seed!(6521) - - for _ in 1:100 - n_dims = rand(2:10) - μ₁ = fill(0, n_dims) - Σ₁ = Array(I(n_dims) * 1.0) - Σ₂ = Array(I(n_dims) * 1.0) - r = 2 * 2 - x = randn(n_dims) - μ₂ = (x / sqrt(x' * x)) * r * .99 - @test intersects(μ₁, μ₂, Σ₁, Σ₂) - end - - for _ in 1:100 - n_dims = rand(2:10) - μ₁ = fill(0, n_dims) - Σ₁ = Array(I(n_dims) * 1.0) - Σ₂ = Array(I(n_dims) * 1.0) - # radius of hypersphere - r = 2 - x = randn(n_dims) - # distance is twice the radius plus 1% - μ₂ = (x / sqrt(x' * x)) * 2 * r * 1.01 - @test intersects(μ₁, μ₂, Σ₁, Σ₂) == false - end - end - - @safetestset "remove redundant chains" begin - using ParameterSpacePartitions, Test, LinearAlgebra - using ParameterSpacePartitions: Chain - using ParameterSpacePartitions: intersects, remove_redundant_chains! - using ParameterSpacePartitions: get_group_indices, group_by_pattern - using ParameterSpacePartitions: merge_chains! - using Statistics - - chains = [ - Chain(1, [.1,.3], 1, 1), - Chain(2, [.1,.3,], 1, 1), - Chain(3, [.1,.3,], 1, 1), - Chain(4, [.1,.3,], 2, 1), - Chain(5, [.1,.3,], 2, 1), - - ] - points = [rand(2) for _ in 1:100] - push!(chains[1].all_parms, points...) - push!(chains[2].all_parms, points...) - points = [rand(2) .+ 2 for _ in 1:100] - push!(chains[3].all_parms, points...) - push!(chains[4].all_parms, points...) - push!(chains[5].all_parms, points...) - - c_indices = group_by_pattern(chains) - @test c_indices[1] == [1,2,3] - @test c_indices[2] == [4,5] - - indices = map(c -> get_group_indices(chains, c), c_indices) - indices = vcat(indices...) - remove_redundant_chains!(chains, indices) - - @test length(chains) == 3 - @test chains[1].chain_id == 1 - @test chains[2].chain_id == 3 - @test chains[3].chain_id == 4 - end -end - diff --git a/test/runtests.jl b/test/runtests.jl index ad1a668..1f0ba78 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,4 @@ using Test, SafeTestsets include("volume_tests.jl") -include("sampler_tests.jl") -include("intersection_tests.jl") \ No newline at end of file +include("sampler_tests.jl") \ No newline at end of file