diff --git a/Project.toml b/Project.toml index b044ef2..9215dcf 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,10 @@ name = "ParameterSpacePartitions" uuid = "64931d76-8770-46b0-8423-0a5d3a7d2d72" authors = ["itsdfish"] -version = "0.3.5" +version = "0.3.7" [deps] +ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" @@ -17,6 +18,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" ThreadedIterables = "11d239b0-c0b9-11e8-1935-d5cfa53abb03" [compat] +ComponentArrays = "0.11.9" ConcreteStructs = "v0.2.3" DataFrames = "1.3.0" Distributions = "0.23.0,0.24.0,v0.25.37" diff --git a/src/ParameterSpacePartitions.jl b/src/ParameterSpacePartitions.jl index c5f0b92..b488269 100644 --- a/src/ParameterSpacePartitions.jl +++ b/src/ParameterSpacePartitions.jl @@ -1,6 +1,6 @@ module ParameterSpacePartitions using Requires, Distributions, ConcreteStructs, LinearAlgebra - using ThreadedIterables, SpecialFunctions + using ThreadedIterables, SpecialFunctions, ComponentArrays export find_partitions, adapt!, diff --git a/src/volume.jl b/src/volume.jl index f6a04e4..744bf1a 100644 --- a/src/volume.jl +++ b/src/volume.jl @@ -7,19 +7,26 @@ volume_hypersphere(n, r=1) = π^(n / 2) / gamma(1 + n / 2) * r^n Computes the log volume of a hypersphere # Arguments + - `n`: the number of dimensions + +# Keywords + - `r=1`: the radius """ log_vol_hypersphere(n; r=1) = (n / 2) * log(π) - loggamma(1 + n / 2) + n * log(r) """ - log_vol_ellipsiod(n, r=1) + log_vol_ellipsiod(n, cov_mat; r = 1) -Computes the log volume of a hypersphere +Computes the log volume of a ellipsoid # Arguments - `n`: the number of dimensions - `cov_mat`: covariance matrix + +# Keywords + - `r=1`: the radius """ function log_vol_ellipsiod(n, cov_mat; r = 1) @@ -29,13 +36,16 @@ function log_vol_ellipsiod(n, cov_mat; r = 1) end """ - volume_ellipsoid(n, cov_mat, r=1) + volume_ellipsoid(n, cov_mat; r=1) Computes the volume of an ellipsoid # Arguments - `n`: the number of dimensions - `cov_mat`: covariance matrix + +# Keywords + - `r=1`: the radius """ function volume_ellipsoid(n, cov_mat; r = 1) @@ -63,7 +73,7 @@ function sample_ellipsoid(μ, n_dims, cov_mat) end """ - sample_ellipsoid_surface(μ, n_dims, cov_mat) + sample_ellipsoid_surface(μ, n_dims, cov_mat; r=1) Samples a point uniformly from the surface of an ellipsoid @@ -72,18 +82,44 @@ Samples a point uniformly from the surface of an ellipsoid - `μ`: mean of ellipsoid on each dimension - `n_dims`: the number of dimensions - `cov_mat`: covariance matrix + +# Keywords + - `r=1`: radius of underlying unit hypersphere """ -function sample_ellipsoid_surface(μ, n_dims, cov_mat, r=1) +function sample_ellipsoid_surface(μ, n_dims, cov_mat; r=1) x = random_position(r, n_dims) return μ .+ sqrt((n_dims + 2) * cov_mat) * x end """ - bias_correction(model, p_fun, points, target, bounds, n_sim, args...; kwargs...) + bias_correction( + model, + p_fun, + points, + target, + bounds, + args...; + parm_names = Symbol.("p", 1:size(points,2)), + n_sim=10_000, + kwargs... + ) Uses the hit or miss technique to adjust volume estimates based on ellipsoids. +# Arguments + +- `model`: a model function that returns predictions given a vector of parameters +- `p_fun`: a function that that evaluates the qualitative data pattern +- `points`: a p x n matrix of sampled points where p is the number of parameters and n is the number of samples +- `bounds`: a vector of tuples representing (lowerbound, upperbound) for each dimension in +- `args...`: arguments passed to `model` and `p_fun` + +# Keywords + +- `parm_names`: parameter names with default `p1`, `p2,`... `pn` +- `n_sim=10_000`: number of simulations for hit or miss adjustment +- `kwargs...`: optional keyword arguments for `model` and `p_fun` """ function bias_correction( model, @@ -92,11 +128,13 @@ function bias_correction( target, bounds, args...; + parm_names = Symbol.("p", 1:size(points,2)), n_sim=10_000, kwargs... ) - μ = mean(points, dims=1)[:] + _μ = mean(points, dims=1)[:] + μ = ComponentArray(_μ, make_axis(parm_names)) cov_mat = cov(points) n = length(μ) @@ -155,6 +193,7 @@ function estimate_volume( bounds, args...; n_sim = 10_000, + parm_names = Symbol.("p", 1:size(points,2)), kwargs... ) @@ -166,6 +205,7 @@ function estimate_volume( bounds, args...; n_sim, + parm_names, kwargs... ) @@ -173,4 +213,6 @@ function estimate_volume( volume = volume_ellipsoid(cov_mat) volume *= cf return volume -end \ No newline at end of file +end + +make_axis(symbols) = Axis(NamedTuple(symbols .=> eachindex(symbols))) \ No newline at end of file diff --git a/src/volume_df.jl b/src/volume_df.jl index a4141ec..f2d08e8 100644 --- a/src/volume_df.jl +++ b/src/volume_df.jl @@ -1,5 +1,34 @@ using DataFrames +""" + estimate_volume( + model, + p_fun, + df::SubDataFrame, + bounds, + args...; + n_sim = 10_000, + parm_names, + kwargs... + ) + + +Estimate volume of region with an eillipsoid and hit or miss bias adjustment. + +# Arguments +- `model`: a model function that returns predictions given a vector of parameters +- `p_fun`: a function that that evaluates the qualitative data pattern +- `df::SubDataFrame`: a SubDataFrame containing MCMC samples in columns identities by `parm_names` + and column containing data patterns named `pattern` +- `bounds`: a vector of tuples representing (lowerbound, upperbound) for each dimension in +- `args...`: arguments passed to `model` and `p_fun` + +# Keywords + +- `n_sim=10_000`: the number of samples for hit or miss bias adjustment +- `parm_names`: a vector of symbols representing parameter columns in `df` +- `kwargs...`: additional keyword arguments passed to `model` or `p_fun` +""" function estimate_volume( model, p_fun, @@ -19,8 +48,9 @@ function estimate_volume( p_fun, points, target, - bounds, + bounds, args...; + parm_names, n_sim, kwargs... ) diff --git a/temp/temp copy 2.jl b/temp/temp copy 2.jl deleted file mode 100644 index ff0c6da..0000000 --- a/temp/temp copy 2.jl +++ /dev/null @@ -1,54 +0,0 @@ -cd(@__DIR__) -using Pkg -Pkg.activate("..") -using Revise, Distributions, ParameterSpacePartitions -using ParameterSpacePartitions.TestModels -using LinearAlgebra, Random, DataFrames - -#Random.seed!(54545) -# dimensions of the hypbercue -n_dims = 7 -# number of partitions -n_obj = 50 -# distribution of radii -r_dist = () -> rand(Uniform(.15, .25)) -# number of starting points -n_start = 10 - -# partition boundaries -hyperspheres = set_locations(r_dist, n_dims, n_obj) -bounds = fill((0, 1), n_dims) - -sample(bounds) = map(b -> rand(Uniform(b...)), bounds) - -init_parms = map(_ -> sample(bounds), 1:n_start) - -options = Options(; - radius = .25, - bounds, - n_iters = 5000, - parallel = false, - init_parms, - λ = .2, - t_rate = .40, - adapt_radius! = no_adaption!, - trace_on = false -) - -results = find_partitions( - model, - p_fun, - options, - hyperspheres; -) - -df = DataFrame(results) - -parm_names = Symbol.("p",1:n_dims) -transform!( - df, - :parms => identity => parm_names -) - -groups = groupby(df, :chain_id) -mean_accept = combine(groups, :acceptance=>mean=>:mean) diff --git a/temp/temp copy.jl b/temp/temp copy.jl deleted file mode 100644 index d9be139..0000000 --- a/temp/temp copy.jl +++ /dev/null @@ -1,56 +0,0 @@ -cd(@__DIR__) -using Pkg -Pkg.activate("..") -using ParameterSpacePartitions -using ParameterSpacePartitions.TestModels -using Revise, Random, Distributions -using DataFrames, LinearAlgebra -#Random.seed!(2001) - -# dimensions of the hypbercue -n_dims = 2 -# number of partitions -n_part = 50 -# number of starting points -n_start = 1 - - -# partition boundaries -polytopes = [Polytope(rand(n_dims)) for i in 1:n_part] -bounds = fill((0, 1), n_dims) - -sample(bounds) = map(b -> rand(Uniform(b...)), bounds) - -init_parms = map(_ -> sample(bounds), 1:n_start) - -options = Options(; - radius = .1, - bounds, - n_iters = 500, - parallel = false, - init_parms, - λ = .2 -) - -results = find_partitions( - model, - p_fun, - options, - polytopes -) - -df = DataFrame(results) -transform!( - df, - :parms => identity => [:p1, :p2] -) - -groups = groupby(df, :pattern) - -groups = groupby(df, :chain_id) -mean_accept = combine(groups, :acceptance=>mean=>:mean) - - -groups = groupby(df, :pattern) - -@df df scatter(:p1, :p2, group=:pattern, leg=false) \ No newline at end of file diff --git a/temp/temp.jl b/temp/temp.jl deleted file mode 100644 index dfc172d..0000000 --- a/temp/temp.jl +++ /dev/null @@ -1,76 +0,0 @@ -# https://cookierobotics.com/007/ -cd(@__DIR__) -using Pkg -Pkg.activate("..") -using Revise, Distributions, ParameterSpacePartitions -using ParameterSpacePartitions.TestModels -using Random, DataFrames, StatsBase -#Random.seed!(554) - -# dimensions of the hypbercue -n_dims = 3 -# partitions per dimension -n_part = 2 -# number of starting points -n_start = 1 - -# partition boundaries -bounds = fill((0, 1), n_dims) -p_bounds = range(0, 1, n_part + 1) -hypercube = HyperCube(p_bounds) - -sample(bounds) = map(b -> rand(Uniform(b...)), bounds) - -init_parms = map(_ -> sample(bounds), 1:n_start) - -options = Options(; - radius = .40, - bounds, - n_iters = 500, - parallel = false, - init_parms, - t_rate = .3, - λ = .2 -) - -results = find_partitions( - model, - p_fun, - options, - hypercube -) - -df = DataFrame(results) - -transform!( - df, - :parms => identity => [:p1, :p2, :p3] -) - -transform!(df, :pattern => denserank => :pattern_id) - -groups = groupby(df, :chain_id) -mean_accept = combine(groups, :acceptance=>mean=>:mean) - -groups = groupby(df, :pattern) -summary = combine(groups, - :p1 => minimum, :p1 => maximum, - :p2 => minimum, :p2 => maximum, - :p3 => minimum, :p3 => maximum -) |> sort - -using GLMakie, ColorSchemes - -# transform pattern into integer id -transform!(df, :pattern => denserank => :pattern_id) - -fig, ax, scat = scatter( - df.p1, - df.p2, - df.p3, - color = df.pattern_id, - axis = (;type=Axis3), - markersize = 5000, - colormap = ColorSchemes.seaborn_deep6.colors, - grid = true -) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 7673a30..2eb5c9c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,613 +1,4 @@ using SafeTestsets -@safetestset "Process New Patterns" begin - using ParameterSpacePartitions, Test - import ParameterSpacePartitions: process_new_patterns!, Chain, is_new - - all_patterns = [[1,2],[1,3]] - patterns = [[1,2],[1,4]] - - chains = [Chain([.3,.3], [1,2], .2),Chain([.3,.3], [1,3], .2)] - chains[1].acceptance[1] = false - chains[2].acceptance[1] = false - - parms = [[.3,.4],[.3,.5]] - options = Options(; - radius = .2, - bounds = [(0,1),(0,1)], - n_iters = 100, - init_parms = [[.2,.3]] - ) - - @test !is_new(all_patterns, patterns[1]) - @test is_new(all_patterns, patterns[2]) - - process_new_patterns!(all_patterns, patterns, parms, chains, options) - - @test length(chains) == 3 - @test chains[3].pattern == [1,4] -end - -@safetestset "in_bounds" begin - using ParameterSpacePartitions, Test - import ParameterSpacePartitions: in_bounds, Chain - - bounds = [(0,1),(0,1)] - parms = [.5,.5] - in_bounds(parms, bounds) - @test in_bounds(parms, bounds) - - bounds = [(0,1),(0,1)] - parms = [.5,-1] - @test !in_bounds(parms, bounds) - - bounds = [(0,1),(0,1)] - parms = [1.2,.5] - @test !in_bounds(parms, bounds) -end - -@safetestset "update_position!" begin - using ParameterSpacePartitions - using Test - import ParameterSpacePartitions: update_position!, Chain - - parms = [.4] - proposal = [.3] - pattern = [4,3] - bounds = [(0,1)] - chain = Chain(parms, pattern, .3) - update_position!(chain, proposal, pattern, bounds) - @test chain.parms == [.3] - @test chain.pattern == [4,3] - @test chain.acceptance[2] == true - - parms = [.4] - proposal = [.3] - pattern = [4,3] - chain = Chain(parms, pattern, .3) - update_position!(chain, proposal, [4,5], bounds) - @test chain.parms == [.4] - @test chain.pattern == [4,3] - @test chain.acceptance[2] == false - - - parms = [.4] - proposal = [-.3] - pattern = [4,3] - chain = Chain(parms, pattern, .3) - update_position!(chain, proposal, pattern, bounds) - @test chain.parms == [.4] - @test chain.pattern == [4,3] - @test chain.acceptance[2] == false - -end - -@safetestset "Cube" begin - using ParameterSpacePartitions - using ParameterSpacePartitions.TestModels - using Test, Random, Distributions - using DataFrames - - Random.seed!(541) - # dimensions of the hypbercue - n_dims = 3 - # partitions per dimension - n_part = 4 - # number of starting points - n_start = 1 - - # partition boundaries - p_bounds = range(0, 1, n_part + 1) - hypercube = HyperCube(p_bounds) - bounds = fill((0, 1), n_dims) - - sample(bounds) = map(b -> rand(Uniform(b...)), bounds) - - init_parms = map(_ -> sample(bounds), 1:n_start) - - options = Options(; - radius = .18, - bounds, - n_iters = 100, - parallel = false, - init_parms - ) - - results = find_partitions( - model, - p_fun, - options, - hypercube - ) - - df = DataFrame(results) - transform!( - df, - :parms => identity => [:p1, :p2, :p3] - ) - - groups = groupby(df, :pattern) - combine(groups, - :p1 => minimum, :p1 => maximum, - :p2 => minimum, :p2 => maximum, - :p3 => minimum, :p3 => maximum - ) |> sort - - @test length(groups) == n_part^n_dims -end - - -@safetestset "5D Hypercube" begin - using ParameterSpacePartitions - using ParameterSpacePartitions.TestModels - using Test, Random, Distributions - using DataFrames - - Random.seed!(103) - # partitions per dimension - n_part = 4 - # dimensions of hypbercue - n_dims = 5 - # number of starting points - n_start = 1 - - # bounds of the hypercube - bounds = fill((0, 1), n_dims) - # partition boundaries - p_bounds = range(0, 1, n_part + 1) - hypercube = HyperCube(p_bounds) - - sample(bounds) = map(b -> rand(Uniform(b...)), bounds) - - init_parms = map(_ -> sample(bounds), 1:n_start) - - options = Options(; - radius = .15, - bounds, - n_iters = 100, - parallel = false, - init_parms - ) - - results = find_partitions( - model, - p_fun, - options, - hypercube - ) - - df = DataFrame(results) - transform!( - df, - :parms => identity => [:p1,:p2,:p3,:p4,:p5] - ) - - groups = groupby(df, :pattern) - combine(groups, - :p1 => minimum, :p1 => maximum, - :p2 => minimum, :p2 => maximum, - :p3 => minimum, :p3 => maximum, - :p4 => minimum, :p4 => maximum, - :p5 => minimum, :p5 => maximum - ) |> sort - - @test length(groups) == n_part^n_dims -end - -@safetestset "5D Polytope" begin - using ParameterSpacePartitions - using ParameterSpacePartitions.TestModels - using Test, Random, Distributions - using DataFrames, LinearAlgebra - Random.seed!(778) - - # dimensions of the hypbercue - n_dims = 5 - # number of partitions - n_part = 50 - # number of starting points - n_start = 1 - - - # partition boundaries - polytopes = [Polytope(rand(n_dims)) for i in 1:n_part] - bounds = fill((0, 1), n_dims) - - sample(bounds) = map(b -> rand(Uniform(b...)), bounds) - - init_parms = map(_ -> sample(bounds), 1:n_start) - - options = Options(; - radius = .1, - bounds, - n_iters = 500, - parallel = false, - init_parms - ) - - results = find_partitions( - model, - p_fun, - options, - polytopes - ) - - df = DataFrame(results) - transform!( - df, - :parms => identity => [:p1, :p2, :p3, :p4, :p5] - ) - - groups = groupby(df, :pattern) - @test length(groups) == n_part -end - -@safetestset "HyperSpheres" begin - using Test, ParameterSpacePartitions - using ParameterSpacePartitions.TestModels - using LinearAlgebra, Random, DataFrames, Distributions - Random.seed!(383) - - # dimensions of the hypbercue - n_dims = 7 - # number of hyperspheres - n_obj = 50 - # distribution of radii - r_dist = () -> rand(Uniform(.2, .25)) - # number of starting points - n_start = 10 - - # partition boundaries - hyperspheres = set_locations(r_dist, n_dims, n_obj) - bounds = fill((0, 1), n_dims) - - sample(bounds) = map(b -> rand(Uniform(b...)), bounds) - - init_parms = map(_ -> sample(bounds), 1:n_start) - - options = Options(; - radius = .15, - bounds, - n_iters = 5000, - init_parms, - λ = .05 - ) - - results = find_partitions( - model, - p_fun, - options, - hyperspheres - ) - - df = DataFrame(results) - groups = groupby(df, :pattern) - - @test length(groups) == (n_obj + 1) -end - -@safetestset "Adapt Radius 5D HyperCube" begin - using ParameterSpacePartitions - using ParameterSpacePartitions.TestModels - using Test, Random, Distributions - using DataFrames - - Random.seed!(587) - # partitions per dimension - n_part = 4 - # dimensions of hypbercue - n_dims = 5 - # number of starting points - n_start = 1 - - # bounds of the hypercube - bounds = fill((0, 1), n_dims) - # partition boundaries - p_bounds = range(0, 1, n_part + 1) - hypercube = HyperCube(p_bounds) - - sample(bounds) = map(b -> rand(Uniform(b...)), bounds) - - init_parms = map(_ -> sample(bounds), 1:n_start) - - t_rate = .4 - - options = Options(; - radius = .40, - bounds, - n_iters = 500, - parallel = false, - init_parms, - t_rate, - λ = .0 - ) - - results = find_partitions( - model, - p_fun, - options, - hypercube - ) - - df = DataFrame(results) - transform!( - df, - :parms => identity => [:p1,:p2,:p3,:p4,:p5] - ) - - groups = groupby(df, :chain_id) - mean_accept = combine(groups, :acceptance=>mean=>:mean) - # show that initial radius results in poor acceptance - @test mean(mean_accept.mean) < .10 - - options = Options(; - radius = .40, - bounds, - n_iters = 500, - parallel = false, - init_parms, - t_rate, - λ = .2 - ) - - results = find_partitions( - model, - p_fun, - options, - hypercube - ) - - df = DataFrame(results) - transform!( - df, - :parms => identity => [:p1,:p2,:p3,:p4,:p5] - ) - - groups = groupby(df, :chain_id) - mean_accept = combine(groups, :acceptance=>mean=>:mean) - @test mean(mean_accept.mean) ≈ t_rate atol = .03 - @test std(mean_accept.mean) < .05 -end - -@safetestset "Adapt Radius 5D Polytope" begin - using ParameterSpacePartitions - using ParameterSpacePartitions.TestModels - using Test, Random, Distributions - using DataFrames, LinearAlgebra - Random.seed!(2002) - - # dimensions of the hypbercue - n_dims = 5 - # number of partitions - n_part = 50 - # number of starting points - n_start = 1 - - - # partition boundaries - polytopes = [Polytope(rand(n_dims)) for i in 1:n_part] - bounds = fill((0, 1), n_dims) - - sample(bounds) = map(b -> rand(Uniform(b...)), bounds) - - init_parms = map(_ -> sample(bounds), 1:n_start) - - options = Options(; - radius = .5, - bounds, - n_iters = 500, - parallel = false, - init_parms, - adapt_radius! = no_adaption! - ) - - results = find_partitions( - model, - p_fun, - options, - polytopes - ) - - df = DataFrame(results) - transform!( - df, - :parms => identity => [:p1, :p2, :p3, :p4, :p5] - ) - - groups = groupby(df, :pattern) - @test length(groups) == n_part - - groups = groupby(df, :chain_id) - mean_accept = combine(groups, :acceptance=>mean=>:mean) - # show that initial radius results in poor acceptance - @test mean(mean_accept.mean) < .10 - - t_rate = .4 - - options = Options(; - radius = 1.0, - bounds, - n_iters = 1000, - parallel = false, - init_parms, - λ = 0.2, - t_rate - ) - - results = find_partitions( - model, - p_fun, - options, - polytopes - ) - - df = DataFrame(results) - transform!( - df, - :parms => identity => [:p1, :p2, :p3, :p4, :p5] - ) - - groups = groupby(df, :pattern) - @test length(groups) == n_part - - groups = groupby(df, :chain_id) - mean_accept = combine(groups, :acceptance=>mean=>:mean) - # show that initial radius results in poor acceptance - @test mean(mean_accept.mean) ≈ t_rate atol = .04 - @test std(mean_accept.mean) < .05 -end - -@safetestset "2D Ellipsoid Volume" begin - using Test, QHull, Random, Distributions, LinearAlgebra - using ParameterSpacePartitions: volume_ellipsoid, sample_ellipsoid_surface - Random.seed!(633) - n = 2 - μ = fill(0.0, n) - - x = randn(n, n) - cov_mat = x' * x - fun = sample_ellipsoid_surface - x = map(_ -> fun(μ, n, cov_mat), 1:10_000) - points = reduce(vcat, transpose.(x)) - - v1 = volume_ellipsoid(n, cov_mat) - - hull = chull(points) - v2 = hull.volume - - @test v1 ≈ v2 rtol = .03 - - x = randn(n, n) - cov_mat = x' * x - fun = sample_ellipsoid_surface - x = map(_ -> fun(μ, n, cov_mat), 1:10_000) - points = reduce(vcat, transpose.(x)) - - v1 = volume_ellipsoid(n, cov_mat) - - hull = chull(points) - v2 = hull.volume - - @test v1 ≈ v2 rtol = .03 -end - -@safetestset "3D Ellipsoid Volume" begin - using Test, QHull, Random, Distributions, LinearAlgebra - using ParameterSpacePartitions: volume_ellipsoid, sample_ellipsoid_surface - Random.seed!(584) - n = 3 - μ = fill(0.0, n) - - x = randn(n, n) - cov_mat = x' * x - fun = sample_ellipsoid_surface - x = map(_ -> fun(μ, n, cov_mat), 1:10_000) - points = reduce(vcat, transpose.(x)) - - v1 = volume_ellipsoid(n, cov_mat) - - hull = chull(points) - v2 = hull.volume - - @test v1 ≈ v2 rtol = .03 - - x = randn(n, n) - cov_mat = x' * x - fun = sample_ellipsoid_surface - x = map(_ -> fun(μ, n, cov_mat), 1:10_000) - points = reduce(vcat, transpose.(x)) - - v1 = volume_ellipsoid(n, cov_mat) - - hull = chull(points) - v2 = hull.volume - - @test v1 ≈ v2 rtol = .03 -end - -@safetestset "4D Ellipsoid Volume" begin - using Test, QHull, Random, Distributions, LinearAlgebra - using ParameterSpacePartitions: volume_ellipsoid, sample_ellipsoid_surface - Random.seed!(28805) - n = 4 - μ = fill(0.0, n) - - x = randn(n, n) - cov_mat = x' * x - fun = sample_ellipsoid_surface - x = map(_ -> fun(μ, n, cov_mat), 1:10_000) - points = reduce(vcat, transpose.(x)) - - v1 = volume_ellipsoid(n, cov_mat) - - hull = chull(points) - v2 = hull.volume - - @test v1 ≈ v2 rtol = .03 -end - -@safetestset "3D Odd Shape Volume" begin - - using Test, Distributions, ParameterSpacePartitions - using ParameterSpacePartitions.TestModels - using Random, DataFrames, StatsBase - include("volume_functions.jl") - - Random.seed!(2215) - - c = ( - # number of shapes - n_shapes = 2, - # number of simulations for hit or miss - n_sim = 10_000, - # number of dimensions - n_dims = 3, - # number of partitions per dimension - n_part = 5, - # number of cells for each shapes - n_cells = [10,20], - # number of starting points for the algorithm - n_start = 1, - ) - - ratios = map(_ -> volume_sim(c), 1:10) - - tv1 = c.n_cells[1] * (1 / c.n_part)^c.n_dims - tv2 = c.n_cells[2] * (1 / c.n_part)^c.n_dims - true_ratio = tv1 / tv2 - ratio = mean(ratios) - @test ratio ≈ true_ratio rtol = .2 -end - -@safetestset "5D Odd Shape Volume" begin - - using Test, Distributions, ParameterSpacePartitions - using ParameterSpacePartitions.TestModels - using Random, DataFrames, StatsBase - include("volume_functions.jl") - - Random.seed!(844) - - c = ( - # number of shapes - n_shapes = 2, - # number of simulations for hit or miss - n_sim = 10_000, - # number of dimensions - n_dims = 5, - # number of partitions per dimension - n_part = 5, - # number of cells for each shapes - n_cells = [10,20], - # number of starting points for the algorithm - n_start = 1, - ) - - ratios = map(_ -> volume_sim(c), 1:10) - - tv1 = c.n_cells[1] * (1 / c.n_part)^c.n_dims - tv2 = c.n_cells[2] * (1 / c.n_part)^c.n_dims - true_ratio = tv1 / tv2 - ratio = mean(ratios) - @test ratio ≈ true_ratio rtol = .5 -end \ No newline at end of file +include("sampler_tests.jl") +include("volume_tests.jl") \ No newline at end of file diff --git a/test/sampler_tests.jl b/test/sampler_tests.jl new file mode 100644 index 0000000..a2d6e70 --- /dev/null +++ b/test/sampler_tests.jl @@ -0,0 +1,456 @@ +@safetestset "Process New Patterns" begin + using ParameterSpacePartitions, Test + import ParameterSpacePartitions: process_new_patterns!, Chain, is_new + + all_patterns = [[1,2],[1,3]] + patterns = [[1,2],[1,4]] + + chains = [Chain([.3,.3], [1,2], .2),Chain([.3,.3], [1,3], .2)] + chains[1].acceptance[1] = false + chains[2].acceptance[1] = false + + parms = [[.3,.4],[.3,.5]] + options = Options(; + radius = .2, + bounds = [(0,1),(0,1)], + n_iters = 100, + init_parms = [[.2,.3]] + ) + + @test !is_new(all_patterns, patterns[1]) + @test is_new(all_patterns, patterns[2]) + + process_new_patterns!(all_patterns, patterns, parms, chains, options) + + @test length(chains) == 3 + @test chains[3].pattern == [1,4] +end + +@safetestset "in_bounds" begin + using ParameterSpacePartitions, Test + import ParameterSpacePartitions: in_bounds, Chain + + bounds = [(0,1),(0,1)] + parms = [.5,.5] + in_bounds(parms, bounds) + @test in_bounds(parms, bounds) + + bounds = [(0,1),(0,1)] + parms = [.5,-1] + @test !in_bounds(parms, bounds) + + bounds = [(0,1),(0,1)] + parms = [1.2,.5] + @test !in_bounds(parms, bounds) +end + +@safetestset "update_position!" begin + using ParameterSpacePartitions + using Test + import ParameterSpacePartitions: update_position!, Chain + + parms = [.4] + proposal = [.3] + pattern = [4,3] + bounds = [(0,1)] + chain = Chain(parms, pattern, .3) + update_position!(chain, proposal, pattern, bounds) + @test chain.parms == [.3] + @test chain.pattern == [4,3] + @test chain.acceptance[2] == true + + parms = [.4] + proposal = [.3] + pattern = [4,3] + chain = Chain(parms, pattern, .3) + update_position!(chain, proposal, [4,5], bounds) + @test chain.parms == [.4] + @test chain.pattern == [4,3] + @test chain.acceptance[2] == false + + + parms = [.4] + proposal = [-.3] + pattern = [4,3] + chain = Chain(parms, pattern, .3) + update_position!(chain, proposal, pattern, bounds) + @test chain.parms == [.4] + @test chain.pattern == [4,3] + @test chain.acceptance[2] == false + +end + +@safetestset "Cube" begin + using ParameterSpacePartitions + using ParameterSpacePartitions.TestModels + using Test, Random, Distributions + using DataFrames + + Random.seed!(541) + # dimensions of the hypbercue + n_dims = 3 + # partitions per dimension + n_part = 4 + # number of starting points + n_start = 1 + + # partition boundaries + p_bounds = range(0, 1, n_part + 1) + hypercube = HyperCube(p_bounds) + bounds = fill((0, 1), n_dims) + + sample(bounds) = map(b -> rand(Uniform(b...)), bounds) + + init_parms = map(_ -> sample(bounds), 1:n_start) + + options = Options(; + radius = .18, + bounds, + n_iters = 100, + parallel = false, + init_parms + ) + + results = find_partitions( + model, + p_fun, + options, + hypercube + ) + + df = DataFrame(results) + transform!( + df, + :parms => identity => [:p1, :p2, :p3] + ) + + groups = groupby(df, :pattern) + combine(groups, + :p1 => minimum, :p1 => maximum, + :p2 => minimum, :p2 => maximum, + :p3 => minimum, :p3 => maximum + ) |> sort + + @test length(groups) == n_part^n_dims +end + + +@safetestset "5D Hypercube" begin + using ParameterSpacePartitions + using ParameterSpacePartitions.TestModels + using Test, Random, Distributions + using DataFrames + + Random.seed!(103) + # partitions per dimension + n_part = 4 + # dimensions of hypbercue + n_dims = 5 + # number of starting points + n_start = 1 + + # bounds of the hypercube + bounds = fill((0, 1), n_dims) + # partition boundaries + p_bounds = range(0, 1, n_part + 1) + hypercube = HyperCube(p_bounds) + + sample(bounds) = map(b -> rand(Uniform(b...)), bounds) + + init_parms = map(_ -> sample(bounds), 1:n_start) + + options = Options(; + radius = .15, + bounds, + n_iters = 100, + parallel = false, + init_parms + ) + + results = find_partitions( + model, + p_fun, + options, + hypercube + ) + + df = DataFrame(results) + transform!( + df, + :parms => identity => [:p1,:p2,:p3,:p4,:p5] + ) + + groups = groupby(df, :pattern) + combine(groups, + :p1 => minimum, :p1 => maximum, + :p2 => minimum, :p2 => maximum, + :p3 => minimum, :p3 => maximum, + :p4 => minimum, :p4 => maximum, + :p5 => minimum, :p5 => maximum + ) |> sort + + @test length(groups) == n_part^n_dims +end + +@safetestset "5D Polytope" begin + using ParameterSpacePartitions + using ParameterSpacePartitions.TestModels + using Test, Random, Distributions + using DataFrames, LinearAlgebra + Random.seed!(778) + + # dimensions of the hypbercue + n_dims = 5 + # number of partitions + n_part = 50 + # number of starting points + n_start = 1 + + + # partition boundaries + polytopes = [Polytope(rand(n_dims)) for i in 1:n_part] + bounds = fill((0, 1), n_dims) + + sample(bounds) = map(b -> rand(Uniform(b...)), bounds) + + init_parms = map(_ -> sample(bounds), 1:n_start) + + options = Options(; + radius = .1, + bounds, + n_iters = 500, + parallel = false, + init_parms + ) + + results = find_partitions( + model, + p_fun, + options, + polytopes + ) + + df = DataFrame(results) + transform!( + df, + :parms => identity => [:p1, :p2, :p3, :p4, :p5] + ) + + groups = groupby(df, :pattern) + @test length(groups) == n_part +end + +@safetestset "HyperSpheres" begin + using Test, ParameterSpacePartitions + using ParameterSpacePartitions.TestModels + using LinearAlgebra, Random, DataFrames, Distributions + Random.seed!(383) + + # dimensions of the hypbercue + n_dims = 7 + # number of hyperspheres + n_obj = 50 + # distribution of radii + r_dist = () -> rand(Uniform(.2, .25)) + # number of starting points + n_start = 10 + + # partition boundaries + hyperspheres = set_locations(r_dist, n_dims, n_obj) + bounds = fill((0, 1), n_dims) + + sample(bounds) = map(b -> rand(Uniform(b...)), bounds) + + init_parms = map(_ -> sample(bounds), 1:n_start) + + options = Options(; + radius = .15, + bounds, + n_iters = 5000, + init_parms, + λ = .05 + ) + + results = find_partitions( + model, + p_fun, + options, + hyperspheres + ) + + df = DataFrame(results) + groups = groupby(df, :pattern) + + @test length(groups) == (n_obj + 1) +end + +@safetestset "Adapt Radius 5D HyperCube" begin + using ParameterSpacePartitions + using ParameterSpacePartitions.TestModels + using Test, Random, Distributions + using DataFrames + + Random.seed!(587) + # partitions per dimension + n_part = 4 + # dimensions of hypbercue + n_dims = 5 + # number of starting points + n_start = 1 + + # bounds of the hypercube + bounds = fill((0, 1), n_dims) + # partition boundaries + p_bounds = range(0, 1, n_part + 1) + hypercube = HyperCube(p_bounds) + + sample(bounds) = map(b -> rand(Uniform(b...)), bounds) + + init_parms = map(_ -> sample(bounds), 1:n_start) + + t_rate = .4 + + options = Options(; + radius = .40, + bounds, + n_iters = 500, + parallel = false, + init_parms, + t_rate, + λ = .0 + ) + + results = find_partitions( + model, + p_fun, + options, + hypercube + ) + + df = DataFrame(results) + transform!( + df, + :parms => identity => [:p1,:p2,:p3,:p4,:p5] + ) + + groups = groupby(df, :chain_id) + mean_accept = combine(groups, :acceptance=>mean=>:mean) + # show that initial radius results in poor acceptance + @test mean(mean_accept.mean) < .10 + + options = Options(; + radius = .40, + bounds, + n_iters = 500, + parallel = false, + init_parms, + t_rate, + λ = .2 + ) + + results = find_partitions( + model, + p_fun, + options, + hypercube + ) + + df = DataFrame(results) + transform!( + df, + :parms => identity => [:p1,:p2,:p3,:p4,:p5] + ) + + groups = groupby(df, :chain_id) + mean_accept = combine(groups, :acceptance=>mean=>:mean) + @test mean(mean_accept.mean) ≈ t_rate atol = .03 + @test std(mean_accept.mean) < .05 +end + +@safetestset "Adapt Radius 5D Polytope" begin + using ParameterSpacePartitions + using ParameterSpacePartitions.TestModels + using Test, Random, Distributions + using DataFrames, LinearAlgebra + Random.seed!(2002) + + # dimensions of the hypbercue + n_dims = 5 + # number of partitions + n_part = 50 + # number of starting points + n_start = 1 + + + # partition boundaries + polytopes = [Polytope(rand(n_dims)) for i in 1:n_part] + bounds = fill((0, 1), n_dims) + + sample(bounds) = map(b -> rand(Uniform(b...)), bounds) + + init_parms = map(_ -> sample(bounds), 1:n_start) + + options = Options(; + radius = .5, + bounds, + n_iters = 500, + parallel = false, + init_parms, + adapt_radius! = no_adaption! + ) + + results = find_partitions( + model, + p_fun, + options, + polytopes + ) + + df = DataFrame(results) + transform!( + df, + :parms => identity => [:p1, :p2, :p3, :p4, :p5] + ) + + groups = groupby(df, :pattern) + @test length(groups) == n_part + + groups = groupby(df, :chain_id) + mean_accept = combine(groups, :acceptance=>mean=>:mean) + # show that initial radius results in poor acceptance + @test mean(mean_accept.mean) < .10 + + t_rate = .4 + + options = Options(; + radius = 1.0, + bounds, + n_iters = 1000, + parallel = false, + init_parms, + λ = 0.2, + t_rate + ) + + results = find_partitions( + model, + p_fun, + options, + polytopes + ) + + df = DataFrame(results) + transform!( + df, + :parms => identity => [:p1, :p2, :p3, :p4, :p5] + ) + + groups = groupby(df, :pattern) + @test length(groups) == n_part + + groups = groupby(df, :chain_id) + mean_accept = combine(groups, :acceptance=>mean=>:mean) + # show that initial radius results in poor acceptance + @test mean(mean_accept.mean) ≈ t_rate atol = .04 + @test std(mean_accept.mean) < .05 +end \ No newline at end of file diff --git a/test/volume_tests.jl b/test/volume_tests.jl new file mode 100644 index 0000000..ef22b8e --- /dev/null +++ b/test/volume_tests.jl @@ -0,0 +1,249 @@ +@safetestset "2D Ellipsoid Volume" begin + using Test, QHull, Random, Distributions, LinearAlgebra + using ParameterSpacePartitions: volume_ellipsoid, sample_ellipsoid_surface + Random.seed!(633) + n = 2 + μ = fill(0.0, n) + + x = randn(n, n) + cov_mat = x' * x + fun = sample_ellipsoid_surface + x = map(_ -> fun(μ, n, cov_mat), 1:10_000) + points = reduce(vcat, transpose.(x)) + + v1 = volume_ellipsoid(n, cov_mat) + + hull = chull(points) + v2 = hull.volume + + @test v1 ≈ v2 rtol = .03 + + x = randn(n, n) + cov_mat = x' * x + fun = sample_ellipsoid_surface + x = map(_ -> fun(μ, n, cov_mat), 1:10_000) + points = reduce(vcat, transpose.(x)) + + v1 = volume_ellipsoid(n, cov_mat) + + hull = chull(points) + v2 = hull.volume + + @test v1 ≈ v2 rtol = .03 +end + +@safetestset "3D Ellipsoid Volume" begin + using Test, QHull, Random, Distributions, LinearAlgebra + using ParameterSpacePartitions: volume_ellipsoid, sample_ellipsoid_surface + Random.seed!(584) + n = 3 + μ = fill(0.0, n) + + x = randn(n, n) + cov_mat = x' * x + fun = sample_ellipsoid_surface + x = map(_ -> fun(μ, n, cov_mat), 1:10_000) + points = reduce(vcat, transpose.(x)) + + v1 = volume_ellipsoid(n, cov_mat) + + hull = chull(points) + v2 = hull.volume + + @test v1 ≈ v2 rtol = .03 + + x = randn(n, n) + cov_mat = x' * x + fun = sample_ellipsoid_surface + x = map(_ -> fun(μ, n, cov_mat), 1:10_000) + points = reduce(vcat, transpose.(x)) + + v1 = volume_ellipsoid(n, cov_mat) + + hull = chull(points) + v2 = hull.volume + + @test v1 ≈ v2 rtol = .03 +end + +@safetestset "4D Ellipsoid Volume" begin + using Test, QHull, Random, Distributions, LinearAlgebra + using ParameterSpacePartitions: volume_ellipsoid, sample_ellipsoid_surface + Random.seed!(28805) + n = 4 + μ = fill(0.0, n) + + x = randn(n, n) + cov_mat = x' * x + fun = sample_ellipsoid_surface + x = map(_ -> fun(μ, n, cov_mat), 1:10_000) + points = reduce(vcat, transpose.(x)) + + v1 = volume_ellipsoid(n, cov_mat) + + hull = chull(points) + v2 = hull.volume + + @test v1 ≈ v2 rtol = .03 +end + +@safetestset "3D Odd Shape Volume" begin + + using Test, Distributions, ParameterSpacePartitions + using ParameterSpacePartitions.TestModels + using Random, DataFrames, StatsBase + include("volume_functions.jl") + + Random.seed!(2215) + + c = ( + # number of shapes + n_shapes = 2, + # number of simulations for hit or miss + n_sim = 10_000, + # number of dimensions + n_dims = 3, + # number of partitions per dimension + n_part = 5, + # number of cells for each shapes + n_cells = [10,20], + # number of starting points for the algorithm + n_start = 1, + ) + + ratios = map(_ -> volume_sim(c), 1:10) + + tv1 = c.n_cells[1] * (1 / c.n_part)^c.n_dims + tv2 = c.n_cells[2] * (1 / c.n_part)^c.n_dims + true_ratio = tv1 / tv2 + ratio = mean(ratios) + @test ratio ≈ true_ratio rtol = .2 +end + +@safetestset "5D Odd Shape Volume" begin + + using Test, Distributions, ParameterSpacePartitions + using ParameterSpacePartitions.TestModels + using Random, DataFrames, StatsBase + include("volume_functions.jl") + + Random.seed!(844) + + c = ( + # number of shapes + n_shapes = 2, + # number of simulations for hit or miss + n_sim = 10_000, + # number of dimensions + n_dims = 5, + # number of partitions per dimension + n_part = 5, + # number of cells for each shapes + n_cells = [10,20], + # number of starting points for the algorithm + n_start = 1, + ) + + ratios = map(_ -> volume_sim(c), 1:10) + + tv1 = c.n_cells[1] * (1 / c.n_part)^c.n_dims + tv2 = c.n_cells[2] * (1 / c.n_part)^c.n_dims + true_ratio = tv1 / tv2 + ratio = mean(ratios) + @test ratio ≈ true_ratio rtol = .5 +end + +@safetestset "Volume 5D Polytope" begin + + using Test, ParameterSpacePartitions + using ParameterSpacePartitions.TestModels + using Random, DataFrames, Distributions + Random.seed!(8492) + + # dimensions of the hypbercue + n_dims = 5 + # number of partitions + n_part = 5 + # number of starting points + n_start = 1 + + #points = [round.(rand(n_dims), digits=3) for _ in 1:n_part] + + points = [ + [0.659, 0.547, 0.842, 0.349, 0.667], + [0.121, 0.758, 0.141, 0.4, 0.955], + [0.232, 0.696, 0.999, 0.655, 0.257], + [0.312, 0.057, 0.977, 0.025, 0.154], + [0.191, 0.782, 0.051, 0.697, 0.653], + ] + + # partition boundaries + polytopes = Polytope.(points) + bounds = fill((0, 1), n_dims) + + sample(bounds) = map(b -> rand(Uniform(b...)), bounds) + + init_parms = map(_ -> sample(bounds), 1:n_start) + + options = Options(; + radius = .1, + bounds, + n_iters = 10_000, + parallel = false, + init_parms, + λ = .2, + t_rate = .3, + ) + + results = find_partitions( + model, + p_fun, + options, + polytopes + ) + + parm_names = Symbol.("p", 1:n_dims) + df = DataFrame(results) + transform!( + df, + :parms => identity => parm_names + ) + + groups = groupby(df, :chain_id) + mean_accept = combine(groups, :acceptance=>mean=>:mean) + + groups = groupby(df, :pattern) + + df_volume = combine( + groups, + x -> estimate_volume( + model, + p_fun, + x, + bounds, + polytopes; + parm_names, + ) + ) + + df_volume.volume = df_volume.x1 / sum(df_volume.x1) + std(df_volume.volume) + # based on 10M SMC samples + true_volume = [ + 0.40331 + 0.10566 + 0.17438 + 0.0848 + 0.23185 + ] + true_df = DataFrame(pattern = 1:n_part, volume = true_volume) + abs_error = abs.(true_volume .- df_volume.volume) + max_abs_error = maximum(abs_error) + mean_abs_error = mean(abs_error) + sort!(df_volume, :volume) + sort!(true_df, :volume) + r = cor(df_volume.pattern, true_df.pattern) + @test max_abs_error < .03 + @test mean_abs_error < .015 + @test r > .8 +end \ No newline at end of file