From b2dc51972229bf99b3983a0dd06a30be8084c9cd Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 25 Nov 2020 02:33:35 +0000 Subject: [PATCH 01/14] initial structure of Parameter Distributions, and their constraints --- src/ParameterDistributions.jl | 256 ++++++++++++++++++++++++++++++++++ 1 file changed, 256 insertions(+) create mode 100644 src/ParameterDistributions.jl diff --git a/src/ParameterDistributions.jl b/src/ParameterDistributions.jl new file mode 100644 index 000000000..444760d12 --- /dev/null +++ b/src/ParameterDistributions.jl @@ -0,0 +1,256 @@ +module ParameterDistribution + +## Imports +using Distributions +using StatsBase +using Random + +## Exports + +#objects +export ParameterDistribution + +#functions +export set_distribution, get_distribution, sample_distribution +export transform_real_to_prior, transform_prior_to_real +#export apply_units_to_real + + +## Objects +# for the Distribution +abstract type ParameterDistributionType end + +""" + Parameterized <: ParameterDistributionType + +A distribution constructed from a parametrized formula (e.g Julia Distributions.jl) +""" +struct Parameterized <: ParameterDistributionType + distribution::Distribution +end + +""" + Samples{FT<:AbstractFloat} <: ParameterDistributionType + +A distribution comprised of only samples +""" +struct Samples{FT<:AbstractFloat} <: ParameterDistributionType + distribution_samples::Array{FT} +end + +# """ +# KernelDensityEstimate <: ParameterDistributionType +# +# A distribution comprised of a Kernel Density Estimate of samples +# """ +# struct KernelDensityEstimate <: ParameterDistributionType +# # Not implemented +# end + +# For the transforms +abstract type ConstraintType end + +""" + NoConstrain <: ConstraintType + +The mapping x -> x +""" + +struct NoConstraint <: ConstraintType end + +""" + BoundedBelow{FT <: AbstractFloat} <: ConstraintType + +The mapping x -> ln(x - lower_bound{FT}); from a space bounded below to an unbounded space. +""" +struct BoundedBelow{FT <: AbstractFloat} <: + lower_bound::FT +end + +""" + BoundedAbove{FT <: AbstractFloat} <: ConstraintType + +The mapping to x -> ln( upper_bound{FT} - x); from a space bounded above to an unbounded space. +""" +struct BoundedAbove{FT <: AbstractFloat} + upper_bound::FT +end + +""" + Bounded{FT <: AbstractFloat} <: ConstraintType + +The mapping x -> ln( (x - lower_bound) / (upper_bound - x) ); from a bounded space to an unbounded space. +""" +struct Bounded{FT <: AbstractFloat} + lower_bound::FT + upper_bound::FT +end + +function Bounded(lower::FT, upper::FT) where {FT <: AbstractFloat} + @assert lower_bound < upper_bound + Bounded{FT}(lower_bound,upper_bound) +end + +""" + struct ParameterDistribution + +Structure to hold a parameter distribution +""" +struct ParameterDistribution{FT<:AbstractFloat, + PDType <: ParameterDistributionType, + CType <: ConstraintType} + distribution::PDType + constraints::Array{CType} + name::String +end + +function ParameterDistribution(parameter_distribution::PDType, + constraints::Array{CType}, + name::String) where {FT <: AbstractFloat, + PDType <: ParameterDistributionType, + CType <: ConstraintType} + #assert that total parameter dimensions matches number of constraints + length(parameter_distribution) == length(constraints) || + throw(DimensionMismatch("There must be one constraint per parameter")) + + ParameterDistribution{FT,PDType,CType}(parameter_distributions, constraints, name) +end + +""" + struct ParameterDistributions + +Structure to hold an array of ParameterDistribution's +""" +struct ParameterDistributions + parameter_distribution::Array{ParameterDistribution} +end + + + + +## Functions + +""" + get_names(pds::ParameterDistributions) + +Returns a list of ParameterDistribution names +""" +function get_name(pds::ParameterDistributions) + return [get_name(pd) for pd in pds] +end + +function get_name(pd::ParameterDistribution) + return pd.name +end + +""" + get_distributions(pds::ParameterDistributions) + +Returns a `Dict` of `ParameterDistribution` distributions by name, (unless sample type) +""" +function get_distribution(pds::ParameterDistributions) + return Dict{String,Any}(get_name(pd) => get_distribution(pd) for pd in pds) +end + +function get_distribution(pd::ParameterDistribution{FT,Samples,CType}) where {FT <: AbstractFloat, + CType <: ConstraintType} + return "Contains samples only" +end +function get_distribution(pd::ParameterDistribution{FT,Parameterized,CType}) where {FT <: AbstractFloat, + CType <: ConstraintType} + return pd.distribution.distribution +end + +""" + function sample_distribution(pds::ParameterDistributions) + +Draws samples from the parameter distributions +""" +function sample_distribution(pds::ParameterDistributions) + return sample_distribution(pds,1) +end + +function sample_distribution(pds::ParameterDistributions, n_samples::IT) where {IT <: Integer} + return Dict{String,Any}(get_name(pd) => sample_distribution(pd, n_samples) for pd in pds) +end + +function sample_distribution(pd::ParameterDistribution) + return sample_distribution(pd,1) +end + +function sample_distribution(pd::ParameterDistribution, n_samples::IT) where {IT <: Integer} + return sample_distribution(pd.distribution) +end + +function sample_distribution(d::Samples, n_samples::IT) where {IT <: Integer} + total_samples = size(d.distribution_samples)[1] + samples_idx = sample(collect(1:total_samples), n_samples ; replace=false) + return d.distribution_samples[samples_idx,:] +end + +function sample_distribution(d::Parameterized,n_samples::IT) + return rand(pd.distirbution.distribution, n_samples) +end + + +#transforms keyword template +""" + transform_real_to_prior(pds::ParameterDistributions, x::Array{AbstractFloat}) + +apply the transformation to map real parameter x into the unbounded prior space +""" +function transform_real_to_prior(pds::ParameterDistributions, x::Array{FT}) where {FT <: AbstractFloat} + return Dict{String,Any}(get_name(pd) => transform_real_to_prior(pd,x) for pd in pds) +end + +function transform_real_to_prior(pd::ParameterDistribution, x::Array{FT}) where {FT <: AbstractFloat} + return [transform_real_to_prior(constraint,x) for constraint in pd.constraints] +end + +function transform_real_to_prior(c::NoConstraints , x::Array{FT}) where {FT <: AbstractFloat} + return x +end +function transform_real_to_prior(c::BoundedBelow, x::Array{FT}) where {FT <: AbstractFloat} + return ln(x - c.lower_bound) +end +function transform_real_to_prior(c::BoundedAbove, x::Array{FT}) where {FT <: AbstractFloat} + return ln(c.upper_bound - x) +end + +function transform_real_to_prior(c::Bounded, x::Array{FT}) where {FT <: AbstractFloat} + return ln((c.upper_bound - x) / (x - c.lower_bound)) +end + + + +""" + transform_prior_to_real(pds::ParameterDistributions, x::Array{AbstractFloat}) + +apply the transformation to map parameter x into the real space +""" +function transform_prior_to_real(pds::ParameterDistributions, x::Array{FT}) where {FT <: AbstractFloat} + return Dict{String,Any}(get_name(pd) => transform_prior_to_real(pd,x) for pd in pds) +end + +function transform_prior_to_real(pd::ParameterDistributions, x::Array{FT}) where {FT <: AbstractFloat} + return [transform_prior_to_real(constraint,x) for constraint in pd.constraints] +end + + +function transform_prior_to_real(c::NoConstraints , x::Array{FT}) where {FT <: AbstractFloat} + return x +end +function transform_prior_to_real(c::BoundedBelow, x::Array{FT}) where {FT <: AbstractFloat} + return exp(x) + c.lower_bound +end +function transform_prior_to_real(c::BoundedAbove, x::Array{FT}) where {FT <: AbstractFloat} + return c.upper_bound - exp(x) +end + +function transform_prior_to_real(c::Bounded, x::Array{FT}) where {FT <: AbstractFloat} + return (c.upper_bound * exp(x) + c.lower_bound) / (exp(x) + 1) +end + + + +end # of module From 2c3d420e4154c9725c13c6105b8670830a2fffd4 Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 25 Nov 2020 11:43:37 +0000 Subject: [PATCH 02/14] properly split parameters in transforms --- src/ParameterDistributions.jl | 56 +++++++++++++++++++++++++---------- 1 file changed, 40 insertions(+), 16 deletions(-) diff --git a/src/ParameterDistributions.jl b/src/ParameterDistributions.jl index 444760d12..1cd8e6b25 100644 --- a/src/ParameterDistributions.jl +++ b/src/ParameterDistributions.jl @@ -193,18 +193,32 @@ function sample_distribution(d::Parameterized,n_samples::IT) end -#transforms keyword template +#apply transforms +function length(c::CType) where {CType <: ConstraintType} + return 1 +end +function length(carray::Array{CType}) where {CType <: ConstraintType} + return size(carray)[1] +end + """ transform_real_to_prior(pds::ParameterDistributions, x::Array{AbstractFloat}) -apply the transformation to map real parameter x into the unbounded prior space +Apply the transformation to map real parameter x into the unbounded prior space """ -function transform_real_to_prior(pds::ParameterDistributions, x::Array{FT}) where {FT <: AbstractFloat} - return Dict{String,Any}(get_name(pd) => transform_real_to_prior(pd,x) for pd in pds) +function transform_real_to_prior(pds::ParameterDistributions, xarray::Array{FT}) where {FT <: AbstractFloat} + #split xarray into chunks to be fed into each distribution + clen = [length(pd.constraints) for pd in pds] # e.g [2,3,1,3] + cumulative_clen = [sum(clen[1:i]) for i = 1:size(clen)[1]] # e.g [1 1:2, 3:5, 6, 7:9 ] + x_idx = Dict(Integer,Array{Integer})(i => collect(cumulative_cl[i-1]+1:cumulative_cl[i]) + for i in 2:size(cumulative_clen)[1]) + x_idx[1] = collect(1:cumulative_clen[1]) + + return Dict{String,Any}(get_name(pd) => transform_real_to_prior(pd,xarray[x_idx[i]]) for (i,pd) in enumerate(pds)) end -function transform_real_to_prior(pd::ParameterDistribution, x::Array{FT}) where {FT <: AbstractFloat} - return [transform_real_to_prior(constraint,x) for constraint in pd.constraints] +function transform_real_to_prior(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: AbstractFloat} + return [transform_real_to_prior(constraint,xarray[i]) for (i,constraint) in enumerate(pd.constraints)] end function transform_real_to_prior(c::NoConstraints , x::Array{FT}) where {FT <: AbstractFloat} @@ -224,30 +238,40 @@ end """ - transform_prior_to_real(pds::ParameterDistributions, x::Array{AbstractFloat}) + transform_prior_to_real(pds::ParameterDistributions, xarray::Array{AbstractFloat}) -apply the transformation to map parameter x into the real space +Apply the transformation to map parameter x into the real space """ -function transform_prior_to_real(pds::ParameterDistributions, x::Array{FT}) where {FT <: AbstractFloat} - return Dict{String,Any}(get_name(pd) => transform_prior_to_real(pd,x) for pd in pds) +function transform_prior_to_real(pds::ParameterDistributions, + xarray::Array{FT}) where {FT <: AbstractFloat} + + #split xarray into chunks to be fed into each distribution + clen = [length(pd.constraints) for pd in pds] # e.g [2,3,1,3] + cumulative_clen = [sum(clen[1:i]) for i = 1:size(clen)[1]] # e.g [1 1:2, 3:5, 6, 7:9 ] + x_idx = Dict(Integer,Array{Integer})(i => collect(cumulative_cl[i-1]+1:cumulative_cl[i]) + for i in 2:size(cumulative_clen)[1]) + x_idx[1] = collect(1:cumulative_clen[1]) + + return Dict{String,Any}(get_name(pd) => transform_prior_to_real(pd,xarray[x_idx[i]]) + for (i,pd) in enumerate(pds)) end -function transform_prior_to_real(pd::ParameterDistributions, x::Array{FT}) where {FT <: AbstractFloat} - return [transform_prior_to_real(constraint,x) for constraint in pd.constraints] +function transform_prior_to_real(pd::ParameterDistributions, xarray::Array{FT}) where {FT <: AbstractFloat} + return [transform_prior_to_real(constraint,xarray[i]) for (i,constraint) in enumerate(pd.constraints)] end -function transform_prior_to_real(c::NoConstraints , x::Array{FT}) where {FT <: AbstractFloat} +function transform_prior_to_real(c::NoConstraints , x::FT) where {FT <: AbstractFloat} return x end -function transform_prior_to_real(c::BoundedBelow, x::Array{FT}) where {FT <: AbstractFloat} +function transform_prior_to_real(c::BoundedBelow, x::FT) where {FT <: AbstractFloat} return exp(x) + c.lower_bound end -function transform_prior_to_real(c::BoundedAbove, x::Array{FT}) where {FT <: AbstractFloat} +function transform_prior_to_real(c::BoundedAbove, x::FT) where {FT <: AbstractFloat} return c.upper_bound - exp(x) end -function transform_prior_to_real(c::Bounded, x::Array{FT}) where {FT <: AbstractFloat} +function transform_prior_to_real(c::Bounded, x::FT) where {FT <: AbstractFloat} return (c.upper_bound * exp(x) + c.lower_bound) / (exp(x) + 1) end From 6b85f3af1685ee29a512b51690748728fe78771d Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 25 Nov 2020 11:52:21 +0000 Subject: [PATCH 03/14] docstrings --- src/ParameterDistributions.jl | 43 ++++++++++++++++++++++++++++------- 1 file changed, 35 insertions(+), 8 deletions(-) diff --git a/src/ParameterDistributions.jl b/src/ParameterDistributions.jl index 1cd8e6b25..2665c4e7d 100644 --- a/src/ParameterDistributions.jl +++ b/src/ParameterDistributions.jl @@ -53,7 +53,7 @@ abstract type ConstraintType end """ NoConstrain <: ConstraintType -The mapping x -> x +No constraint. """ struct NoConstraint <: ConstraintType end @@ -61,7 +61,7 @@ struct NoConstraint <: ConstraintType end """ BoundedBelow{FT <: AbstractFloat} <: ConstraintType -The mapping x -> ln(x - lower_bound{FT}); from a space bounded below to an unbounded space. +A lower bound constraint. """ struct BoundedBelow{FT <: AbstractFloat} <: lower_bound::FT @@ -70,7 +70,7 @@ end """ BoundedAbove{FT <: AbstractFloat} <: ConstraintType -The mapping to x -> ln( upper_bound{FT} - x); from a space bounded above to an unbounded space. +And upper bound constraint """ struct BoundedAbove{FT <: AbstractFloat} upper_bound::FT @@ -79,7 +79,7 @@ end """ Bounded{FT <: AbstractFloat} <: ConstraintType -The mapping x -> ln( (x - lower_bound) / (upper_bound - x) ); from a bounded space to an unbounded space. +Both a lower and upper bound constraint. """ struct Bounded{FT <: AbstractFloat} lower_bound::FT @@ -87,7 +87,7 @@ struct Bounded{FT <: AbstractFloat} end function Bounded(lower::FT, upper::FT) where {FT <: AbstractFloat} - @assert lower_bound < upper_bound + lower_bound < upper_bound || throw(DomainError("upper bound must be greater than lower bound")) Bounded{FT}(lower_bound,upper_bound) end @@ -204,7 +204,7 @@ end """ transform_real_to_prior(pds::ParameterDistributions, x::Array{AbstractFloat}) -Apply the transformation to map real parameter x into the unbounded prior space +Apply the transformation to map real (and possibly constrained) parameters `xarray` into the unbounded prior space """ function transform_real_to_prior(pds::ParameterDistributions, xarray::Array{FT}) where {FT <: AbstractFloat} #split xarray into chunks to be fed into each distribution @@ -221,16 +221,30 @@ function transform_real_to_prior(pd::ParameterDistribution, xarray::Array{FT}) w return [transform_real_to_prior(constraint,xarray[i]) for (i,constraint) in enumerate(pd.constraints)] end +""" +No constraint mapping x -> x +""" function transform_real_to_prior(c::NoConstraints , x::Array{FT}) where {FT <: AbstractFloat} return x end + +""" +Bounded below -> unbounded, use mapping x -> ln(x - lower_bound) +""" function transform_real_to_prior(c::BoundedBelow, x::Array{FT}) where {FT <: AbstractFloat} return ln(x - c.lower_bound) end + +""" +Bounded above -> unbounded, use mapping x -> ln(upper_bound - x) +""" function transform_real_to_prior(c::BoundedAbove, x::Array{FT}) where {FT <: AbstractFloat} return ln(c.upper_bound - x) end +""" +Bounded -> unbounded, use mapping x -> ln((x - lower_bound) / (upper_bound - x) +""" function transform_real_to_prior(c::Bounded, x::Array{FT}) where {FT <: AbstractFloat} return ln((c.upper_bound - x) / (x - c.lower_bound)) end @@ -240,7 +254,7 @@ end """ transform_prior_to_real(pds::ParameterDistributions, xarray::Array{AbstractFloat}) -Apply the transformation to map parameter x into the real space +Apply the transformation to map parameters `xarray` from the unbounded space into (possibly constrained) real space """ function transform_prior_to_real(pds::ParameterDistributions, xarray::Array{FT}) where {FT <: AbstractFloat} @@ -260,17 +274,30 @@ function transform_prior_to_real(pd::ParameterDistributions, xarray::Array{FT}) return [transform_prior_to_real(constraint,xarray[i]) for (i,constraint) in enumerate(pd.constraints)] end - +""" +No constraint mapping x -> x +""" function transform_prior_to_real(c::NoConstraints , x::FT) where {FT <: AbstractFloat} return x end + +""" +Unbounded -> bounded below, use mapping x -> exp(x) + lower_bound +""" function transform_prior_to_real(c::BoundedBelow, x::FT) where {FT <: AbstractFloat} return exp(x) + c.lower_bound end + +""" +Unbounded -> bounded above, use mapping x -> upper_bound - exp(x) +""" function transform_prior_to_real(c::BoundedAbove, x::FT) where {FT <: AbstractFloat} return c.upper_bound - exp(x) end +""" +Unbounded -> bounded, use mapping x -> (upper_bound * exp(x) + lower_bound) / (exp(x) + 1) +""" function transform_prior_to_real(c::Bounded, x::FT) where {FT <: AbstractFloat} return (c.upper_bound * exp(x) + c.lower_bound) / (exp(x) + 1) end From 828ef55ce0c858378de79ee827758cbaf2850a97 Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 25 Nov 2020 16:59:30 +0000 Subject: [PATCH 04/14] added Unit tests to objects and functions --- src/CalibrateEmulateSample.jl | 1 + src/ParameterDistributions.jl | 162 ++++++++++++------------ test/ParameterDistributions/runtests.jl | 132 +++++++++++++++++++ test/runtests.jl | 3 +- 4 files changed, 217 insertions(+), 81 deletions(-) create mode 100644 test/ParameterDistributions/runtests.jl diff --git a/src/CalibrateEmulateSample.jl b/src/CalibrateEmulateSample.jl index 0871ff4e9..0a8720c70 100644 --- a/src/CalibrateEmulateSample.jl +++ b/src/CalibrateEmulateSample.jl @@ -5,6 +5,7 @@ using Distributions, Statistics, LinearAlgebra, DocStringExtensions # No internal deps, light external deps include("Observations.jl") include("Priors.jl") +include("ParameterDistributions.jl") # No internal deps, heavy external deps include("EKP.jl") diff --git a/src/ParameterDistributions.jl b/src/ParameterDistributions.jl index 2665c4e7d..34137356f 100644 --- a/src/ParameterDistributions.jl +++ b/src/ParameterDistributions.jl @@ -1,4 +1,4 @@ -module ParameterDistribution +module ParameterDistributionsStorage ## Imports using Distributions @@ -8,13 +8,14 @@ using Random ## Exports #objects -export ParameterDistribution +export Parameterized, Samples +export NoConstraint, BoundedBelow, BoundedAbove, Bounded +export ParameterDistribution, ParameterDistributions #functions -export set_distribution, get_distribution, sample_distribution +export get_name, get_distribution +export sample_distribution export transform_real_to_prior, transform_prior_to_real -#export apply_units_to_real - ## Objects # for the Distribution @@ -30,28 +31,20 @@ struct Parameterized <: ParameterDistributionType end """ - Samples{FT<:AbstractFloat} <: ParameterDistributionType + Samples{FT<:Real} <: ParameterDistributionType A distribution comprised of only samples """ -struct Samples{FT<:AbstractFloat} <: ParameterDistributionType +struct Samples{FT<:Real} <: ParameterDistributionType distribution_samples::Array{FT} end -# """ -# KernelDensityEstimate <: ParameterDistributionType -# -# A distribution comprised of a Kernel Density Estimate of samples -# """ -# struct KernelDensityEstimate <: ParameterDistributionType -# # Not implemented -# end # For the transforms abstract type ConstraintType end """ - NoConstrain <: ConstraintType + NoConstraint <: ConstraintType No constraint. """ @@ -59,36 +52,56 @@ No constraint. struct NoConstraint <: ConstraintType end """ - BoundedBelow{FT <: AbstractFloat} <: ConstraintType + BoundedBelow{FT <: Real} <: ConstraintType A lower bound constraint. """ -struct BoundedBelow{FT <: AbstractFloat} <: +struct BoundedBelow{FT <: Real} <: ConstraintType lower_bound::FT end - """ - BoundedAbove{FT <: AbstractFloat} <: ConstraintType + BoundedAbove{FT <: Real} <: ConstraintType And upper bound constraint """ -struct BoundedAbove{FT <: AbstractFloat} +struct BoundedAbove{FT <: Real} <: ConstraintType upper_bound::FT end """ - Bounded{FT <: AbstractFloat} <: ConstraintType + Bounded{FT <: Real} <: ConstraintType Both a lower and upper bound constraint. """ -struct Bounded{FT <: AbstractFloat} +struct Bounded{FT <: Real} <: ConstraintType lower_bound::FT upper_bound::FT + Bounded(lower_bound::FT,upper_bound::FT) where {FT <: Real} = upper_bound <= lower_bound ? throw(DomainError("upper bound must be greater than lower bound")) : new{FT}(lower_bound, upper_bound) end -function Bounded(lower::FT, upper::FT) where {FT <: AbstractFloat} - lower_bound < upper_bound || throw(DomainError("upper bound must be greater than lower bound")) - Bounded{FT}(lower_bound,upper_bound) +""" + len(c::Array{CType}) + +The number of constraints, each constraint has length 1. +""" +function len(c::CType) where {CType <: ConstraintType} + return 1 +end + +function len(carray::Array{CType}) where {CType <: ConstraintType} + return size(carray)[1] +end +""" + len(d::ParametrizedDistributionType) + +The number of dimensions of the parameter space +""" +function len(d::Parameterized) + return length(d.distribution) +end + +function len(d::Samples) + return size(d.distribution_samples)[2] end """ @@ -96,25 +109,18 @@ end Structure to hold a parameter distribution """ -struct ParameterDistribution{FT<:AbstractFloat, - PDType <: ParameterDistributionType, - CType <: ConstraintType} +struct ParameterDistribution{PDType <: ParameterDistributionType, CType <: ConstraintType} distribution::PDType constraints::Array{CType} name::String -end -function ParameterDistribution(parameter_distribution::PDType, - constraints::Array{CType}, - name::String) where {FT <: AbstractFloat, - PDType <: ParameterDistributionType, - CType <: ConstraintType} - #assert that total parameter dimensions matches number of constraints - length(parameter_distribution) == length(constraints) || - throw(DimensionMismatch("There must be one constraint per parameter")) - - ParameterDistribution{FT,PDType,CType}(parameter_distributions, constraints, name) + ParameterDistribution(parameter_distribution::PDType, + constraints::Array{CType}, + name::String) where {PDType <: ParameterDistributionType, + CType <: ConstraintType} = !(len(parameter_distribution) == +len(constraints)) ? throw(DimensionMismatch("There must be one constraint per parameter")) : new{PDType,CType}(parameter_distribution, constraints, name) end + """ struct ParameterDistributions @@ -122,7 +128,7 @@ end Structure to hold an array of ParameterDistribution's """ struct ParameterDistributions - parameter_distribution::Array{ParameterDistribution} + parameter_distributions::Array{ParameterDistribution} end @@ -131,12 +137,12 @@ end ## Functions """ - get_names(pds::ParameterDistributions) + get_name(pds::ParameterDistributions) Returns a list of ParameterDistribution names """ function get_name(pds::ParameterDistributions) - return [get_name(pd) for pd in pds] + return [get_name(pd) for pd in pds.parameter_distributions] end function get_name(pd::ParameterDistribution) @@ -149,16 +155,18 @@ end Returns a `Dict` of `ParameterDistribution` distributions by name, (unless sample type) """ function get_distribution(pds::ParameterDistributions) - return Dict{String,Any}(get_name(pd) => get_distribution(pd) for pd in pds) + return Dict{String,Any}(get_name(pd) => get_distribution(pd) for pd in pds.parameter_distributions) +end + +function get_distribution(pd::ParameterDistribution) + return get_distribution(pd.distribution) end -function get_distribution(pd::ParameterDistribution{FT,Samples,CType}) where {FT <: AbstractFloat, - CType <: ConstraintType} +function get_distribution(d::Samples) return "Contains samples only" end -function get_distribution(pd::ParameterDistribution{FT,Parameterized,CType}) where {FT <: AbstractFloat, - CType <: ConstraintType} - return pd.distribution.distribution +function get_distribution(d::Parameterized) + return d.distribution end """ @@ -171,7 +179,7 @@ function sample_distribution(pds::ParameterDistributions) end function sample_distribution(pds::ParameterDistributions, n_samples::IT) where {IT <: Integer} - return Dict{String,Any}(get_name(pd) => sample_distribution(pd, n_samples) for pd in pds) + return Dict{String,Any}(get_name(pd) => sample_distribution(pd, n_samples) for pd in pds.parameter_distributions) end function sample_distribution(pd::ParameterDistribution) @@ -179,126 +187,120 @@ function sample_distribution(pd::ParameterDistribution) end function sample_distribution(pd::ParameterDistribution, n_samples::IT) where {IT <: Integer} - return sample_distribution(pd.distribution) + return sample_distribution(pd.distribution,n_samples) end function sample_distribution(d::Samples, n_samples::IT) where {IT <: Integer} total_samples = size(d.distribution_samples)[1] - samples_idx = sample(collect(1:total_samples), n_samples ; replace=false) - return d.distribution_samples[samples_idx,:] + samples_idx = sample(collect(1:total_samples), n_samples; replace=false) + return d.distribution_samples[samples_idx, :] end -function sample_distribution(d::Parameterized,n_samples::IT) - return rand(pd.distirbution.distribution, n_samples) +function sample_distribution(d::Parameterized,n_samples::IT) where {IT <: Integer} + return rand(d.distribution, n_samples) end #apply transforms -function length(c::CType) where {CType <: ConstraintType} - return 1 -end -function length(carray::Array{CType}) where {CType <: ConstraintType} - return size(carray)[1] -end """ - transform_real_to_prior(pds::ParameterDistributions, x::Array{AbstractFloat}) + transform_real_to_prior(pds::ParameterDistributions, x::Array{Real}) Apply the transformation to map real (and possibly constrained) parameters `xarray` into the unbounded prior space """ -function transform_real_to_prior(pds::ParameterDistributions, xarray::Array{FT}) where {FT <: AbstractFloat} +function transform_real_to_prior(pds::ParameterDistributions, xarray::Array{FT}) where {FT <: Real} #split xarray into chunks to be fed into each distribution - clen = [length(pd.constraints) for pd in pds] # e.g [2,3,1,3] + clen = [len(pd.constraints) for pd in pds.parameter_distributions] # e.g [2,3,1,3] cumulative_clen = [sum(clen[1:i]) for i = 1:size(clen)[1]] # e.g [1 1:2, 3:5, 6, 7:9 ] - x_idx = Dict(Integer,Array{Integer})(i => collect(cumulative_cl[i-1]+1:cumulative_cl[i]) + x_idx = Dict(Integer,Array{Integer})(i => collect(cumulative_cl[i - 1] + 1:cumulative_cl[i]) for i in 2:size(cumulative_clen)[1]) x_idx[1] = collect(1:cumulative_clen[1]) - return Dict{String,Any}(get_name(pd) => transform_real_to_prior(pd,xarray[x_idx[i]]) for (i,pd) in enumerate(pds)) + return Dict{String,Any}(get_name(pd) => transform_real_to_prior(pd,xarray[x_idx[i]]) for (i,pd) in enumerate(pds.parameter_distributions)) end -function transform_real_to_prior(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: AbstractFloat} +function transform_real_to_prior(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real} return [transform_real_to_prior(constraint,xarray[i]) for (i,constraint) in enumerate(pd.constraints)] end """ No constraint mapping x -> x """ -function transform_real_to_prior(c::NoConstraints , x::Array{FT}) where {FT <: AbstractFloat} +function transform_real_to_prior(c::NoConstraint , x::Array{FT}) where {FT <: Real} return x end """ Bounded below -> unbounded, use mapping x -> ln(x - lower_bound) """ -function transform_real_to_prior(c::BoundedBelow, x::Array{FT}) where {FT <: AbstractFloat} +function transform_real_to_prior(c::BoundedBelow, x::Array{FT}) where {FT <: Real} return ln(x - c.lower_bound) end """ Bounded above -> unbounded, use mapping x -> ln(upper_bound - x) """ -function transform_real_to_prior(c::BoundedAbove, x::Array{FT}) where {FT <: AbstractFloat} +function transform_real_to_prior(c::BoundedAbove, x::Array{FT}) where {FT <: Real} return ln(c.upper_bound - x) end """ Bounded -> unbounded, use mapping x -> ln((x - lower_bound) / (upper_bound - x) """ -function transform_real_to_prior(c::Bounded, x::Array{FT}) where {FT <: AbstractFloat} +function transform_real_to_prior(c::Bounded, x::Array{FT}) where {FT <: Real} return ln((c.upper_bound - x) / (x - c.lower_bound)) end """ - transform_prior_to_real(pds::ParameterDistributions, xarray::Array{AbstractFloat}) + transform_prior_to_real(pds::ParameterDistributions, xarray::Array{Real}) Apply the transformation to map parameters `xarray` from the unbounded space into (possibly constrained) real space """ function transform_prior_to_real(pds::ParameterDistributions, - xarray::Array{FT}) where {FT <: AbstractFloat} + xarray::Array{FT}) where {FT <: Real} #split xarray into chunks to be fed into each distribution - clen = [length(pd.constraints) for pd in pds] # e.g [2,3,1,3] + clen = [len(pd.constraints) for pd in pds] # e.g [2,3,1,3] cumulative_clen = [sum(clen[1:i]) for i = 1:size(clen)[1]] # e.g [1 1:2, 3:5, 6, 7:9 ] x_idx = Dict(Integer,Array{Integer})(i => collect(cumulative_cl[i-1]+1:cumulative_cl[i]) for i in 2:size(cumulative_clen)[1]) x_idx[1] = collect(1:cumulative_clen[1]) return Dict{String,Any}(get_name(pd) => transform_prior_to_real(pd,xarray[x_idx[i]]) - for (i,pd) in enumerate(pds)) + for (i,pd) in enumerate(pds.parameter_distributions)) end -function transform_prior_to_real(pd::ParameterDistributions, xarray::Array{FT}) where {FT <: AbstractFloat} +function transform_prior_to_real(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real} return [transform_prior_to_real(constraint,xarray[i]) for (i,constraint) in enumerate(pd.constraints)] end """ No constraint mapping x -> x """ -function transform_prior_to_real(c::NoConstraints , x::FT) where {FT <: AbstractFloat} +function transform_prior_to_real(c::NoConstraint , x::FT) where {FT <: Real} return x end """ Unbounded -> bounded below, use mapping x -> exp(x) + lower_bound """ -function transform_prior_to_real(c::BoundedBelow, x::FT) where {FT <: AbstractFloat} +function transform_prior_to_real(c::BoundedBelow, x::FT) where {FT <: Real} return exp(x) + c.lower_bound end """ Unbounded -> bounded above, use mapping x -> upper_bound - exp(x) """ -function transform_prior_to_real(c::BoundedAbove, x::FT) where {FT <: AbstractFloat} +function transform_prior_to_real(c::BoundedAbove, x::FT) where {FT <: Real} return c.upper_bound - exp(x) end """ Unbounded -> bounded, use mapping x -> (upper_bound * exp(x) + lower_bound) / (exp(x) + 1) """ -function transform_prior_to_real(c::Bounded, x::FT) where {FT <: AbstractFloat} +function transform_prior_to_real(c::Bounded, x::FT) where {FT <: Real} return (c.upper_bound * exp(x) + c.lower_bound) / (exp(x) + 1) end diff --git a/test/ParameterDistributions/runtests.jl b/test/ParameterDistributions/runtests.jl new file mode 100644 index 000000000..3e389db79 --- /dev/null +++ b/test/ParameterDistributions/runtests.jl @@ -0,0 +1,132 @@ +using Test +using Distributions +using Random + +using CalibrateEmulateSample.ParameterDistributionsStorage + +@testset "ParameterDistributions" begin + @testset "ParameterDistributionType" begin + # Tests for the ParameterDistributionType + d = Parameterized(Gamma(2.0, 0.8)) + @test d.distribution == Gamma(2.0,0.8) + + d = Samples([[1 2 3]; [4 5 6]]) + @test d.distribution_samples == [[1.0 2.0 3.0]; [4.0 5.0 6.0]] + end + @testset "ConstraintType" begin + # Tests for the ConstraintType + c = BoundedBelow(0.2) + @test c.lower_bound == 0.2 + + c = BoundedAbove(0.2) + @test c.upper_bound == 0.2 + + c = Bounded(-0.1,0.2) + @test c.lower_bound == -0.1 + @test c.upper_bound == 0.2 + @test_throws DomainError Bounded(0.2,-0.1) + end + + @testset "ParameterDistribution(s)" begin + # Tests for the ParameterDistribution + d = Parameterized(MvNormal(4,0.1)) + c = [NoConstraint(), + BoundedBelow(-1.0), + BoundedAbove(0.4), + Bounded(-0.1,0.2)] + + name = "constrained_mvnormal" + u = ParameterDistribution(d,c,name) + @test u.distribution == d + @test u.constraints == c + @test u.name == name + @test_throws DimensionMismatch ParameterDistribution(d,c[1:3],name) + + # Tests for the ParameterDistribuions + d1 = Parameterized(MvNormal(4,0.1)) + c1 = [NoConstraint(), + BoundedBelow(-1.0), + BoundedAbove(0.4), + Bounded(-0.1,0.2)] + name1 = "constrained_mvnormal" + u1 = ParameterDistribution(d1,c1,name1) + + d2 = Samples([[1.0 3.0]; [5.0 7.0]; [9.0 11.0]; [13.0 15.0]]) + c2 = [Bounded(10,15), + NoConstraint()] + name2 = "constrained_sampled" + u2 = ParameterDistribution(d2,c2,name2) + + u = ParameterDistributions([u1,u2]) + @test u.parameter_distributions == [u1,u2] + end + + @testset "get/sample functions" begin + # setup for the tests: + d1 = Parameterized(MvNormal(4,0.1)) + c1 = [NoConstraint(), + BoundedBelow(-1.0), + BoundedAbove(0.4), + Bounded(-0.1,0.2)] + name1 = "constrained_mvnormal" + u1 = ParameterDistribution(d1,c1,name1) + + d2 = Samples([[1.0 3.0]; [5.0 7.0]; [9.0 11.0]; [13.0 15.0]]) + c2 = [Bounded(10,15), + NoConstraint()] + name2 = "constrained_sampled" + u2 = ParameterDistribution(d2,c2,name2) + + u = ParameterDistributions([u1,u2]) + + # Tests for get_name + @test get_name(u1) == name1 + @test get_name(u) == [name1, name2] + + # Tests for get_distribution + @test get_distribution(d1) == MvNormal(4,0.1) + @test get_distribution(u1) == MvNormal(4,0.1) + @test typeof(get_distribution(d2)) <: String + @test typeof(get_distribution(u2)) <: String + + d = get_distribution(u) + @test d[name1] == MvNormal(4,0.1) + @test typeof(d[name2]) <: String + + # Tests for sample distribution + Random.seed!(2020) + s1 = rand(MvNormal(4,0.1),1) + Random.seed!(2020) + @test sample_distribution(u1) == s1 + + Random.seed!(2020) + s1 = rand(MvNormal(4,0.1),3) + Random.seed!(2020) + @test sample_distribution(u1,3) == s1 + + Random.seed!(2020) + idx = sample(collect(1:size(d2.distribution_samples)[1]),1) + s2 = d2.distribution_samples[idx, :] + Random.seed!(2020) + @test sample_distribution(u2) == s2 + + Random.seed!(2020) + idx = sample(collect(1:size(d2.distribution_samples)[1]), 3; replace=false) + s2 = d2.distribution_samples[idx, :] + Random.seed!(2020) + @test sample_distribution(u2,3) == s2 + + Random.seed!(2020) + s1 = sample_distribution(u1,3) + s2 = sample_distribution(u2,3) + Random.seed!(2020) + s = sample_distribution(u,3) + @test s[name1] == s1 + @test s[name2] == s2 + end + + @testset "transform functions" begin + # Tests for transforms + end + +end diff --git a/test/runtests.jl b/test/runtests.jl index fbe348d7b..d1842c556 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -29,7 +29,8 @@ end "GPEmulator", "MCMC", "Observations", - "Utilities"] + "Utilities", + "ParameterDistributions"] if all_tests || has_submodule(submodule) || "CalibrateEmulateSample" in ARGS include_test(submodule) From 10a1614faefbeedd80eee92b11f8906a40f66b0e Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 25 Nov 2020 17:28:16 +0000 Subject: [PATCH 05/14] added StatsBase --- Project.toml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7218ea773..61583b419 100644 --- a/Project.toml +++ b/Project.toml @@ -11,12 +11,14 @@ GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [extras] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Test"] From efef03f5ca1aa64b15cffc4bb40e8c12e887c22b Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 25 Nov 2020 17:39:27 +0000 Subject: [PATCH 06/14] include StatsBase --- test/ParameterDistributions/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/ParameterDistributions/runtests.jl b/test/ParameterDistributions/runtests.jl index 3e389db79..cd1588eb6 100644 --- a/test/ParameterDistributions/runtests.jl +++ b/test/ParameterDistributions/runtests.jl @@ -1,5 +1,6 @@ using Test using Distributions +using StatsBase using Random using CalibrateEmulateSample.ParameterDistributionsStorage From 16e64068e923f550fffa059461102891c4c35295 Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 25 Nov 2020 18:12:08 +0000 Subject: [PATCH 07/14] ensured uniqueness of StatsBase sample function --- src/ParameterDistributions.jl | 2 +- test/ParameterDistributions/runtests.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ParameterDistributions.jl b/src/ParameterDistributions.jl index 34137356f..5838609fd 100644 --- a/src/ParameterDistributions.jl +++ b/src/ParameterDistributions.jl @@ -192,7 +192,7 @@ end function sample_distribution(d::Samples, n_samples::IT) where {IT <: Integer} total_samples = size(d.distribution_samples)[1] - samples_idx = sample(collect(1:total_samples), n_samples; replace=false) + samples_idx = StatsBase.sample(collect(1:total_samples), n_samples; replace=false) return d.distribution_samples[samples_idx, :] end diff --git a/test/ParameterDistributions/runtests.jl b/test/ParameterDistributions/runtests.jl index cd1588eb6..fab45719c 100644 --- a/test/ParameterDistributions/runtests.jl +++ b/test/ParameterDistributions/runtests.jl @@ -106,13 +106,13 @@ using CalibrateEmulateSample.ParameterDistributionsStorage @test sample_distribution(u1,3) == s1 Random.seed!(2020) - idx = sample(collect(1:size(d2.distribution_samples)[1]),1) + idx = StatsBase.sample(collect(1:size(d2.distribution_samples)[1]),1) s2 = d2.distribution_samples[idx, :] Random.seed!(2020) @test sample_distribution(u2) == s2 Random.seed!(2020) - idx = sample(collect(1:size(d2.distribution_samples)[1]), 3; replace=false) + idx = StatsBase.sample(collect(1:size(d2.distribution_samples)[1]), 3; replace=false) s2 = d2.distribution_samples[idx, :] Random.seed!(2020) @test sample_distribution(u2,3) == s2 From 04ffee5d094086b09aadb94f533b702b88d18106 Mon Sep 17 00:00:00 2001 From: odunbar Date: Thu, 26 Nov 2020 15:07:48 +0000 Subject: [PATCH 08/14] runtests for transforms --- src/ParameterDistributions.jl | 31 ++++++----- test/ParameterDistributions/runtests.jl | 72 ++++++++++++++++++++----- 2 files changed, 73 insertions(+), 30 deletions(-) diff --git a/src/ParameterDistributions.jl b/src/ParameterDistributions.jl index 5838609fd..e573b212a 100644 --- a/src/ParameterDistributions.jl +++ b/src/ParameterDistributions.jl @@ -212,11 +212,11 @@ function transform_real_to_prior(pds::ParameterDistributions, xarray::Array{FT}) #split xarray into chunks to be fed into each distribution clen = [len(pd.constraints) for pd in pds.parameter_distributions] # e.g [2,3,1,3] cumulative_clen = [sum(clen[1:i]) for i = 1:size(clen)[1]] # e.g [1 1:2, 3:5, 6, 7:9 ] - x_idx = Dict(Integer,Array{Integer})(i => collect(cumulative_cl[i - 1] + 1:cumulative_cl[i]) + x_idx = Dict{Integer,Array{Integer}}(i => collect(cumulative_clen[i - 1] + 1:cumulative_clen[i]) for i in 2:size(cumulative_clen)[1]) x_idx[1] = collect(1:cumulative_clen[1]) - return Dict{String,Any}(get_name(pd) => transform_real_to_prior(pd,xarray[x_idx[i]]) for (i,pd) in enumerate(pds.parameter_distributions)) + return cat([transform_real_to_prior(pd,xarray[x_idx[i]]) for (i,pd) in enumerate(pds.parameter_distributions)]...,dims=1) end function transform_real_to_prior(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real} @@ -226,29 +226,29 @@ end """ No constraint mapping x -> x """ -function transform_real_to_prior(c::NoConstraint , x::Array{FT}) where {FT <: Real} +function transform_real_to_prior(c::NoConstraint , x::FT) where {FT <: Real} return x end """ -Bounded below -> unbounded, use mapping x -> ln(x - lower_bound) +Bounded below -> unbounded, use mapping x -> log(x - lower_bound) """ -function transform_real_to_prior(c::BoundedBelow, x::Array{FT}) where {FT <: Real} - return ln(x - c.lower_bound) +function transform_real_to_prior(c::BoundedBelow, x::FT) where {FT <: Real} + return log(x - c.lower_bound) end """ -Bounded above -> unbounded, use mapping x -> ln(upper_bound - x) +Bounded above -> unbounded, use mapping x -> log(upper_bound - x) """ -function transform_real_to_prior(c::BoundedAbove, x::Array{FT}) where {FT <: Real} - return ln(c.upper_bound - x) +function transform_real_to_prior(c::BoundedAbove, x::FT) where {FT <: Real} + return log(c.upper_bound - x) end """ -Bounded -> unbounded, use mapping x -> ln((x - lower_bound) / (upper_bound - x) +Bounded -> unbounded, use mapping x -> log((x - lower_bound) / (upper_bound - x) """ -function transform_real_to_prior(c::Bounded, x::Array{FT}) where {FT <: Real} - return ln((c.upper_bound - x) / (x - c.lower_bound)) +function transform_real_to_prior(c::Bounded, x::FT) where {FT <: Real} + return log( (x - c.lower_bound) / (c.upper_bound - x)) end @@ -262,14 +262,13 @@ function transform_prior_to_real(pds::ParameterDistributions, xarray::Array{FT}) where {FT <: Real} #split xarray into chunks to be fed into each distribution - clen = [len(pd.constraints) for pd in pds] # e.g [2,3,1,3] + clen = [len(pd.constraints) for pd in pds.parameter_distributions] # e.g [2,3,1,3] cumulative_clen = [sum(clen[1:i]) for i = 1:size(clen)[1]] # e.g [1 1:2, 3:5, 6, 7:9 ] - x_idx = Dict(Integer,Array{Integer})(i => collect(cumulative_cl[i-1]+1:cumulative_cl[i]) + x_idx = Dict{Integer,Array{Integer}}(i => collect(cumulative_clen[i-1]+1:cumulative_clen[i]) for i in 2:size(cumulative_clen)[1]) x_idx[1] = collect(1:cumulative_clen[1]) - return Dict{String,Any}(get_name(pd) => transform_prior_to_real(pd,xarray[x_idx[i]]) - for (i,pd) in enumerate(pds.parameter_distributions)) + return cat([transform_prior_to_real(pd,xarray[x_idx[i]]) for (i,pd) in enumerate(pds.parameter_distributions)]...,dims=1) end function transform_prior_to_real(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real} diff --git a/test/ParameterDistributions/runtests.jl b/test/ParameterDistributions/runtests.jl index fab45719c..649369550 100644 --- a/test/ParameterDistributions/runtests.jl +++ b/test/ParameterDistributions/runtests.jl @@ -11,8 +11,8 @@ using CalibrateEmulateSample.ParameterDistributionsStorage d = Parameterized(Gamma(2.0, 0.8)) @test d.distribution == Gamma(2.0,0.8) - d = Samples([[1 2 3]; [4 5 6]]) - @test d.distribution_samples == [[1.0 2.0 3.0]; [4.0 5.0 6.0]] + d = Samples([1 2 3; 4 5 6]) + @test d.distribution_samples == [1.0 2.0 3.0; 4.0 5.0 6.0] end @testset "ConstraintType" begin # Tests for the ConstraintType @@ -52,7 +52,7 @@ using CalibrateEmulateSample.ParameterDistributionsStorage name1 = "constrained_mvnormal" u1 = ParameterDistribution(d1,c1,name1) - d2 = Samples([[1.0 3.0]; [5.0 7.0]; [9.0 11.0]; [13.0 15.0]]) + d2 = Samples([1.0 3.0; 5.0 7.0; 9.0 11.0; 13.0 15.0]) c2 = [Bounded(10,15), NoConstraint()] name2 = "constrained_sampled" @@ -72,7 +72,7 @@ using CalibrateEmulateSample.ParameterDistributionsStorage name1 = "constrained_mvnormal" u1 = ParameterDistribution(d1,c1,name1) - d2 = Samples([[1.0 3.0]; [5.0 7.0]; [9.0 11.0]; [13.0 15.0]]) + d2 = Samples([1.0 3.0; 5.0 7.0; 9.0 11.0; 13.0 15.0]) c2 = [Bounded(10,15), NoConstraint()] name2 = "constrained_sampled" @@ -95,39 +95,83 @@ using CalibrateEmulateSample.ParameterDistributionsStorage @test typeof(d[name2]) <: String # Tests for sample distribution - Random.seed!(2020) + seed=2020 + Random.seed!(seed) s1 = rand(MvNormal(4,0.1),1) - Random.seed!(2020) + Random.seed!(seed) @test sample_distribution(u1) == s1 - Random.seed!(2020) + Random.seed!(seed) s1 = rand(MvNormal(4,0.1),3) - Random.seed!(2020) + Random.seed!(seed) @test sample_distribution(u1,3) == s1 - Random.seed!(2020) + Random.seed!(seed) idx = StatsBase.sample(collect(1:size(d2.distribution_samples)[1]),1) s2 = d2.distribution_samples[idx, :] - Random.seed!(2020) + Random.seed!(seed) @test sample_distribution(u2) == s2 - Random.seed!(2020) + Random.seed!(seed) idx = StatsBase.sample(collect(1:size(d2.distribution_samples)[1]), 3; replace=false) s2 = d2.distribution_samples[idx, :] - Random.seed!(2020) + Random.seed!(seed) @test sample_distribution(u2,3) == s2 - Random.seed!(2020) + Random.seed!(seed) s1 = sample_distribution(u1,3) s2 = sample_distribution(u2,3) - Random.seed!(2020) + Random.seed!(seed) s = sample_distribution(u,3) @test s[name1] == s1 @test s[name2] == s2 end @testset "transform functions" begin + #setup for the tests + d1 = Parameterized(MvNormal(4,0.1)) + c1 = [NoConstraint(), + BoundedBelow(-1.0), + BoundedAbove(0.4), + Bounded(-0.1,0.2)] + name1 = "constrained_mvnormal" + u1 = ParameterDistribution(d1,c1,name1) + + d2 = Samples([1.0 3.0; 5.0 7.0; 9.0 11.0; 13.0 15.0]) + c2 = [Bounded(10,15), + NoConstraint()] + name2 = "constrained_sampled" + u2 = ParameterDistribution(d2,c2,name2) + + u = ParameterDistributions([u1,u2]) + + x_unbd = rand(MvNormal(6,3), 1000) #6 x 1000 # Tests for transforms + # prior to real + x_real_noconstraint = map(x -> transform_prior_to_real(NoConstraint(), x), x_unbd[1,:]) + @test x_real_noconstraint == x_unbd[1,:] + x_real_below = map(x -> transform_prior_to_real(BoundedBelow(30), x), x_unbd[1,:]) + @test all(x_real_below .> 30) + x_real_above = map(x -> transform_prior_to_real(BoundedAbove(-10), x), x_unbd[1,:]) + @test all(x_real_above .< -10) + x_real_bounded = map(x -> transform_prior_to_real(Bounded(-2,-1), x), x_unbd[1,:]) + @test all([all(x_real_bounded .> -2), all(x_real_bounded .< -1)]) + + # prior to real + @test isapprox(x_unbd[1,:] - map(x -> transform_real_to_prior(NoConstraint(), x), x_real_noconstraint), zeros(size(x_unbd)[2]); atol=1e-6) + @test isapprox(x_unbd[1,:] - map(x -> transform_real_to_prior(BoundedBelow(30), x), x_real_below), zeros(size(x_unbd)[2]) ; atol=1e-6) + @test isapprox(x_unbd[1,:] - map(x -> transform_real_to_prior(BoundedAbove(-10), x), x_real_above), zeros(size(x_unbd)[2]); atol=1e-6) + @test isapprox(x_unbd[1,:] - map(x -> transform_real_to_prior(Bounded(-2,-1), x), x_real_bounded), zeros(size(x_unbd)[2]); atol=1e-6) + + x_real_constrained1 = mapslices(x -> transform_prior_to_real(u1,x), x_unbd[1:4,:]; dims=1) + @test isapprox(x_unbd[1:4,:] - mapslices(x -> transform_real_to_prior(u1,x), x_real_constrained1; dims=1), zeros(size(x_unbd[1:4,:])) ; atol=1e-6) + x_real_constrained2 = mapslices(x -> transform_prior_to_real(u2,x), x_unbd[5:6,:]; dims=1) + @test isapprox(x_unbd[5:6,:] - mapslices(x -> transform_real_to_prior(u2,x), x_real_constrained2; dims=1), zeros(size(x_unbd[5:6,:])); atol=1e-6) + + x_real = mapslices(x -> transform_prior_to_real(u,x), x_unbd; dims=1) + x_unbd_tmp = mapslices(x -> transform_real_to_prior(u,x), x_real; dims=1) + @test isapprox(x_unbd - x_unbd_tmp,zeros(size(x_unbd));atol=1e-6) + end end From 8775299e1f31948cf7f4034e578bf5cf51702023 Mon Sep 17 00:00:00 2001 From: odunbar Date: Thu, 26 Nov 2020 22:39:43 +0000 Subject: [PATCH 09/14] remove Revise from Project toml --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 61583b419..9d5914d6d 100644 --- a/Project.toml +++ b/Project.toml @@ -11,7 +11,6 @@ GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" From 5e324de8c94914031969ca2128ea5919a0c550e9 Mon Sep 17 00:00:00 2001 From: odunbar Date: Fri, 27 Nov 2020 00:28:51 +0000 Subject: [PATCH 10/14] more definition on storage of parameters in columns --- src/ParameterDistributions.jl | 51 ++++++++++++++++--------- test/ParameterDistributions/runtests.jl | 29 ++++++++------ 2 files changed, 50 insertions(+), 30 deletions(-) diff --git a/src/ParameterDistributions.jl b/src/ParameterDistributions.jl index e573b212a..be6b209d2 100644 --- a/src/ParameterDistributions.jl +++ b/src/ParameterDistributions.jl @@ -23,7 +23,7 @@ abstract type ParameterDistributionType end """ Parameterized <: ParameterDistributionType - + A distribution constructed from a parametrized formula (e.g Julia Distributions.jl) """ struct Parameterized <: ParameterDistributionType @@ -33,10 +33,14 @@ end """ Samples{FT<:Real} <: ParameterDistributionType -A distribution comprised of only samples +A distribution comprised of only samples, stored as columns of parameters """ struct Samples{FT<:Real} <: ParameterDistributionType - distribution_samples::Array{FT} + distribution_samples::Array{FT,2} #parameters are columns + Samples(distribution_samples::Array{FT,2}; params_are_columns=true) where {FT <: Real} = params_are_columns ? new{FT}(distribution_samples) : new{FT}(permutedims(distribution_samples,[2,1])) + #Distinguish 1 sample of an ND parameter or N samples of 1D parameter, and store as 2D array + Samples(distribution_samples::Array{FT,1}; params_are_columns=true) where {FT <: Real} = params_are_columns ? new{FT}(reshape(distribution_samples,1,:)) : new{FT}(reshape(distribution_samples,:,1)) + end @@ -91,16 +95,26 @@ end function len(carray::Array{CType}) where {CType <: ConstraintType} return size(carray)[1] end + """ - len(d::ParametrizedDistributionType) + dimension(d::ParametrizedDistributionType) The number of dimensions of the parameter space """ -function len(d::Parameterized) +function dimension(d::Parameterized) return length(d.distribution) end -function len(d::Samples) +function dimension(d::Samples) + return size(d.distribution_samples)[1] +end + +""" + n_samples(d::Samples) + +The number of samples in the array +""" +function n_samples(d::Samples) return size(d.distribution_samples)[2] end @@ -117,7 +131,7 @@ struct ParameterDistribution{PDType <: ParameterDistributionType, CType <: Const ParameterDistribution(parameter_distribution::PDType, constraints::Array{CType}, name::String) where {PDType <: ParameterDistributionType, - CType <: ConstraintType} = !(len(parameter_distribution) == + CType <: ConstraintType} = !(dimension(parameter_distribution) == len(constraints)) ? throw(DimensionMismatch("There must be one constraint per parameter")) : new{PDType,CType}(parameter_distribution, constraints, name) end @@ -132,8 +146,6 @@ struct ParameterDistributions end - - ## Functions """ @@ -178,26 +190,27 @@ function sample_distribution(pds::ParameterDistributions) return sample_distribution(pds,1) end -function sample_distribution(pds::ParameterDistributions, n_samples::IT) where {IT <: Integer} - return Dict{String,Any}(get_name(pd) => sample_distribution(pd, n_samples) for pd in pds.parameter_distributions) +function sample_distribution(pds::ParameterDistributions, n_draws::IT) where {IT <: Integer} + return cat([sample_distribution(pd,n_draws) for pd in pds.parameter_distributions]...,dims=1) end function sample_distribution(pd::ParameterDistribution) return sample_distribution(pd,1) end -function sample_distribution(pd::ParameterDistribution, n_samples::IT) where {IT <: Integer} - return sample_distribution(pd.distribution,n_samples) +function sample_distribution(pd::ParameterDistribution, n_draws::IT) where {IT <: Integer} + return sample_distribution(pd.distribution,n_draws) end -function sample_distribution(d::Samples, n_samples::IT) where {IT <: Integer} - total_samples = size(d.distribution_samples)[1] - samples_idx = StatsBase.sample(collect(1:total_samples), n_samples; replace=false) - return d.distribution_samples[samples_idx, :] +function sample_distribution(d::Samples, n_draws::IT) where {IT <: Integer} + n_stored_samples = n_samples(d) + samples_idx = StatsBase.sample(collect(1:n_stored_samples), n_draws) + return d.distribution_samples[:,samples_idx] + end -function sample_distribution(d::Parameterized,n_samples::IT) where {IT <: Integer} - return rand(d.distribution, n_samples) +function sample_distribution(d::Parameterized, n_draws::IT) where {IT <: Integer} + return rand(d.distribution, n_draws) end diff --git a/test/ParameterDistributions/runtests.jl b/test/ParameterDistributions/runtests.jl index 649369550..af20b2698 100644 --- a/test/ParameterDistributions/runtests.jl +++ b/test/ParameterDistributions/runtests.jl @@ -10,9 +10,16 @@ using CalibrateEmulateSample.ParameterDistributionsStorage # Tests for the ParameterDistributionType d = Parameterized(Gamma(2.0, 0.8)) @test d.distribution == Gamma(2.0,0.8) - + d = Samples([1 2 3; 4 5 6]) - @test d.distribution_samples == [1.0 2.0 3.0; 4.0 5.0 6.0] + @test d.distribution_samples == [1 2 3; 4 5 6] + d = Samples([1 4; 2 5; 3 6]; params_are_columns=false) + @test d.distribution_samples == [1 2 3; 4 5 6] + d = Samples([1, 2, 3, 4, 5, 6]) + @test d.distribution_samples == [1 2 3 4 5 6] + d = Samples([1, 2, 3, 4, 5, 6]; params_are_columns=false) + @test d.distribution_samples == reshape([1, 2, 3, 4, 5, 6],:,1) + end @testset "ConstraintType" begin # Tests for the ConstraintType @@ -52,7 +59,7 @@ using CalibrateEmulateSample.ParameterDistributionsStorage name1 = "constrained_mvnormal" u1 = ParameterDistribution(d1,c1,name1) - d2 = Samples([1.0 3.0; 5.0 7.0; 9.0 11.0; 13.0 15.0]) + d2 = Samples([1.0 3.0 5.0 7.0; 9.0 11.0 13.0 15.0]) c2 = [Bounded(10,15), NoConstraint()] name2 = "constrained_sampled" @@ -72,7 +79,8 @@ using CalibrateEmulateSample.ParameterDistributionsStorage name1 = "constrained_mvnormal" u1 = ParameterDistribution(d1,c1,name1) - d2 = Samples([1.0 3.0; 5.0 7.0; 9.0 11.0; 13.0 15.0]) + d2 = Samples([1 2 3 4 + 5 6 7 8]) c2 = [Bounded(10,15), NoConstraint()] name2 = "constrained_sampled" @@ -107,14 +115,14 @@ using CalibrateEmulateSample.ParameterDistributionsStorage @test sample_distribution(u1,3) == s1 Random.seed!(seed) - idx = StatsBase.sample(collect(1:size(d2.distribution_samples)[1]),1) - s2 = d2.distribution_samples[idx, :] + idx = StatsBase.sample(collect(1:size(d2.distribution_samples)[2]),1) + s2 = d2.distribution_samples[:, idx] Random.seed!(seed) @test sample_distribution(u2) == s2 Random.seed!(seed) - idx = StatsBase.sample(collect(1:size(d2.distribution_samples)[1]), 3; replace=false) - s2 = d2.distribution_samples[idx, :] + idx = StatsBase.sample(collect(1:size(d2.distribution_samples)[2]), 3) + s2 = d2.distribution_samples[:, idx] Random.seed!(seed) @test sample_distribution(u2,3) == s2 @@ -123,8 +131,7 @@ using CalibrateEmulateSample.ParameterDistributionsStorage s2 = sample_distribution(u2,3) Random.seed!(seed) s = sample_distribution(u,3) - @test s[name1] == s1 - @test s[name2] == s2 + @test s == cat([s1,s2]...,dims=1) end @testset "transform functions" begin @@ -137,7 +144,7 @@ using CalibrateEmulateSample.ParameterDistributionsStorage name1 = "constrained_mvnormal" u1 = ParameterDistribution(d1,c1,name1) - d2 = Samples([1.0 3.0; 5.0 7.0; 9.0 11.0; 13.0 15.0]) + d2 = Samples([1.0 3.0 5.0 7.0; 9.0 11.0 13.0 15.0]) c2 = [Bounded(10,15), NoConstraint()] name2 = "constrained_sampled" From f1c4bfef9e03fbf5bd42c7fb7c59129247219d52 Mon Sep 17 00:00:00 2001 From: odunbar Date: Sat, 28 Nov 2020 02:57:38 +0000 Subject: [PATCH 11/14] Removed ParameterDistributions and renamed files --- src/CalibrateEmulateSample.jl | 2 +- ...tributions.jl => ParameterDistribution.jl} | 115 +++++++----------- .../runtests.jl | 45 ++++--- test/runtests.jl | 2 +- 4 files changed, 72 insertions(+), 92 deletions(-) rename src/{ParameterDistributions.jl => ParameterDistribution.jl} (63%) rename test/{ParameterDistributions => ParameterDistribution}/runtests.jl (86%) diff --git a/src/CalibrateEmulateSample.jl b/src/CalibrateEmulateSample.jl index 0a8720c70..247b6b9e8 100644 --- a/src/CalibrateEmulateSample.jl +++ b/src/CalibrateEmulateSample.jl @@ -5,7 +5,7 @@ using Distributions, Statistics, LinearAlgebra, DocStringExtensions # No internal deps, light external deps include("Observations.jl") include("Priors.jl") -include("ParameterDistributions.jl") +include("ParameterDistribution.jl") # No internal deps, heavy external deps include("EKP.jl") diff --git a/src/ParameterDistributions.jl b/src/ParameterDistribution.jl similarity index 63% rename from src/ParameterDistributions.jl rename to src/ParameterDistribution.jl index be6b209d2..44dffe0f2 100644 --- a/src/ParameterDistributions.jl +++ b/src/ParameterDistribution.jl @@ -1,4 +1,4 @@ -module ParameterDistributionsStorage +module ParameterDistributionStorage ## Imports using Distributions @@ -10,7 +10,7 @@ using Random #objects export Parameterized, Samples export NoConstraint, BoundedBelow, BoundedAbove, Bounded -export ParameterDistribution, ParameterDistributions +export ParameterDistribution #functions export get_name, get_distribution @@ -121,57 +121,59 @@ end """ struct ParameterDistribution -Structure to hold a parameter distribution +Structure to hold a parameter distribution, always stored as an array of distributions """ -struct ParameterDistribution{PDType <: ParameterDistributionType, CType <: ConstraintType} - distribution::PDType +struct ParameterDistribution{PDType <: ParameterDistributionType, CType <: ConstraintType, ST <: AbstractString} + distributions::Array{PDType} constraints::Array{CType} - name::String - - ParameterDistribution(parameter_distribution::PDType, - constraints::Array{CType}, - name::String) where {PDType <: ParameterDistributionType, - CType <: ConstraintType} = !(dimension(parameter_distribution) == -len(constraints)) ? throw(DimensionMismatch("There must be one constraint per parameter")) : new{PDType,CType}(parameter_distribution, constraints, name) -end + names::Array{ST} + + function ParameterDistribution(parameter_distributions::Union{PDType,Array{PDType}}, + constraints::Union{CType, Array{CType}, Array}, + names::Union{ST, Array{ST}}) where {PDType <: ParameterDistributionType, + CType <: ConstraintType, + ST <: AbstractString} + + parameter_distributions = isa(parameter_distributions, PDType) ? [parameter_distributions] : parameter_distributions + constraints = isa(constraints, Union{<:ConstraintType,Array{<:ConstraintType}}) ? [constraints] : constraints #to calc n_constraints_per_dist + names = isa(names, ST) ? [names] : names + + n_parameter_per_dist = [dimension(pd) for pd in parameter_distributions] + n_constraints_per_dist = [len(c) for c in constraints] + n_dists = length(parameter_distributions) + n_names = length(names) + if !(n_parameter_per_dist == n_constraints_per_dist) + throw(DimensionMismatch("There must be one constraint per parameter in a distribution, use NoConstraint() type if no constraint is required")) + elseif !(n_dists == n_names) + throw(DimensionMismatch("There must be one name per parameter distribution")) + else + constraints = cat(constraints...,dims=1) + + new{PDType,ConstraintType,ST}(parameter_distributions, constraints, names) + end + end - -""" - struct ParameterDistributions - -Structure to hold an array of ParameterDistribution's -""" -struct ParameterDistributions - parameter_distributions::Array{ParameterDistribution} end - + ## Functions """ - get_name(pds::ParameterDistributions) + get_name(pd::ParameterDistribution) Returns a list of ParameterDistribution names """ -function get_name(pds::ParameterDistributions) - return [get_name(pd) for pd in pds.parameter_distributions] -end - function get_name(pd::ParameterDistribution) - return pd.name + return pd.names end """ - get_distributions(pds::ParameterDistributions) + get_distribution(pd::ParameterDistribution) Returns a `Dict` of `ParameterDistribution` distributions by name, (unless sample type) """ -function get_distribution(pds::ParameterDistributions) - return Dict{String,Any}(get_name(pd) => get_distribution(pd) for pd in pds.parameter_distributions) -end - function get_distribution(pd::ParameterDistribution) - return get_distribution(pd.distribution) + return Dict{String,Any}(pd.names[i] => get_distribution(d) for (i,d) in enumerate(pd.distributions)) end function get_distribution(d::Samples) @@ -182,24 +184,16 @@ function get_distribution(d::Parameterized) end """ - function sample_distribution(pds::ParameterDistributions) + function sample_distribution(pd::ParameterDistribution) Draws samples from the parameter distributions """ -function sample_distribution(pds::ParameterDistributions) - return sample_distribution(pds,1) -end - -function sample_distribution(pds::ParameterDistributions, n_draws::IT) where {IT <: Integer} - return cat([sample_distribution(pd,n_draws) for pd in pds.parameter_distributions]...,dims=1) -end - function sample_distribution(pd::ParameterDistribution) return sample_distribution(pd,1) end function sample_distribution(pd::ParameterDistribution, n_draws::IT) where {IT <: Integer} - return sample_distribution(pd.distribution,n_draws) + return cat([sample_distribution(d,n_draws) for d in pd.distributions]...,dims=1) end function sample_distribution(d::Samples, n_draws::IT) where {IT <: Integer} @@ -217,23 +211,13 @@ end #apply transforms """ - transform_real_to_prior(pds::ParameterDistributions, x::Array{Real}) + transform_real_to_prior(pd::ParameterDistribution, x::Array{Real}) Apply the transformation to map real (and possibly constrained) parameters `xarray` into the unbounded prior space """ -function transform_real_to_prior(pds::ParameterDistributions, xarray::Array{FT}) where {FT <: Real} - #split xarray into chunks to be fed into each distribution - clen = [len(pd.constraints) for pd in pds.parameter_distributions] # e.g [2,3,1,3] - cumulative_clen = [sum(clen[1:i]) for i = 1:size(clen)[1]] # e.g [1 1:2, 3:5, 6, 7:9 ] - x_idx = Dict{Integer,Array{Integer}}(i => collect(cumulative_clen[i - 1] + 1:cumulative_clen[i]) - for i in 2:size(cumulative_clen)[1]) - x_idx[1] = collect(1:cumulative_clen[1]) - - return cat([transform_real_to_prior(pd,xarray[x_idx[i]]) for (i,pd) in enumerate(pds.parameter_distributions)]...,dims=1) -end - function transform_real_to_prior(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real} - return [transform_real_to_prior(constraint,xarray[i]) for (i,constraint) in enumerate(pd.constraints)] + #split xarray into chunks + return cat([transform_real_to_prior(c,xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1) end """ @@ -267,25 +251,12 @@ end """ - transform_prior_to_real(pds::ParameterDistributions, xarray::Array{Real}) + transform_prior_to_real(pd::ParameterDistribution, xarray::Array{Real}) Apply the transformation to map parameters `xarray` from the unbounded space into (possibly constrained) real space """ -function transform_prior_to_real(pds::ParameterDistributions, - xarray::Array{FT}) where {FT <: Real} - - #split xarray into chunks to be fed into each distribution - clen = [len(pd.constraints) for pd in pds.parameter_distributions] # e.g [2,3,1,3] - cumulative_clen = [sum(clen[1:i]) for i = 1:size(clen)[1]] # e.g [1 1:2, 3:5, 6, 7:9 ] - x_idx = Dict{Integer,Array{Integer}}(i => collect(cumulative_clen[i-1]+1:cumulative_clen[i]) - for i in 2:size(cumulative_clen)[1]) - x_idx[1] = collect(1:cumulative_clen[1]) - - return cat([transform_prior_to_real(pd,xarray[x_idx[i]]) for (i,pd) in enumerate(pds.parameter_distributions)]...,dims=1) -end - function transform_prior_to_real(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real} - return [transform_prior_to_real(constraint,xarray[i]) for (i,constraint) in enumerate(pd.constraints)] + return [transform_prior_to_real(c,xarray[i]) for (i,c) in enumerate(pd.constraints)] end """ diff --git a/test/ParameterDistributions/runtests.jl b/test/ParameterDistribution/runtests.jl similarity index 86% rename from test/ParameterDistributions/runtests.jl rename to test/ParameterDistribution/runtests.jl index af20b2698..3041fca7f 100644 --- a/test/ParameterDistributions/runtests.jl +++ b/test/ParameterDistribution/runtests.jl @@ -3,9 +3,9 @@ using Distributions using StatsBase using Random -using CalibrateEmulateSample.ParameterDistributionsStorage +using CalibrateEmulateSample.ParameterDistributionStorage -@testset "ParameterDistributions" begin +@testset "ParameterDistribution" begin @testset "ParameterDistributionType" begin # Tests for the ParameterDistributionType d = Parameterized(Gamma(2.0, 0.8)) @@ -36,6 +36,15 @@ using CalibrateEmulateSample.ParameterDistributionsStorage end @testset "ParameterDistribution(s)" begin + + d = Parameterized(Normal(0,1)) + c = NoConstraint() + name = "unconstrained_normal" + u = ParameterDistribution(d,c,name) + @test u.distributions == [d] + @test u.constraints == [c] + @test u.names == [name] + # Tests for the ParameterDistribution d = Parameterized(MvNormal(4,0.1)) c = [NoConstraint(), @@ -45,12 +54,12 @@ using CalibrateEmulateSample.ParameterDistributionsStorage name = "constrained_mvnormal" u = ParameterDistribution(d,c,name) - @test u.distribution == d + @test u.distributions == [d] @test u.constraints == c - @test u.name == name + @test u.names == [name] @test_throws DimensionMismatch ParameterDistribution(d,c[1:3],name) - - # Tests for the ParameterDistribuions + + # Tests for the ParameterDistribution d1 = Parameterized(MvNormal(4,0.1)) c1 = [NoConstraint(), BoundedBelow(-1.0), @@ -65,8 +74,10 @@ using CalibrateEmulateSample.ParameterDistributionsStorage name2 = "constrained_sampled" u2 = ParameterDistribution(d2,c2,name2) - u = ParameterDistributions([u1,u2]) - @test u.parameter_distributions == [u1,u2] + u = ParameterDistribution([d1,d2], [c1,c2], [name1,name2]) + @test u.distributions == [d1,d2] + @test u.constraints == cat([c1,c2]...,dims=1) + @test u.names == [name1,name2] end @testset "get/sample functions" begin @@ -79,24 +90,22 @@ using CalibrateEmulateSample.ParameterDistributionsStorage name1 = "constrained_mvnormal" u1 = ParameterDistribution(d1,c1,name1) - d2 = Samples([1 2 3 4 - 5 6 7 8]) - c2 = [Bounded(10,15), - NoConstraint()] + d2 = Samples([1 2 3 4]) + c2 = [Bounded(10,15)] name2 = "constrained_sampled" u2 = ParameterDistribution(d2,c2,name2) - u = ParameterDistributions([u1,u2]) - + u = ParameterDistribution([d1,d2], [c1,c2], [name1,name2]) + # Tests for get_name - @test get_name(u1) == name1 + @test get_name(u1) == [name1] @test get_name(u) == [name1, name2] # Tests for get_distribution @test get_distribution(d1) == MvNormal(4,0.1) - @test get_distribution(u1) == MvNormal(4,0.1) + @test get_distribution(u1)[name1] == MvNormal(4,0.1) @test typeof(get_distribution(d2)) <: String - @test typeof(get_distribution(u2)) <: String + @test typeof(get_distribution(u2)[name2]) <: String d = get_distribution(u) @test d[name1] == MvNormal(4,0.1) @@ -150,7 +159,7 @@ using CalibrateEmulateSample.ParameterDistributionsStorage name2 = "constrained_sampled" u2 = ParameterDistribution(d2,c2,name2) - u = ParameterDistributions([u1,u2]) + u = ParameterDistribution([d1,d2], [c1,c2], [name1,name2]) x_unbd = rand(MvNormal(6,3), 1000) #6 x 1000 # Tests for transforms diff --git a/test/runtests.jl b/test/runtests.jl index d1842c556..ed3c4f384 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -30,7 +30,7 @@ end "MCMC", "Observations", "Utilities", - "ParameterDistributions"] + "ParameterDistribution"] if all_tests || has_submodule(submodule) || "CalibrateEmulateSample" in ARGS include_test(submodule) From 9a6d7e56721c121e56fe138df1d8b58a790896ae Mon Sep 17 00:00:00 2001 From: odunbar Date: Sat, 28 Nov 2020 02:59:36 +0000 Subject: [PATCH 12/14] renamed test --- test/ParameterDistribution/runtests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/ParameterDistribution/runtests.jl b/test/ParameterDistribution/runtests.jl index 3041fca7f..7a94efd7f 100644 --- a/test/ParameterDistribution/runtests.jl +++ b/test/ParameterDistribution/runtests.jl @@ -5,7 +5,7 @@ using Random using CalibrateEmulateSample.ParameterDistributionStorage -@testset "ParameterDistribution" begin +@testset "ParameterDistributionStorage" begin @testset "ParameterDistributionType" begin # Tests for the ParameterDistributionType d = Parameterized(Gamma(2.0, 0.8)) @@ -35,7 +35,7 @@ using CalibrateEmulateSample.ParameterDistributionStorage @test_throws DomainError Bounded(0.2,-0.1) end - @testset "ParameterDistribution(s)" begin + @testset "ParameterDistribution" begin d = Parameterized(Normal(0,1)) c = NoConstraint() From 029ff80bb981a7158a053013bee8efe31476579b Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 2 Dec 2020 20:52:38 +0000 Subject: [PATCH 13/14] changed transform names --- src/ParameterDistribution.jl | 30 +++++++++++++------------- test/ParameterDistribution/runtests.jl | 28 ++++++++++++------------ 2 files changed, 29 insertions(+), 29 deletions(-) diff --git a/src/ParameterDistribution.jl b/src/ParameterDistribution.jl index 44dffe0f2..629623f09 100644 --- a/src/ParameterDistribution.jl +++ b/src/ParameterDistribution.jl @@ -15,7 +15,7 @@ export ParameterDistribution #functions export get_name, get_distribution export sample_distribution -export transform_real_to_prior, transform_prior_to_real +export transform_constrained_to_unconstrained, transform_unconstrained_to_constrained ## Objects # for the Distribution @@ -211,79 +211,79 @@ end #apply transforms """ - transform_real_to_prior(pd::ParameterDistribution, x::Array{Real}) + transform_constrained_to_unconstrained(pd::ParameterDistribution, x::Array{Real}) Apply the transformation to map real (and possibly constrained) parameters `xarray` into the unbounded prior space """ -function transform_real_to_prior(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real} +function transform_constrained_to_unconstrained(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real} #split xarray into chunks - return cat([transform_real_to_prior(c,xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1) + return cat([transform_constrained_to_unconstrained(c,xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1) end """ No constraint mapping x -> x """ -function transform_real_to_prior(c::NoConstraint , x::FT) where {FT <: Real} +function transform_constrained_to_unconstrained(c::NoConstraint , x::FT) where {FT <: Real} return x end """ Bounded below -> unbounded, use mapping x -> log(x - lower_bound) """ -function transform_real_to_prior(c::BoundedBelow, x::FT) where {FT <: Real} +function transform_constrained_to_unconstrained(c::BoundedBelow, x::FT) where {FT <: Real} return log(x - c.lower_bound) end """ Bounded above -> unbounded, use mapping x -> log(upper_bound - x) """ -function transform_real_to_prior(c::BoundedAbove, x::FT) where {FT <: Real} +function transform_constrained_to_unconstrained(c::BoundedAbove, x::FT) where {FT <: Real} return log(c.upper_bound - x) end """ Bounded -> unbounded, use mapping x -> log((x - lower_bound) / (upper_bound - x) """ -function transform_real_to_prior(c::Bounded, x::FT) where {FT <: Real} +function transform_constrained_to_unconstrained(c::Bounded, x::FT) where {FT <: Real} return log( (x - c.lower_bound) / (c.upper_bound - x)) end """ - transform_prior_to_real(pd::ParameterDistribution, xarray::Array{Real}) + transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{Real}) Apply the transformation to map parameters `xarray` from the unbounded space into (possibly constrained) real space """ -function transform_prior_to_real(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real} - return [transform_prior_to_real(c,xarray[i]) for (i,c) in enumerate(pd.constraints)] +function transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real} + return [transform_unconstrained_to_constrained(c,xarray[i]) for (i,c) in enumerate(pd.constraints)] end """ No constraint mapping x -> x """ -function transform_prior_to_real(c::NoConstraint , x::FT) where {FT <: Real} +function transform_unconstrained_to_constrained(c::NoConstraint , x::FT) where {FT <: Real} return x end """ Unbounded -> bounded below, use mapping x -> exp(x) + lower_bound """ -function transform_prior_to_real(c::BoundedBelow, x::FT) where {FT <: Real} +function transform_unconstrained_to_constrained(c::BoundedBelow, x::FT) where {FT <: Real} return exp(x) + c.lower_bound end """ Unbounded -> bounded above, use mapping x -> upper_bound - exp(x) """ -function transform_prior_to_real(c::BoundedAbove, x::FT) where {FT <: Real} +function transform_unconstrained_to_constrained(c::BoundedAbove, x::FT) where {FT <: Real} return c.upper_bound - exp(x) end """ Unbounded -> bounded, use mapping x -> (upper_bound * exp(x) + lower_bound) / (exp(x) + 1) """ -function transform_prior_to_real(c::Bounded, x::FT) where {FT <: Real} +function transform_unconstrained_to_constrained(c::Bounded, x::FT) where {FT <: Real} return (c.upper_bound * exp(x) + c.lower_bound) / (exp(x) + 1) end diff --git a/test/ParameterDistribution/runtests.jl b/test/ParameterDistribution/runtests.jl index 7a94efd7f..7185bc3a3 100644 --- a/test/ParameterDistribution/runtests.jl +++ b/test/ParameterDistribution/runtests.jl @@ -164,28 +164,28 @@ using CalibrateEmulateSample.ParameterDistributionStorage x_unbd = rand(MvNormal(6,3), 1000) #6 x 1000 # Tests for transforms # prior to real - x_real_noconstraint = map(x -> transform_prior_to_real(NoConstraint(), x), x_unbd[1,:]) + x_real_noconstraint = map(x -> transform_unconstrained_to_constrained(NoConstraint(), x), x_unbd[1,:]) @test x_real_noconstraint == x_unbd[1,:] - x_real_below = map(x -> transform_prior_to_real(BoundedBelow(30), x), x_unbd[1,:]) + x_real_below = map(x -> transform_unconstrained_to_constrained(BoundedBelow(30), x), x_unbd[1,:]) @test all(x_real_below .> 30) - x_real_above = map(x -> transform_prior_to_real(BoundedAbove(-10), x), x_unbd[1,:]) + x_real_above = map(x -> transform_unconstrained_to_constrained(BoundedAbove(-10), x), x_unbd[1,:]) @test all(x_real_above .< -10) - x_real_bounded = map(x -> transform_prior_to_real(Bounded(-2,-1), x), x_unbd[1,:]) + x_real_bounded = map(x -> transform_unconstrained_to_constrained(Bounded(-2,-1), x), x_unbd[1,:]) @test all([all(x_real_bounded .> -2), all(x_real_bounded .< -1)]) # prior to real - @test isapprox(x_unbd[1,:] - map(x -> transform_real_to_prior(NoConstraint(), x), x_real_noconstraint), zeros(size(x_unbd)[2]); atol=1e-6) - @test isapprox(x_unbd[1,:] - map(x -> transform_real_to_prior(BoundedBelow(30), x), x_real_below), zeros(size(x_unbd)[2]) ; atol=1e-6) - @test isapprox(x_unbd[1,:] - map(x -> transform_real_to_prior(BoundedAbove(-10), x), x_real_above), zeros(size(x_unbd)[2]); atol=1e-6) - @test isapprox(x_unbd[1,:] - map(x -> transform_real_to_prior(Bounded(-2,-1), x), x_real_bounded), zeros(size(x_unbd)[2]); atol=1e-6) + @test isapprox(x_unbd[1,:] - map(x -> transform_constrained_to_unconstrained(NoConstraint(), x), x_real_noconstraint), zeros(size(x_unbd)[2]); atol=1e-6) + @test isapprox(x_unbd[1,:] - map(x -> transform_constrained_to_unconstrained(BoundedBelow(30), x), x_real_below), zeros(size(x_unbd)[2]) ; atol=1e-6) + @test isapprox(x_unbd[1,:] - map(x -> transform_constrained_to_unconstrained(BoundedAbove(-10), x), x_real_above), zeros(size(x_unbd)[2]); atol=1e-6) + @test isapprox(x_unbd[1,:] - map(x -> transform_constrained_to_unconstrained(Bounded(-2,-1), x), x_real_bounded), zeros(size(x_unbd)[2]); atol=1e-6) - x_real_constrained1 = mapslices(x -> transform_prior_to_real(u1,x), x_unbd[1:4,:]; dims=1) - @test isapprox(x_unbd[1:4,:] - mapslices(x -> transform_real_to_prior(u1,x), x_real_constrained1; dims=1), zeros(size(x_unbd[1:4,:])) ; atol=1e-6) - x_real_constrained2 = mapslices(x -> transform_prior_to_real(u2,x), x_unbd[5:6,:]; dims=1) - @test isapprox(x_unbd[5:6,:] - mapslices(x -> transform_real_to_prior(u2,x), x_real_constrained2; dims=1), zeros(size(x_unbd[5:6,:])); atol=1e-6) + x_real_constrained1 = mapslices(x -> transform_unconstrained_to_constrained(u1,x), x_unbd[1:4,:]; dims=1) + @test isapprox(x_unbd[1:4,:] - mapslices(x -> transform_constrained_to_unconstrained(u1,x), x_real_constrained1; dims=1), zeros(size(x_unbd[1:4,:])) ; atol=1e-6) + x_real_constrained2 = mapslices(x -> transform_unconstrained_to_constrained(u2,x), x_unbd[5:6,:]; dims=1) + @test isapprox(x_unbd[5:6,:] - mapslices(x -> transform_constrained_to_unconstrained(u2,x), x_real_constrained2; dims=1), zeros(size(x_unbd[5:6,:])); atol=1e-6) - x_real = mapslices(x -> transform_prior_to_real(u,x), x_unbd; dims=1) - x_unbd_tmp = mapslices(x -> transform_real_to_prior(u,x), x_real; dims=1) + x_real = mapslices(x -> transform_unconstrained_to_constrained(u,x), x_unbd; dims=1) + x_unbd_tmp = mapslices(x -> transform_constrained_to_unconstrained(u,x), x_real; dims=1) @test isapprox(x_unbd - x_unbd_tmp,zeros(size(x_unbd));atol=1e-6) end From cc07d5e67c7c3b9e39ef499a60c2c4a74dd3ab99 Mon Sep 17 00:00:00 2001 From: odunbar Date: Thu, 3 Dec 2020 12:58:20 +0000 Subject: [PATCH 14/14] Constraints store only the transform function now, not the bounds --- src/ParameterDistribution.jl | 189 +++++++++++++++---------- test/ParameterDistribution/runtests.jl | 98 ++++++------- 2 files changed, 161 insertions(+), 126 deletions(-) diff --git a/src/ParameterDistribution.jl b/src/ParameterDistribution.jl index 629623f09..a9466d2cd 100644 --- a/src/ParameterDistribution.jl +++ b/src/ParameterDistribution.jl @@ -9,12 +9,13 @@ using Random #objects export Parameterized, Samples -export NoConstraint, BoundedBelow, BoundedAbove, Bounded export ParameterDistribution +export Constraint #functions export get_name, get_distribution export sample_distribution +export no_constraint, bounded_below, bounded_above, bounded export transform_constrained_to_unconstrained, transform_unconstrained_to_constrained ## Objects @@ -46,41 +47,68 @@ end # For the transforms abstract type ConstraintType end - """ - NoConstraint <: ConstraintType + Constraint <: ConstraintType -No constraint. +Contains two functions to map between constrained and unconstrained spaces. """ -struct NoConstraint <: ConstraintType end +struct Constraint <: ConstraintType + constrained_to_unconstrained::Function + unconstrained_to_constrained::Function +end + """ - BoundedBelow{FT <: Real} <: ConstraintType + no_constraint() -A lower bound constraint. +Constructs a Constraint with no constraints, enforced by maps x -> x and x -> x. """ -struct BoundedBelow{FT <: Real} <: ConstraintType - lower_bound::FT +function no_constraint() + c_to_u = (x -> x) + u_to_c = (x -> x) + return Constraint(c_to_u,u_to_c) end + """ - BoundedAbove{FT <: Real} <: ConstraintType + bounded_below(lower_bound::FT) where {FT <: Real} -And upper bound constraint +Constructs a Constraint with provided lower bound, enforced by maps x -> log(x - lower_bound) and x -> exp(x) + lower_bound. """ -struct BoundedAbove{FT <: Real} <: ConstraintType - upper_bound::FT +function bounded_below(lower_bound::FT) where {FT <: Real} + c_to_u = ( x -> log(x - lower_bound) ) + u_to_c = ( x -> exp(x) + lower_bound ) + return Constraint(c_to_u, u_to_c) end +""" + bounded_above(upper_bound::FT) where {FT <: Real} + +Constructs a Constraint with provided upper bound, enforced by maps x -> log(upper_bound - x) and x -> upper_bound - exp(x). +""" +function bounded_above(upper_bound::FT) where {FT<:Real} + c_to_u = ( x -> log(upper_bound - x) ) + u_to_c = ( x -> upper_bound - exp(x) ) + return Constraint(c_to_u, u_to_c) +end + + """ Bounded{FT <: Real} <: ConstraintType -Both a lower and upper bound constraint. +Constructs a Constraint with provided upper and lower bounds, enforced by maps +x -> log((x - lower_bound) / (upper_bound - x)) +and +x -> (upper_bound * exp(x) + lower_bound) / (exp(x) + 1) + """ -struct Bounded{FT <: Real} <: ConstraintType - lower_bound::FT - upper_bound::FT - Bounded(lower_bound::FT,upper_bound::FT) where {FT <: Real} = upper_bound <= lower_bound ? throw(DomainError("upper bound must be greater than lower bound")) : new{FT}(lower_bound, upper_bound) +function bounded(lower_bound::FT, upper_bound::FT) where {FT <: Real} + if (upper_bound <= lower_bound) + throw(DomainError("upper bound must be greater than lower bound")) + end + c_to_u = ( x -> log( (x - lower_bound) / (upper_bound - x)) ) + u_to_c = ( x -> (upper_bound * exp(x) + lower_bound) / (exp(x) + 1) ) + return Constraint(c_to_u, u_to_c) end """ @@ -97,7 +125,7 @@ function len(carray::Array{CType}) where {CType <: ConstraintType} end """ - dimension(d::ParametrizedDistributionType) + dimension(d<:ParametrizedDistributionType) The number of dimensions of the parameter space """ @@ -156,6 +184,7 @@ struct ParameterDistribution{PDType <: ParameterDistributionType, CType <: Const end + ## Functions """ @@ -213,79 +242,85 @@ end """ transform_constrained_to_unconstrained(pd::ParameterDistribution, x::Array{Real}) -Apply the transformation to map real (and possibly constrained) parameters `xarray` into the unbounded prior space +Apply the transformation to map (possibly constrained) parameters `xarray` into the unconstrained space """ function transform_constrained_to_unconstrained(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real} #split xarray into chunks - return cat([transform_constrained_to_unconstrained(c,xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1) -end - -""" -No constraint mapping x -> x -""" -function transform_constrained_to_unconstrained(c::NoConstraint , x::FT) where {FT <: Real} - return x +# return cat([transform_constrained_to_unconstrained(c,xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1) + return cat([c.constrained_to_unconstrained(xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1) end -""" -Bounded below -> unbounded, use mapping x -> log(x - lower_bound) -""" -function transform_constrained_to_unconstrained(c::BoundedBelow, x::FT) where {FT <: Real} - return log(x - c.lower_bound) -end - -""" -Bounded above -> unbounded, use mapping x -> log(upper_bound - x) -""" -function transform_constrained_to_unconstrained(c::BoundedAbove, x::FT) where {FT <: Real} - return log(c.upper_bound - x) -end - -""" -Bounded -> unbounded, use mapping x -> log((x - lower_bound) / (upper_bound - x) -""" -function transform_constrained_to_unconstrained(c::Bounded, x::FT) where {FT <: Real} - return log( (x - c.lower_bound) / (c.upper_bound - x)) -end - - - """ transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{Real}) -Apply the transformation to map parameters `xarray` from the unbounded space into (possibly constrained) real space +Apply the transformation to map parameters `xarray` from the unconstrained space into (possibly constrained) space """ function transform_unconstrained_to_constrained(pd::ParameterDistribution, xarray::Array{FT}) where {FT <: Real} - return [transform_unconstrained_to_constrained(c,xarray[i]) for (i,c) in enumerate(pd.constraints)] +# return [transform_unconstrained_to_constrained(c,xarray[i]) for (i,c) in enumerate(pd.constraints)] + return cat([c.unconstrained_to_constrained(xarray[i]) for (i,c) in enumerate(pd.constraints)]...,dims=1) end -""" -No constraint mapping x -> x -""" -function transform_unconstrained_to_constrained(c::NoConstraint , x::FT) where {FT <: Real} - return x -end -""" -Unbounded -> bounded below, use mapping x -> exp(x) + lower_bound -""" -function transform_unconstrained_to_constrained(c::BoundedBelow, x::FT) where {FT <: Real} - return exp(x) + c.lower_bound -end -""" -Unbounded -> bounded above, use mapping x -> upper_bound - exp(x) -""" -function transform_unconstrained_to_constrained(c::BoundedAbove, x::FT) where {FT <: Real} - return c.upper_bound - exp(x) -end -""" -Unbounded -> bounded, use mapping x -> (upper_bound * exp(x) + lower_bound) / (exp(x) + 1) -""" -function transform_unconstrained_to_constrained(c::Bounded, x::FT) where {FT <: Real} - return (c.upper_bound * exp(x) + c.lower_bound) / (exp(x) + 1) -end +# """ +# No constraint mapping x -> x +# """ +# function transform_constrained_to_unconstrained(c::NoConstraint , x::FT) where {FT <: Real} +# return x +# end + +# """ +# Bounded below -> unbounded, use mapping x -> log(x - lower_bound) +# """ +# function transform_constrained_to_unconstrained(c::BoundedBelow, x::FT) where {FT <: Real} +# return log(x - c.lower_bound) +# end + +# """ +# Bounded above -> unbounded, use mapping x -> log(upper_bound - x) +# """ +# function transform_constrained_to_unconstrained(c::BoundedAbove, x::FT) where {FT <: Real} +# return log(c.upper_bound - x) +# end + +# """ +# Bounded -> unbounded, use mapping x -> log((x - lower_bound) / (upper_bound - x) +# """ +# function transform_constrained_to_unconstrained(c::Bounded, x::FT) where {FT <: Real} +# return log( (x - c.lower_bound) / (c.upper_bound - x)) +# end + + + + +# """ +# No constraint mapping x -> x +# """ +# function transform_unconstrained_to_constrained(c::NoConstraint , x::FT) where {FT <: Real} +# return x +# end + +# """ +# Unbounded -> bounded below, use mapping x -> exp(x) + lower_bound +# """ +# function transform_unconstrained_to_constrained(c::BoundedBelow, x::FT) where {FT <: Real} +# return exp(x) + c.lower_bound +# end + +# """ +# Unbounded -> bounded above, use mapping x -> upper_bound - exp(x) +# """ +# function transform_unconstrained_to_constrained(c::BoundedAbove, x::FT) where {FT <: Real} +# return c.upper_bound - exp(x) +# end + +# """ +# Unbounded -> bounded, use mapping x -> (upper_bound * exp(x) + lower_bound) / (exp(x) + 1) +# """ +# function transform_unconstrained_to_constrained(c::Bounded, x::FT) where {FT <: Real} +# return (c.upper_bound * exp(x) + c.lower_bound) / (exp(x) + 1) +# end diff --git a/test/ParameterDistribution/runtests.jl b/test/ParameterDistribution/runtests.jl index 7185bc3a3..f5fa154b9 100644 --- a/test/ParameterDistribution/runtests.jl +++ b/test/ParameterDistribution/runtests.jl @@ -23,22 +23,36 @@ using CalibrateEmulateSample.ParameterDistributionStorage end @testset "ConstraintType" begin # Tests for the ConstraintType - c = BoundedBelow(0.2) - @test c.lower_bound == 0.2 + # The predefined transforms + c1 = bounded_below(0.2) + @test isapprox(c1.constrained_to_unconstrained(1.0) - (log(1.0 - 0.2)), 0) + @test isapprox(c1.unconstrained_to_constrained(0.0) - (exp(0.0) + 0.2), 0) + + c2 = bounded_above(0.2) + @test isapprox(c2.constrained_to_unconstrained(-1.0) - (log(0.2 - -1.0)), 0) + @test isapprox(c2.unconstrained_to_constrained(10.0) - (0.2 - exp(10.0)), 0) + + + c3 = bounded(-0.1,0.2) + @test isapprox(c3.constrained_to_unconstrained(0.0) - (log( (0.0 - -0.1) / (0.2 - 0.0))) , 0 ) + @test isapprox(c3.unconstrained_to_constrained(1.0) - ((0.2 * exp(1.0) + -0.1) / (exp(1.0) + 1)) , 0) + @test_throws DomainError bounded(0.2,-0.1) + + #an example with user defined invertible transforms + c_to_u = (x -> 3 * x + 14) + u_to_c = (x -> (x - 14) / 3) - c = BoundedAbove(0.2) - @test c.upper_bound == 0.2 + c4 = Constraint(c_to_u, u_to_c) + @test isapprox( c4.constrained_to_unconstrained(5.0) - c_to_u(5.0) , 0) + @test isapprox( c4.unconstrained_to_constrained(5.0) - u_to_c(5.0) , 0) - c = Bounded(-0.1,0.2) - @test c.lower_bound == -0.1 - @test c.upper_bound == 0.2 - @test_throws DomainError Bounded(0.2,-0.1) + end @testset "ParameterDistribution" begin d = Parameterized(Normal(0,1)) - c = NoConstraint() + c = no_constraint() name = "unconstrained_normal" u = ParameterDistribution(d,c,name) @test u.distributions == [d] @@ -47,10 +61,10 @@ using CalibrateEmulateSample.ParameterDistributionStorage # Tests for the ParameterDistribution d = Parameterized(MvNormal(4,0.1)) - c = [NoConstraint(), - BoundedBelow(-1.0), - BoundedAbove(0.4), - Bounded(-0.1,0.2)] + c = [no_constraint(), + bounded_below(-1.0), + bounded_above(0.4), + bounded(-0.1,0.2)] name = "constrained_mvnormal" u = ParameterDistribution(d,c,name) @@ -58,19 +72,20 @@ using CalibrateEmulateSample.ParameterDistributionStorage @test u.constraints == c @test u.names == [name] @test_throws DimensionMismatch ParameterDistribution(d,c[1:3],name) + @test_throws DimensionMismatch ParameterDistribution(d,c, [name,"extra_name"]) # Tests for the ParameterDistribution d1 = Parameterized(MvNormal(4,0.1)) - c1 = [NoConstraint(), - BoundedBelow(-1.0), - BoundedAbove(0.4), - Bounded(-0.1,0.2)] + c1 = [no_constraint(), + bounded_below(-1.0), + bounded_above(0.4), + bounded(-0.1,0.2)] name1 = "constrained_mvnormal" u1 = ParameterDistribution(d1,c1,name1) d2 = Samples([1.0 3.0 5.0 7.0; 9.0 11.0 13.0 15.0]) - c2 = [Bounded(10,15), - NoConstraint()] + c2 = [bounded(10,15), + no_constraint()] name2 = "constrained_sampled" u2 = ParameterDistribution(d2,c2,name2) @@ -78,20 +93,21 @@ using CalibrateEmulateSample.ParameterDistributionStorage @test u.distributions == [d1,d2] @test u.constraints == cat([c1,c2]...,dims=1) @test u.names == [name1,name2] + end @testset "get/sample functions" begin # setup for the tests: d1 = Parameterized(MvNormal(4,0.1)) - c1 = [NoConstraint(), - BoundedBelow(-1.0), - BoundedAbove(0.4), - Bounded(-0.1,0.2)] + c1 = [no_constraint(), + bounded_below(-1.0), + bounded_above(0.4), + bounded(-0.1,0.2)] name1 = "constrained_mvnormal" u1 = ParameterDistribution(d1,c1,name1) d2 = Samples([1 2 3 4]) - c2 = [Bounded(10,15)] + c2 = [bounded(10,15)] name2 = "constrained_sampled" u2 = ParameterDistribution(d2,c2,name2) @@ -146,16 +162,16 @@ using CalibrateEmulateSample.ParameterDistributionStorage @testset "transform functions" begin #setup for the tests d1 = Parameterized(MvNormal(4,0.1)) - c1 = [NoConstraint(), - BoundedBelow(-1.0), - BoundedAbove(0.4), - Bounded(-0.1,0.2)] + c1 = [no_constraint(), + bounded_below(-1.0), + bounded_above(0.4), + bounded(-0.1,0.2)] name1 = "constrained_mvnormal" u1 = ParameterDistribution(d1,c1,name1) d2 = Samples([1.0 3.0 5.0 7.0; 9.0 11.0 13.0 15.0]) - c2 = [Bounded(10,15), - NoConstraint()] + c2 = [bounded(10,15), + no_constraint()] name2 = "constrained_sampled" u2 = ParameterDistribution(d2,c2,name2) @@ -163,30 +179,14 @@ using CalibrateEmulateSample.ParameterDistributionStorage x_unbd = rand(MvNormal(6,3), 1000) #6 x 1000 # Tests for transforms - # prior to real - x_real_noconstraint = map(x -> transform_unconstrained_to_constrained(NoConstraint(), x), x_unbd[1,:]) - @test x_real_noconstraint == x_unbd[1,:] - x_real_below = map(x -> transform_unconstrained_to_constrained(BoundedBelow(30), x), x_unbd[1,:]) - @test all(x_real_below .> 30) - x_real_above = map(x -> transform_unconstrained_to_constrained(BoundedAbove(-10), x), x_unbd[1,:]) - @test all(x_real_above .< -10) - x_real_bounded = map(x -> transform_unconstrained_to_constrained(Bounded(-2,-1), x), x_unbd[1,:]) - @test all([all(x_real_bounded .> -2), all(x_real_bounded .< -1)]) - - # prior to real - @test isapprox(x_unbd[1,:] - map(x -> transform_constrained_to_unconstrained(NoConstraint(), x), x_real_noconstraint), zeros(size(x_unbd)[2]); atol=1e-6) - @test isapprox(x_unbd[1,:] - map(x -> transform_constrained_to_unconstrained(BoundedBelow(30), x), x_real_below), zeros(size(x_unbd)[2]) ; atol=1e-6) - @test isapprox(x_unbd[1,:] - map(x -> transform_constrained_to_unconstrained(BoundedAbove(-10), x), x_real_above), zeros(size(x_unbd)[2]); atol=1e-6) - @test isapprox(x_unbd[1,:] - map(x -> transform_constrained_to_unconstrained(Bounded(-2,-1), x), x_real_bounded), zeros(size(x_unbd)[2]); atol=1e-6) - x_real_constrained1 = mapslices(x -> transform_unconstrained_to_constrained(u1,x), x_unbd[1:4,:]; dims=1) - @test isapprox(x_unbd[1:4,:] - mapslices(x -> transform_constrained_to_unconstrained(u1,x), x_real_constrained1; dims=1), zeros(size(x_unbd[1:4,:])) ; atol=1e-6) + @test isapprox(x_unbd[1:4,:] - mapslices(x -> transform_constrained_to_unconstrained(u1,x), x_real_constrained1; dims=1), zeros(size(x_unbd[1:4,:])); atol=1e-8) x_real_constrained2 = mapslices(x -> transform_unconstrained_to_constrained(u2,x), x_unbd[5:6,:]; dims=1) - @test isapprox(x_unbd[5:6,:] - mapslices(x -> transform_constrained_to_unconstrained(u2,x), x_real_constrained2; dims=1), zeros(size(x_unbd[5:6,:])); atol=1e-6) + @test isapprox(x_unbd[5:6,:] - mapslices(x -> transform_constrained_to_unconstrained(u2,x), x_real_constrained2; dims=1), zeros(size(x_unbd[5:6,:])); atol=1e-8) x_real = mapslices(x -> transform_unconstrained_to_constrained(u,x), x_unbd; dims=1) x_unbd_tmp = mapslices(x -> transform_constrained_to_unconstrained(u,x), x_real; dims=1) - @test isapprox(x_unbd - x_unbd_tmp,zeros(size(x_unbd));atol=1e-6) + @test isapprox(x_unbd - x_unbd_tmp,zeros(size(x_unbd)); atol=1e-8) end