From b8ad6cacaa9fd7d4ffd567ab88417209c4f9768c Mon Sep 17 00:00:00 2001 From: Weiqi Date: Mon, 2 Dec 2019 13:44:38 -0600 Subject: [PATCH 01/10] adding measure evaluation methods and tests --- src/InfiniteOpt.jl | 7 ++ src/MeasureEvalMethod/MeasureEvalMethod.jl | 13 ++ src/MeasureEvalMethod/methods.jl | 116 ++++++++++++++++++ src/measures.jl | 134 +++++++++++++++++++++ test/MeasureEvalMethod/methods.jl | 81 +++++++++++++ test/runtests.jl | 7 ++ 6 files changed, 358 insertions(+) create mode 100644 src/MeasureEvalMethod/MeasureEvalMethod.jl create mode 100644 src/MeasureEvalMethod/methods.jl create mode 100644 test/MeasureEvalMethod/methods.jl diff --git a/src/InfiniteOpt.jl b/src/InfiniteOpt.jl index 039ae5665..2f58fde78 100644 --- a/src/InfiniteOpt.jl +++ b/src/InfiniteOpt.jl @@ -29,6 +29,9 @@ include("utilities.jl") include("TranscriptionOpt/TranscriptionOpt.jl") using .TranscriptionOpt +include("MeasureEvalMethod/MeasureEvalMethod.jl") +using .MeasureEvalMethod + # Export model object datatype export InfiniteModel @@ -81,6 +84,10 @@ Measure, MeasureRef export add_measure, measure, measure_function, measure_data, expand, expand_all_measures! +# Export measure evaluation functions +export generate_measure_data, MC_sampling, Gauss_Legendre, Gauss_Hermite, +Gauss_Laguerre, infinite_transform + # Export transcription datatypes export TranscriptionData, TranscriptionModel diff --git a/src/MeasureEvalMethod/MeasureEvalMethod.jl b/src/MeasureEvalMethod/MeasureEvalMethod.jl new file mode 100644 index 000000000..065efc347 --- /dev/null +++ b/src/MeasureEvalMethod/MeasureEvalMethod.jl @@ -0,0 +1,13 @@ +module MeasureEvalMethod + +# Import the necessary packages. +import FastGaussQuadrature +using ..InfiniteOpt + +# include jl files here +include("methods.jl") +# Export functions here +export generate_measure_data, MC_sampling, Gauss_Legendre, Gauss_Hermite, +Gauss_Laguerre, infinite_transform + +end # end module diff --git a/src/MeasureEvalMethod/methods.jl b/src/MeasureEvalMethod/methods.jl new file mode 100644 index 000000000..5a9a34d55 --- /dev/null +++ b/src/MeasureEvalMethod/methods.jl @@ -0,0 +1,116 @@ +function generate_measure_data(param::InfiniteOpt.ParameterRef, + lb::Float64, ub::Float64, num_supports::Int; + method::Function = MC_sampling, name::String = "", + weight_func::Function = InfiniteOpt._w + )::InfiniteOpt.DiscreteMeasureData + (supports, coeffs) = method(lb, ub, num_supports) + return InfiniteOpt.DiscreteMeasureData(param, coeffs, supports, name, weight_func) +end + +function generate_measure_data(param::Vector{InfiniteOpt.ParameterRef}, + lb::Vector{Float64}, ub::Vector{Float64}, + num_supports::Int; + method::Function = MC_sampling, name::String = "", + weight_func::Function = InfiniteOpt._w + )::InfiniteOpt.MultiDiscreteMeasureData + (supports, coeffs) = method(lb, ub, num_supports) + return InfiniteOpt.MultiDiscreteMeasureData(param, coeffs, supports, name, weight_func) +end + +# MC sampling from uniform distribution over the interval [lb, ub] +function MC_sampling(lb::Float64, ub::Float64, num_supports::Int)::Tuple + if lb == -Inf || ub == Inf + (samples, _) = infinite_transform(lb, ub, num_supports) + else + samples = rand(num_supports) .* (ub - lb) .+ lb + end + return (samples, ones(num_supports) / num_supports * (ub - lb)) +end + +# MC sampling - multi-dim version +function MC_sampling(lb::Vector{Float64}, ub::Vector{Float64}, + num_supports::Int)::Tuple + samples = zeros(length(lb), num_supports) + for i in eachindex(lb) + (samples[i, :], _) = MC_sampling(lb[i], ub[i], num_supports) + end + return (samples, ones(num_supports) / num_supports * prod(ub .- lb)) +end + +# default Gaussian quadrature method for bounded domain +function Gauss_Legendre(lb::Float64, ub::Float64, num_supports::Int)::Tuple + (supports, coeffs) = FastGaussQuadrature.gausslegendre(num_supports) + supports = (ub - lb) / 2 * supports .+ (ub + lb) / 2 + coeffs = (ub - lb) / 2 * coeffs + return (supports, coeffs) +end + +# default Gaussian quadrature method for infinite domain +function Gauss_Hermite(lb::Float64, ub::Float64, num_supports::Int)::Tuple + if lb != -Inf || ub != Inf + error("Lower/upper bound is not infinity. Use other measure evaluation " * + "methods.") + end + (supports, coeffs) = FastGaussQuadrature.gausshermite(num_supports) + coeffs = coeffs .* exp.(supports.^2) + return (supports, coeffs) +end + +# default Gaussian quadrature method for semi-infinite domain +function Gauss_Laguerre(lb::Float64, ub::Float64, num_supports::Int)::Tuple + if ub == Inf + if lb == -Inf + error("The range is infinite. Use other measure evaluation methods.") + end + (supports, coeffs) = FastGaussQuadrature.gausslaguerre(num_supports) + coeffs = coeffs .* exp.(supports) + supports = supports .+ lb + return (supports, coeffs) + elseif lb == -Inf + (supports, coeffs) = FastGaussQuadrature.gausslaguerre(num_supports) + coeffs = coeffs .* exp.(supports) + supports = -supports .+ ub + return (supports, coeffs) + else + error("The range is bounded. Use other measure evaluation methods.") + end +end + +# transform (semi-)infinite domain to finite domain +function infinite_transform(lb::Float64, ub::Float64, num_supports::Int; + sub_method::Function = MC_sampling, + transform_x::Function = _default_x, + transform_dx::Function = _default_dx, + t_lb::Float64 = -convert(Float64, lb == -Inf && ub == Inf), + t_ub::Float64 = 1)::Tuple + if lb != -Inf || ub != Inf + error("The range is not (semi-)infinite. Use evaluation methods for " * + "bounded domains.") + end + (t_supports, t_coeffs) = sub_method(t_lb, t_ub, num_supports) + supports = transform_x.(t_supports, lb, ub) + coeffs = t_coeffs .* transform_dx.(t_supports, lb, ub) + return (supports, coeffs) +end + +function _default_x(t::Float64, lb::Float64, ub::Float64)::Float64 + if lb > -Inf + return lb + t / (1 - t) + elseif ub < Inf + return ub - (1 - t) / t + else + return t / (1 - t^2) + end +end + +function _default_dx(t::Float64, lb::Float64, ub::Float64)::Float64 + if lb > -Inf + return 1 / (1 - t)^2 + elseif ub < Inf + return 1 / t^2 + else + return (1 + t^2) / (1 - t^2)^2 + end +end + +# TODO: consider truncated distribution diff --git a/src/measures.jl b/src/measures.jl index 9397e3cb0..02cd1874b 100644 --- a/src/measures.jl +++ b/src/measures.jl @@ -1,3 +1,6 @@ +const MC = :MC +const GaussLegendre = :GaussLegendre + # Extend Base.copy for new variable types Base.copy(v::MeasureRef, new_model::InfiniteModel) = MeasureRef(new_model, v.index) @@ -408,6 +411,137 @@ function measure(expr::JuMP.AbstractJuMPScalar, return add_measure(model, meas) end +""" + measure(expr::JuMP.AbstractJuMPScalar, + params::Union{ParameterRef, Vector{ParameterRef}}, + lb::Union{Float64, Vector{Float64}}, + ub::Union{Float64, Vector{Float64}}; + eval_method::Function, num_supports::Int, weight_func::Function, + use_existing_supports::Bool = false)::MeasureRef + +Returns a measure reference that evaluates `expr` without using an object of +[`AbstractMeasureData`](@ref) type. Similar to the main [`measure`](@ref) +method, this function aims to implement measures of the form: +``\\int_{p \\in P} expr(p) w(p) dp`` where ``p`` is an infinite parameter (scalar +or vector) and ``w`` is the weight function. This function will serve as a +flexible interface where users only have to provide necessary data about the +integration. Instead of taking an [`AbstractMeasureData`](@ref) object as input, +this function constructs the [`AbstractMeasureData`](@ref) object using some +default numerical integration schemes. + +**Example** +```julia + +``` +""" +# Measure function that takes non-AbstractMeasureData types +function measure(expr::JuMP.AbstractJuMPScalar, + params::Union{ParameterRef, Vector{ParameterRef}, Nothing} = nothing, + lb::Union{Float64, Vector{Float64}} = Float64[], + ub::Union{Float64, Vector{Float64}} = Float64[]; + eval_method::Function = MC_sampling, num_supports::Int = 50, + weight_func::Function = _w, + use_existing_supports::Bool = false)::MeasureRef + + # Default: try to collect all parameters in expr. + if isa(params, Nothing) + params = _all_parameter_refs(expr) + if length(params) == 0 + error("No infinite parameters in the expression.") + end + end + + # Check if the parameters are empty + if isa(params, Vector) + if length(params) == 0 + error("Parameters cannot be empty.") + end + num_params = length(params) + else + num_params = 1 + end + + # Check if the parameters belong to multiple groups + ids = unique(group_id.(params)) + if length(ids) > 1 + error("Multiple groups of parameters in the expression. Need to " * + "specify the parameters to integrate over.") + end + + # Fill in lower bounds and upper bounds if not given + if length(lb) == 0 + params_have_lower_bounds = all(JuMP.has_lower_bound.(params)) + if !params_have_lower_bounds + error("Some parameter(s) do not have lower bounds. Need to manually " * + "input the lower bound values.") + end + lb = JuMP.lower_bound.(params) + end + if length(ub) == 0 + params_have_upper_bounds = all(JuMP.has_upper_bound.(params)) + if !params_have_upper_bounds + error("Some parameter(s) do not have upper bounds. Need to manually " * + "input the upper bound values.") + end + ub = JuMP.upper_bound.(params) + end + + # Check the dimension of lb and ub matches number of parameters + if length(lb) != num_params || length(ub) != num_params + error("Number of parameters do not match number of lower bounds or " * + "upper bounds.") + end + + # Check the input lower bounds and upper bounds are reasonable + for i in eachindex(lb) + if lb[i] >= ub[i] + error("Lower bound is not less than upper bound for parameter $(params[i])") + end + end + # construct AbstractMeasureData as data + if use_existing_supports + supports = support(params) + # TODO: think about how to generate reasonable coefficients for given supports + if num_params == 1 + data = DiscreteMeasureData(params[1], ones(size(supports)), supports) + else + data = MultiDiscreteMeasureData(params, ones(size(supports)), supports) + end + else + data = generate_measure_data(params, lb, ub, num_supports, method = eval_method) + end + + # call measure function to construct the measure + return measure(expr, data) +end + +# expectation measure +# TODO: sample from distribution (Distributions package) +function expect(expr::JuMP.AbstractJuMPScalar, + params::Union{ParameterRef, Vector{ParameterRef}, Nothing} = nothing; + num_supports::Int = 50)::MeasureRef + if use_existing_supports + weight(x) = 1 / length(supports(x)); + else + weight(x) = 1 / num_supports; + end + return measure(expr, params, lb, ub; eval_method = eval_method, + num_supports = num_supports, weight_func = weight, + use_existing_supports = use_existing_supports) +end + +# sum measure +function sum(expr::JuMP.AbstractJuMPScalar, + params::Union{ParameterRef, Vector{ParameterRef}, Nothing} = nothing, + lb::Union{Float64, Vector{Float64}} = Float64[], + ub::Union{Float64, Vector{Float64}} = Float64[]; + eval_method::Function = MC_sampling, num_supports::Int = 50, + use_existing_supports::Bool = false)::MeasureRef + return measure(expr, params, lb, ub; eval_method = eval_method, + num_supports = num_supports, + use_existing_supports = use_existing_supports) +end + """ used_by_constraint(mref::MeasureRef)::Bool diff --git a/test/MeasureEvalMethod/methods.jl b/test/MeasureEvalMethod/methods.jl new file mode 100644 index 000000000..9c99764db --- /dev/null +++ b/test/MeasureEvalMethod/methods.jl @@ -0,0 +1,81 @@ +@testset "Default Monte Carlo Sampling" begin + # test univariate Monte Carlo sampling + @testset "MC_sampling (univariate)" begin + (supports, weights) = MC_sampling(0., 2., 5) + @test length(supports) == 5 + @test length(weights) == 5 + @test weights == 0.4 * ones(5) + @test all(supports .>= 0.) + @test all(supports .<= 2.) + end + + # test multivariate Monte Carlo sampling + @testset "MC_sampling (multivariate)" begin + (supports, weights) = MC_sampling([0., 0.], [1., 4.], 5) + @test size(supports)[2] == 5 + @test length(weights) == 5 + @test weights == 0.8 * ones(5) + @test all(supports[1, :] .>= 0.) + @test all(supports[1, :] .<= 1.) + @test all(supports[2, :] .>= 0.) + @test all(supports[2, :] .<= 4.) + end +end + +@testset "Default Quadrature Methods" begin + # test Gauss-Legendre method for bounded interval + @testset "Gauss_Legendre" begin + (supports, weights) = Gauss_Legendre(1., 5., 5) + (expect_supports, expect_weights) = FGQ.gausslegendre(5) + expect_supports = expect_supports * 2 .+ 3 + expect_weights = expect_weights * 2 + @test supports == expect_supports + @test weights == expect_weights + end + + # test Gauss-Hermite method for infinite interval + @testset "Gauss_Hermite" begin + @test_throws ErrorException Gauss_Hermite(0., Inf, 5) + @test_throws ErrorException Gauss_Hermite(-Inf, 0., 5) + @test_throws ErrorException Gauss_Hermite(Inf, Inf, 5) + @test_throws ErrorException Gauss_Hermite(0., 1., 5) + (supports, weights) = Gauss_Hermite(-Inf, Inf, 5) + @test length(supports) == 5 + (expect_supports, expect_weights) = gausshermite(5) + expect_weights = expect_weights .* exp.(expect_supports.^2) + @test length(expect_weights) == 5 + @test supports == expect_supports + @test weights == expect_weights + end + + # test Gauss-Laguerre method for semi-infinite interval + @testset "Gauss_Laguerre" begin + @test_throws ErrorException Gauss_Laguerre(-Inf, Inf, 5) + @test_throws ErrorException Gauss_Laguerre(0., 1., 5) + (supports, weights) = Gauss_Laguerre(1., Inf, 5) + (expect_supports, expect_weights) = gausslaguerre(5) + expect_weights = expect_weights .* exp.(expect_supports) + expect_supports = expect_supports .+ 1. + @test supports == expect_supports + @test weights == expect_weights + (supports, weights) = Gauss_Laguerre(-Inf, 1., 5) + (expect_supports, expect_weights) = gausslaguerre(5) + expect_weights = expect_weights .* exp.(expect_supports) + expect_supports = -expect_supports .+ 1. + @test supports == expect_supports + @test weights == expect_weights + end + + # test the function that transform (semi-)infinite domain to finite domain + @testset "infinite_transform" begin + @test_throws ErrorException infinite_transform(0., 1., 5) + (supports, weights) = infinite_transform(-Inf, Inf, 5, + sub_method = Gauss_Legendre) + (expect_supports, expect_weights) = gausslegendre(5) + expect_supports = MEM._default_x.(expect_supports, -Inf, Inf) + expect_weights = expect_weights .* MEM._default_dx.(expect_weights, + -Inf, Inf) + @test supports == expect_supports + @test weights == expect_weights + end +end diff --git a/test/runtests.jl b/test/runtests.jl index cf11f43c7..0050a3240 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,9 @@ const JuMPC = JuMP.Containers const IOTO = InfiniteOpt.TranscriptionOpt using Test using Random +const MEM = InfiniteOpt.MeasureEvalMethod +using FastGaussQuadrature +const FGQ = FastGaussQuadrature include("utilities.jl") @@ -37,6 +40,10 @@ println("") println("") @time @testset "Measure Methods" begin include("measures.jl") end println("") +@time @testset "Measure Evaluation Methods" begin + include("MeasureEvalMethod/methods.jl") +end +println("") @time @testset "Objective Methods" begin include("objective.jl") end println("") @time @testset "Constraint Methods" begin include("constraints.jl") end From bc6cd678b6a194e8009bffa9fe43ec1a220d3ac0 Mon Sep 17 00:00:00 2001 From: Weiqi Date: Mon, 2 Dec 2019 22:01:06 -0600 Subject: [PATCH 02/10] bug fixes and more tests --- src/InfiniteOpt.jl | 2 +- src/MeasureEvalMethod/methods.jl | 9 ++++--- src/measures.jl | 19 +++++++++---- test/MeasureEvalMethod/methods.jl | 45 ++++++++++++++++++++++++++++--- 4 files changed, 61 insertions(+), 14 deletions(-) diff --git a/src/InfiniteOpt.jl b/src/InfiniteOpt.jl index 2f58fde78..6049e8412 100644 --- a/src/InfiniteOpt.jl +++ b/src/InfiniteOpt.jl @@ -82,7 +82,7 @@ Measure, MeasureRef # Export measure methods export add_measure, measure, measure_function, measure_data, expand, -expand_all_measures! +expand_all_measures!, expect # Export measure evaluation functions export generate_measure_data, MC_sampling, Gauss_Legendre, Gauss_Hermite, diff --git a/src/MeasureEvalMethod/methods.jl b/src/MeasureEvalMethod/methods.jl index 5a9a34d55..5e1b25324 100644 --- a/src/MeasureEvalMethod/methods.jl +++ b/src/MeasureEvalMethod/methods.jl @@ -20,11 +20,11 @@ end # MC sampling from uniform distribution over the interval [lb, ub] function MC_sampling(lb::Float64, ub::Float64, num_supports::Int)::Tuple if lb == -Inf || ub == Inf - (samples, _) = infinite_transform(lb, ub, num_supports) + return infinite_transform(lb, ub, num_supports) else samples = rand(num_supports) .* (ub - lb) .+ lb + return (samples, ones(num_supports) / num_supports * (ub - lb)) end - return (samples, ones(num_supports) / num_supports * (ub - lb)) end # MC sampling - multi-dim version @@ -82,8 +82,8 @@ function infinite_transform(lb::Float64, ub::Float64, num_supports::Int; transform_x::Function = _default_x, transform_dx::Function = _default_dx, t_lb::Float64 = -convert(Float64, lb == -Inf && ub == Inf), - t_ub::Float64 = 1)::Tuple - if lb != -Inf || ub != Inf + t_ub::Float64 = 1.)::Tuple + if lb != -Inf && ub != Inf error("The range is not (semi-)infinite. Use evaluation methods for " * "bounded domains.") end @@ -114,3 +114,4 @@ function _default_dx(t::Float64, lb::Float64, ub::Float64)::Float64 end # TODO: consider truncated distribution +# TODO: consider adding uniform grids diff --git a/src/measures.jl b/src/measures.jl index 02cd1874b..1ea33c7c6 100644 --- a/src/measures.jl +++ b/src/measures.jl @@ -445,7 +445,7 @@ function measure(expr::JuMP.AbstractJuMPScalar, # Default: try to collect all parameters in expr. if isa(params, Nothing) - params = _all_parameter_refs(expr) + params = collect(_all_parameter_refs(expr)) if length(params) == 0 error("No infinite parameters in the expression.") end @@ -457,6 +457,9 @@ function measure(expr::JuMP.AbstractJuMPScalar, error("Parameters cannot be empty.") end num_params = length(params) + if num_params == 1 + params = params[1] + end else num_params = 1 end @@ -475,7 +478,10 @@ function measure(expr::JuMP.AbstractJuMPScalar, error("Some parameter(s) do not have lower bounds. Need to manually " * "input the lower bound values.") end - lb = JuMP.lower_bound.(params) + lb = collect(JuMP.lower_bound.(params)) + if num_params == 1 + lb = lb[1] + end end if length(ub) == 0 params_have_upper_bounds = all(JuMP.has_upper_bound.(params)) @@ -483,7 +489,10 @@ function measure(expr::JuMP.AbstractJuMPScalar, error("Some parameter(s) do not have upper bounds. Need to manually " * "input the upper bound values.") end - ub = JuMP.upper_bound.(params) + ub = collect(JuMP.upper_bound.(params)) + if num_params == 1 + ub = ub[1] + end end # Check the dimension of lb and ub matches number of parameters @@ -519,7 +528,7 @@ end # TODO: sample from distribution (Distributions package) function expect(expr::JuMP.AbstractJuMPScalar, params::Union{ParameterRef, Vector{ParameterRef}, Nothing} = nothing; - num_supports::Int = 50)::MeasureRef + num_supports::Int = 50)::MeasureRef if use_existing_supports weight(x) = 1 / length(supports(x)); else @@ -531,7 +540,7 @@ function expect(expr::JuMP.AbstractJuMPScalar, end # sum measure -function sum(expr::JuMP.AbstractJuMPScalar, +function Base.sum(expr::JuMP.GenericAffExpr{Float64, }, params::Union{ParameterRef, Vector{ParameterRef}, Nothing} = nothing, lb::Union{Float64, Vector{Float64}} = Float64[], ub::Union{Float64, Vector{Float64}} = Float64[]; diff --git a/test/MeasureEvalMethod/methods.jl b/test/MeasureEvalMethod/methods.jl index 9c99764db..f00643ef3 100644 --- a/test/MeasureEvalMethod/methods.jl +++ b/test/MeasureEvalMethod/methods.jl @@ -7,6 +7,9 @@ @test weights == 0.4 * ones(5) @test all(supports .>= 0.) @test all(supports .<= 2.) + (supports, weights) = MC_sampling(-Inf, Inf, 5) + @test length(supports) == 5 + @test length(weights) == 5 end # test multivariate Monte Carlo sampling @@ -22,7 +25,7 @@ end end -@testset "Default Quadrature Methods" begin +@testset "Quadrature Methods" begin # test Gauss-Legendre method for bounded interval @testset "Gauss_Legendre" begin (supports, weights) = Gauss_Legendre(1., 5., 5) @@ -65,17 +68,51 @@ end @test supports == expect_supports @test weights == expect_weights end +end - # test the function that transform (semi-)infinite domain to finite domain +# test the function that transform (semi-)infinite domain to finite domain +# with the default transform function +@testset "Infinite transform for (semi-)infinite domain" begin @testset "infinite_transform" begin @test_throws ErrorException infinite_transform(0., 1., 5) + # Infinite domain (supports, weights) = infinite_transform(-Inf, Inf, 5, sub_method = Gauss_Legendre) (expect_supports, expect_weights) = gausslegendre(5) - expect_supports = MEM._default_x.(expect_supports, -Inf, Inf) - expect_weights = expect_weights .* MEM._default_dx.(expect_weights, + expect_weights = expect_weights .* MEM._default_dx.(expect_supports, -Inf, Inf) + expect_supports = MEM._default_x.(expect_supports, -Inf, Inf) + @test supports == expect_supports + @test weights == expect_weights + # Lower bounded semi-infinite domain + (supports, weights) = infinite_transform(0., Inf, 5, + sub_method = Gauss_Legendre) + (expect_supports, expect_weights) = gausslegendre(5) + expect_supports = (expect_supports .+ 1) ./ 2 + expect_weights = expect_weights ./ 2 + expect_weights = expect_weights .* MEM._default_dx.(expect_supports, + 0., Inf) + expect_supports = MEM._default_x.(expect_supports, 0., Inf) @test supports == expect_supports @test weights == expect_weights + # Upper bounded semi-infinite domain + (supports, weights) = infinite_transform(-Inf, 0., 5, + sub_method = Gauss_Legendre) + (expect_supports, expect_weights) = gausslegendre(5) + expect_supports = (expect_supports .+ 1) ./ 2 + expect_weights = expect_weights ./ 2 + expect_weights = expect_weights .* MEM._default_dx.(expect_supports, + -Inf, 0.) + expect_supports = MEM._default_x.(expect_supports, -Inf, 0.) + @test supports == expect_supports + @test weights == expect_weights + end +end + +@testset "Data Generation" begin + @testset "generate_measure_data (univariate)" begin + end + + @testset "generate_measure_data (multivariate)" begin end end From 8845245f51ef0d8007584e60acd0afc232f96a33 Mon Sep 17 00:00:00 2001 From: Weiqi Date: Thu, 12 Dec 2019 16:23:20 -0600 Subject: [PATCH 03/10] Adopt SparseAxisArray and add tests --- src/InfiniteOpt.jl | 4 +- src/MeasureEvalMethod/MeasureEvalMethod.jl | 5 +- src/MeasureEvalMethod/methods.jl | 88 +++++++--- src/measures.jl | 179 +++++++++++---------- test/MeasureEvalMethod/methods.jl | 29 +++- test/measures.jl | 90 +++++++++++ 6 files changed, 285 insertions(+), 110 deletions(-) diff --git a/src/InfiniteOpt.jl b/src/InfiniteOpt.jl index 6049e8412..b6db47ca3 100644 --- a/src/InfiniteOpt.jl +++ b/src/InfiniteOpt.jl @@ -82,11 +82,11 @@ Measure, MeasureRef # Export measure methods export add_measure, measure, measure_function, measure_data, expand, -expand_all_measures!, expect +expand_all_measures!, expect, support_sum # Export measure evaluation functions export generate_measure_data, MC_sampling, Gauss_Legendre, Gauss_Hermite, -Gauss_Laguerre, infinite_transform +Gauss_Laguerre, infinite_transform, support_formatting # Export transcription datatypes export TranscriptionData, TranscriptionModel diff --git a/src/MeasureEvalMethod/MeasureEvalMethod.jl b/src/MeasureEvalMethod/MeasureEvalMethod.jl index 065efc347..28fe00a42 100644 --- a/src/MeasureEvalMethod/MeasureEvalMethod.jl +++ b/src/MeasureEvalMethod/MeasureEvalMethod.jl @@ -2,12 +2,15 @@ module MeasureEvalMethod # Import the necessary packages. import FastGaussQuadrature +import JuMP +import Distributions +const JuMPC = JuMP.Containers using ..InfiniteOpt # include jl files here include("methods.jl") # Export functions here export generate_measure_data, MC_sampling, Gauss_Legendre, Gauss_Hermite, -Gauss_Laguerre, infinite_transform +Gauss_Laguerre, infinite_transform, support_formatting end # end module diff --git a/src/MeasureEvalMethod/methods.jl b/src/MeasureEvalMethod/methods.jl index 5e1b25324..78cf0337a 100644 --- a/src/MeasureEvalMethod/methods.jl +++ b/src/MeasureEvalMethod/methods.jl @@ -1,24 +1,44 @@ function generate_measure_data(param::InfiniteOpt.ParameterRef, - lb::Float64, ub::Float64, num_supports::Int; + lb::Union{Number, Nothing}, + ub::Union{Number, Nothing}, num_supports::Int; method::Function = MC_sampling, name::String = "", weight_func::Function = InfiniteOpt._w )::InfiniteOpt.DiscreteMeasureData - (supports, coeffs) = method(lb, ub, num_supports) - return InfiniteOpt.DiscreteMeasureData(param, coeffs, supports, name, weight_func) + set = InfiniteOpt._parameter_set(param) + if isa(set, DistributionSet) + (supports, coeffs) = method(set.distribution, num_supports) + else + (supports, coeffs) = method(lb, ub, num_supports) + end + return InfiniteOpt.DiscreteMeasureData(param, coeffs, supports, + name = name, weight_function = weight_func) end -function generate_measure_data(param::Vector{InfiniteOpt.ParameterRef}, - lb::Vector{Float64}, ub::Vector{Float64}, +#function generate_measure_data(params::JuMPC.SparseAxisArray{<:InfiniteOpt.ParameterRef}, +function generate_measure_data(params::AbstractArray{<:InfiniteOpt.ParameterRef}, + lb::Union{JuMPC.SparseAxisArray, Nothing}, + ub::Union{JuMPC.SparseAxisArray, Nothing}, num_supports::Int; method::Function = MC_sampling, name::String = "", weight_func::Function = InfiniteOpt._w )::InfiniteOpt.MultiDiscreteMeasureData - (supports, coeffs) = method(lb, ub, num_supports) - return InfiniteOpt.MultiDiscreteMeasureData(param, coeffs, supports, name, weight_func) + params = convert(JuMPC.SparseAxisArray, params) + set = 1; + for i in params + set = InfiniteOpt._parameter_set(i) + break + end + if isa(set, DistributionSet) + (supports, coeffs) = method(set.distribution, params, num_supports) + else + (supports, coeffs) = method(lb, ub, num_supports) + end + return InfiniteOpt.DiscreteMeasureData(params, coeffs, supports, + name = name, weight_function = weight_func) end # MC sampling from uniform distribution over the interval [lb, ub] -function MC_sampling(lb::Float64, ub::Float64, num_supports::Int)::Tuple +function MC_sampling(lb::Number, ub::Number, num_supports::Int)::Tuple if lb == -Inf || ub == Inf return infinite_transform(lb, ub, num_supports) else @@ -28,17 +48,47 @@ function MC_sampling(lb::Float64, ub::Float64, num_supports::Int)::Tuple end # MC sampling - multi-dim version -function MC_sampling(lb::Vector{Float64}, ub::Vector{Float64}, +function MC_sampling(lb::JuMPC.SparseAxisArray, ub::JuMPC.SparseAxisArray, num_supports::Int)::Tuple - samples = zeros(length(lb), num_supports) + samples_dict = Dict() for i in eachindex(lb) - (samples[i, :], _) = MC_sampling(lb[i], ub[i], num_supports) + (samples_dict[i], _) = MC_sampling(lb[i], ub[i], num_supports) + end + samples = Array{JuMPC.SparseAxisArray, 1}() + for j in 1:num_supports + append!(samples, [JuMP.Containers.SparseAxisArray(Dict(k => samples_dict[k][j] + for k in eachindex(lb)))]) end return (samples, ones(num_supports) / num_supports * prod(ub .- lb)) + end +# MC sampling - Distribution Set +function MC_sampling(dist::Distributions.UnivariateDistribution, + num_supports::Int)::Tuple + samples = rand(dist, num_supports) + return (samples, ones(num_supports) / num_supports) +end + +function MC_sampling(dist::Distributions.MultivariateDistribution, + params::AbstractArray{<:InfiniteOpt.ParameterRef}, + num_supports::Int)::Tuple + samples_matrix = rand(dist, num_supports) + ordered_pairs = sort(collect(params.data), by=x->x.second.index) + samples_dict = Dict() + for i in eachindex(ordered_pairs) + samples_dict[ordered_pairs[i].first] = samples_matrix[i, :] + end + samples = Array{JuMPC.SparseAxisArray, 1}() + for j in 1:num_supports + append!(samples, [JuMP.Containers.SparseAxisArray( + Dict(k.first => samples_dict[k.first][j] + for k in ordered_pairs))]) + end + return (samples, ones(num_supports) / num_supports) +end # default Gaussian quadrature method for bounded domain -function Gauss_Legendre(lb::Float64, ub::Float64, num_supports::Int)::Tuple +function Gauss_Legendre(lb::Number, ub::Number, num_supports::Int)::Tuple (supports, coeffs) = FastGaussQuadrature.gausslegendre(num_supports) supports = (ub - lb) / 2 * supports .+ (ub + lb) / 2 coeffs = (ub - lb) / 2 * coeffs @@ -46,7 +96,7 @@ function Gauss_Legendre(lb::Float64, ub::Float64, num_supports::Int)::Tuple end # default Gaussian quadrature method for infinite domain -function Gauss_Hermite(lb::Float64, ub::Float64, num_supports::Int)::Tuple +function Gauss_Hermite(lb::Number, ub::Number, num_supports::Int)::Tuple if lb != -Inf || ub != Inf error("Lower/upper bound is not infinity. Use other measure evaluation " * "methods.") @@ -57,7 +107,7 @@ function Gauss_Hermite(lb::Float64, ub::Float64, num_supports::Int)::Tuple end # default Gaussian quadrature method for semi-infinite domain -function Gauss_Laguerre(lb::Float64, ub::Float64, num_supports::Int)::Tuple +function Gauss_Laguerre(lb::Number, ub::Number, num_supports::Int)::Tuple if ub == Inf if lb == -Inf error("The range is infinite. Use other measure evaluation methods.") @@ -77,12 +127,12 @@ function Gauss_Laguerre(lb::Float64, ub::Float64, num_supports::Int)::Tuple end # transform (semi-)infinite domain to finite domain -function infinite_transform(lb::Float64, ub::Float64, num_supports::Int; +function infinite_transform(lb::Number, ub::Number, num_supports::Int; sub_method::Function = MC_sampling, transform_x::Function = _default_x, transform_dx::Function = _default_dx, - t_lb::Float64 = -convert(Float64, lb == -Inf && ub == Inf), - t_ub::Float64 = 1.)::Tuple + t_lb::Number = -convert(Number, lb == -Inf && ub == Inf), + t_ub::Number = 1.)::Tuple if lb != -Inf && ub != Inf error("The range is not (semi-)infinite. Use evaluation methods for " * "bounded domains.") @@ -93,7 +143,7 @@ function infinite_transform(lb::Float64, ub::Float64, num_supports::Int; return (supports, coeffs) end -function _default_x(t::Float64, lb::Float64, ub::Float64)::Float64 +function _default_x(t::Number, lb::Number, ub::Number)::Number if lb > -Inf return lb + t / (1 - t) elseif ub < Inf @@ -103,7 +153,7 @@ function _default_x(t::Float64, lb::Float64, ub::Float64)::Float64 end end -function _default_dx(t::Float64, lb::Float64, ub::Float64)::Float64 +function _default_dx(t::Number, lb::Number, ub::Number)::Number if lb > -Inf return 1 / (1 - t)^2 elseif ub < Inf diff --git a/src/measures.jl b/src/measures.jl index 1ea33c7c6..b6f4dc6c1 100644 --- a/src/measures.jl +++ b/src/measures.jl @@ -1,6 +1,3 @@ -const MC = :MC -const GaussLegendre = :GaussLegendre - # Extend Base.copy for new variable types Base.copy(v::MeasureRef, new_model::InfiniteModel) = MeasureRef(new_model, v.index) @@ -414,8 +411,8 @@ end """ measure(expr::JuMP.AbstractJuMPScalar, params::Union{ParameterRef, Vector{ParameterRef}}, - lb::Union{Float64, Vector{Float64}}, - ub::Union{Float64, Vector{Float64}}; + lb::Union{Number, Vector{Number}}, + ub::Union{Number, Vector{Number}}; eval_method::Function, num_supports::Int, weight_func::Function, use_existing_supports::Bool = false)::MeasureRef @@ -436,119 +433,139 @@ default numerical integration schemes. """ # Measure function that takes non-AbstractMeasureData types function measure(expr::JuMP.AbstractJuMPScalar, - params::Union{ParameterRef, Vector{ParameterRef}, Nothing} = nothing, - lb::Union{Float64, Vector{Float64}} = Float64[], - ub::Union{Float64, Vector{Float64}} = Float64[]; + params::Union{ParameterRef, AbstractArray{<:ParameterRef}, + Nothing} = nothing, + lb::Union{Number, AbstractArray{<:Number}, Nothing} = nothing, + ub::Union{Number, AbstractArray{<:Number}, Nothing} = nothing; eval_method::Function = MC_sampling, num_supports::Int = 50, - weight_func::Function = _w, - use_existing_supports::Bool = false)::MeasureRef + weight_func::Function = _w, name = "", + use_existing_supports::Bool = false, + call_from_expect::Bool = false)::MeasureRef - # Default: try to collect all parameters in expr. if isa(params, Nothing) - params = collect(_all_parameter_refs(expr)) + if isa(expr, MeasureRef) + error("Nested call of measure must specify parameters.") + end + params = _all_parameter_refs(expr) if length(params) == 0 error("No infinite parameters in the expression.") + elseif length(params) == 1 + params = params[1] + else + error("Multiple groups of parameters are in the expression. Need to " * + "specify one group of parameters only.") end end - # Check if the parameters are empty - if isa(params, Vector) - if length(params) == 0 - error("Parameters cannot be empty.") - end + if isa(params, ParameterRef) + num_params = 1 + else num_params = length(params) - if num_params == 1 - params = params[1] + if num_params == 0 + error("No infinite parameter is provided.") + elseif num_params == 1 + temp = 1 + for a in params + temp = a + end + params = temp + else + ids = unique(group_id.(params)) + if length(ids) > 1 + error("Multiple groups of parameters are specified.") + end + params = convert(JuMPC.SparseAxisArray, params) + if isa(lb, AbstractArray) + lb = convert(JuMPC.SparseAxisArray, lb) + end + if isa(ub, AbstractArray) + ub = convert(JuMPC.SparseAxisArray, ub) + end end - else - num_params = 1 end - # Check if the parameters belong to multiple groups - ids = unique(group_id.(params)) - if length(ids) > 1 - error("Multiple groups of parameters in the expression. Need to " * - "specify the parameters to integrate over.") + set = 1; + for a in params + set = _parameter_set(a) + break end - # Fill in lower bounds and upper bounds if not given - if length(lb) == 0 - params_have_lower_bounds = all(JuMP.has_lower_bound.(params)) - if !params_have_lower_bounds - error("Some parameter(s) do not have lower bounds. Need to manually " * - "input the lower bound values.") + if num_params > 1 + if isa(lb, Number) + lb = JuMPC.SparseAxisArray(Dict(k => lb for k in keys(params))) end - lb = collect(JuMP.lower_bound.(params)) - if num_params == 1 - lb = lb[1] + if isa(ub, Number) + ub = JuMPC.SparseAxisArray(Dict(k => ub for k in keys(params))) end end - if length(ub) == 0 - params_have_upper_bounds = all(JuMP.has_upper_bound.(params)) - if !params_have_upper_bounds - error("Some parameter(s) do not have upper bounds. Need to manually " * - "input the upper bound values.") + + if isa(set, IntervalSet) + # Fill in lower bounds and upper bounds if not given + if isa(lb, Nothing) || length(lb) == 0 + lb = JuMP.lower_bound.(params) end - ub = collect(JuMP.upper_bound.(params)) - if num_params == 1 - ub = ub[1] + if isa(ub, Nothing) || length(ub) == 0 + ub = JuMP.upper_bound.(params) end - end - # Check the dimension of lb and ub matches number of parameters - if length(lb) != num_params || length(ub) != num_params - error("Number of parameters do not match number of lower bounds or " * - "upper bounds.") - end + # Check the dimension of lb and ub matches number of parameters + if length(lb) != num_params || length(ub) != num_params + error("Number of parameters do not match number of lower bounds or " * + "upper bounds.") + end - # Check the input lower bounds and upper bounds are reasonable - for i in eachindex(lb) - if lb[i] >= ub[i] - error("Lower bound is not less than upper bound for parameter $(params[i])") + # Check the input lower bounds and upper bounds are reasonable + for i in eachindex(lb) + if lb[i] >= ub[i] + error("Lower bound is not less than upper bound for some parameter." * + " Please check the input lower bounds and upper bounds.") + end end end - # construct AbstractMeasureData as data + + # construct data and return measure if we use existing supports if use_existing_supports - supports = support(params) - # TODO: think about how to generate reasonable coefficients for given supports - if num_params == 1 - data = DiscreteMeasureData(params[1], ones(size(supports)), supports) + support = supports(params) + if !isa(lb, Nothing) && !isa(ub, Nothing) + support = [i for i in support if all(i .>= lb) && all(i .<= ub)] + end + + # TODO: generate reasonable coefficients for given supports + len = length(support) + if call_from_expect + data = DiscreteMeasureData(params, ones(len) ./ len, support, + name = name, weight_function = weight_func) else - data = MultiDiscreteMeasureData(params, ones(size(supports)), supports) + data = DiscreteMeasureData(params, ones(len), support, + name = name, weight_function = weight_func) end - else - data = generate_measure_data(params, lb, ub, num_supports, method = eval_method) + return measure(expr, data) end + # construct AbstractMeasureData as data + data = generate_measure_data(params, lb, ub, num_supports, method = eval_method, + name = name, weight_func = weight_func) + # call measure function to construct the measure return measure(expr, data) end # expectation measure -# TODO: sample from distribution (Distributions package) function expect(expr::JuMP.AbstractJuMPScalar, - params::Union{ParameterRef, Vector{ParameterRef}, Nothing} = nothing; - num_supports::Int = 50)::MeasureRef - if use_existing_supports - weight(x) = 1 / length(supports(x)); - else - weight(x) = 1 / num_supports; - end - return measure(expr, params, lb, ub; eval_method = eval_method, - num_supports = num_supports, weight_func = weight, - use_existing_supports = use_existing_supports) + params::Union{ParameterRef, AbstractArray{<:ParameterRef}, + Nothing} = nothing; + num_supports::Int = 50, + use_existing_supports::Bool = false)::MeasureRef + return measure(expr, params, num_supports = num_supports, + name = "expect", use_existing_supports = use_existing_supports, + call_from_expect = true) end # sum measure -function Base.sum(expr::JuMP.GenericAffExpr{Float64, }, - params::Union{ParameterRef, Vector{ParameterRef}, Nothing} = nothing, - lb::Union{Float64, Vector{Float64}} = Float64[], - ub::Union{Float64, Vector{Float64}} = Float64[]; - eval_method::Function = MC_sampling, num_supports::Int = 50, - use_existing_supports::Bool = false)::MeasureRef - return measure(expr, params, lb, ub; eval_method = eval_method, - num_supports = num_supports, - use_existing_supports = use_existing_supports) +function support_sum(expr::JuMP.AbstractJuMPScalar, + params::Union{ParameterRef, AbstractArray{<:ParameterRef}, Nothing} + = nothing)::MeasureRef + return measure(expr, params, use_existing_supports = true, name = "sum") end """ diff --git a/test/MeasureEvalMethod/methods.jl b/test/MeasureEvalMethod/methods.jl index f00643ef3..40f6a99c1 100644 --- a/test/MeasureEvalMethod/methods.jl +++ b/test/MeasureEvalMethod/methods.jl @@ -14,14 +14,14 @@ # test multivariate Monte Carlo sampling @testset "MC_sampling (multivariate)" begin - (supports, weights) = MC_sampling([0., 0.], [1., 4.], 5) - @test size(supports)[2] == 5 + lb = convert(JuMPC.SparseAxisArray, [0., 0.]) + ub = convert(JuMPC.SparseAxisArray, [1., 4.]) + (supports, weights) = MC_sampling(lb, ub, 5) + @test length(supports) == 5 @test length(weights) == 5 @test weights == 0.8 * ones(5) - @test all(supports[1, :] .>= 0.) - @test all(supports[1, :] .<= 1.) - @test all(supports[2, :] .>= 0.) - @test all(supports[2, :] .<= 4.) + @test all([i >= lb for i in supports]) + @test all([i <= ub for i in supports]) end end @@ -37,6 +37,7 @@ end end # test Gauss-Hermite method for infinite interval + # at the compiling step FGQ.gausshermite creates huge allocation @testset "Gauss_Hermite" begin @test_throws ErrorException Gauss_Hermite(0., Inf, 5) @test_throws ErrorException Gauss_Hermite(-Inf, 0., 5) @@ -44,7 +45,7 @@ end @test_throws ErrorException Gauss_Hermite(0., 1., 5) (supports, weights) = Gauss_Hermite(-Inf, Inf, 5) @test length(supports) == 5 - (expect_supports, expect_weights) = gausshermite(5) + (expect_supports, expect_weights) = FGQ.gausshermite(5) expect_weights = expect_weights .* exp.(expect_supports.^2) @test length(expect_weights) == 5 @test supports == expect_supports @@ -110,9 +111,23 @@ end end @testset "Data Generation" begin + m = InfiniteModel() + @infinite_parameter(m, x in [0., 1.]) + @infinite_parameter(m, y[1:2] in [0., 1.]) @testset "generate_measure_data (univariate)" begin + @test_throws ErrorException generate_measure_data(x, 1., 2., 10) + measure_data = generate_measure_data(x, 0.3, 0.7, 10, method = Gauss_Legendre) + (expect_supports, expect_weights) = Gauss_Legendre(0.3, 0.7, 10) + @test measure_data.supports == expect_supports + @test measure_data.coefficients == expect_weights end @testset "generate_measure_data (multivariate)" begin + lb = convert(JuMPC.SparseAxisArray, [0.3, 0.2]) + ub = convert(JuMPC.SparseAxisArray, [0.5, 0.9]) + measure_data = generate_measure_data(y, lb, ub, 10) + @test length(measure_data.supports) == 10 + @test length(measure_data.coefficients) == 10 + @test sum(measure_data.coefficients) == 0.14 end end diff --git a/test/measures.jl b/test/measures.jl index 234f4ffd0..7e52b2657 100644 --- a/test/measures.jl +++ b/test/measures.jl @@ -557,3 +557,93 @@ end m.meas_in_objective[JuMP.index(mref)] = false end end +@testset "User Definition w/o Measure Data" begin + m = InfiniteModel() + dist1 = Normal(0., 1.) + dist2 = MvNormal([0., 0.], [1. 0.;0. 10.]) + @infinite_parameter(m, 0 <= par <= 1) + @infinite_parameter(m, 0 <= par2 <= 1) + @infinite_parameter(m, 0 <= par3 <= 1) + @infinite_parameter(m, 0 <= pars1[1:2] <= 1, container = SparseAxisArray) + @infinite_parameter(m, 0 <= pars2[("a", "b")] <= 1) + @infinite_parameter(m, 0 <= pars3[1:2] <= 1) + @infinite_parameter(m, rp in dist1) + @infinite_parameter(m, rp2[1:2] in dist2) + @infinite_variable(m, inf(par)) + @infinite_variable(m, inf2(par, par2)) + @infinite_variable(m, inf3(par3)) + @infinite_variable(m, inf4(pars1)) + @infinite_variable(m, inf5(pars1, pars2)) + @infinite_variable(m, inf6(pars2)) + @infinite_variable(m, inf7(pars3)) + @infinite_variable(m, inf8(rp)) + @infinite_variable(m, inf9(rp2)) + @infinite_variable(m, inf10(rp, rp2)) + @hold_variable(m, x) + + # test measure that does not use AbstractMeasureData inputs + @testset "measure (no AbstractMeasureData)" begin + meas1 = measure(inf, num_supports = 5, eval_method = Gauss_Legendre) + (expected_supps, expected_coeffs) = FGQ.gausslegendre(5) + expected_supps = expected_supps .* 0.5 .+ 0.5 + expected_coeffs = expected_coeffs .* 0.5 + @test all(measure_data(meas1).supports .== expected_supps) + @test all(measure_data(meas1).coefficients .== expected_coeffs) + + add_supports(par2, [0.3, 0.7]) + meas2 = measure(inf2, par2, use_existing_supports = true) + @test measure_data(meas2).parameter_ref == par2 + @test measure_data(meas2).supports == [0.3, 0.7] + meas2 = measure(inf2, par2, 0.5, 0.9, use_existing_supports = true) + @test measure_data(meas2).supports == [0.7] + + meas3 = measure(inf4, num_supports = 5) + @test pars1[1] in measure_data(meas3).parameter_ref + @test pars1[2] in measure_data(meas3).parameter_ref + + meas4 = measure(inf5, pars2, num_supports = 5) + @test pars2["a"] in measure_data(meas4).parameter_ref + @test pars2["b"] in measure_data(meas4).parameter_ref + + meas5 = measure(inf6, use_existing_supports = true) + @test measure_data(meas4).supports == measure_data(meas5).supports + + # test errors + @test_throws ErrorException measure(x) + @test_throws ErrorException measure(inf, ParameterRef[]) + @test_throws ErrorException measure(inf2) + @test_throws ErrorException measure(inf2, par, 1., 3.) + @test_throws ErrorException measure(inf2, par, [0., 1.]) + @test_throws ErrorException measure(inf2, par, 0., [1., 1.]) + @test_throws ErrorException measure(inf2, par, 0.5, 0.) + @test_throws ErrorException measure(inf7, use_existing_supports = true) + @test_throws ErrorException measure(meas1) + end + # test support_sum + @testset "support_sum" begin + sum1 = support_sum(inf2, par2) + @test measure_data(sum1).parameter_ref == par2 + @test measure_data(sum1).supports == [0.3, 0.7] + @test name(sum1) == "sum(inf2(par, par2))" + add_supports(pars3[1], [0.3, 0.7]) + add_supports(pars3[2], [0.3, 0.7]) + sum2 = support_sum(inf7) + @test pars3[1] in measure_data(sum2).parameter_ref + @test pars3[2] in measure_data(sum2).parameter_ref + supps = [JuMPC.SparseAxisArray(Dict([((1,),0.3), ((2,), 0.3)])), + JuMPC.SparseAxisArray(Dict([((1,),0.7), ((2,), 0.7)]))] + @test measure_data(sum2).supports == supps + end + # test expectation measure + @testset "expect" begin + expect1 = expect(inf8, num_supports = 5) + expect2 = expect(inf9, num_supports = 5) + @test_throws ErrorException expect(inf10) + expect3 = expect(inf10, rp, use_existing_supports = true) + check1 = expect(inf8, use_existing_supports = true) + check2 = expect(inf9, use_existing_supports = true) + @test measure_data(expect1).supports == measure_data(check1).supports + @test measure_data(expect2).supports == measure_data(check2).supports + @test measure_data(expect3).supports == measure_data(expect1).supports + end +end From e14b978bb7f57ef1b1a2618561727c6b37b45f31 Mon Sep 17 00:00:00 2001 From: Weiqi Date: Fri, 13 Dec 2019 10:58:52 -0600 Subject: [PATCH 04/10] add FastGaussQuadrature dependency --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 83af9e19a..d447d2c2e 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" [compat] julia = "1" From 7e5ba91004333e985b6f04900f649fe0faff4788 Mon Sep 17 00:00:00 2001 From: Weiqi Date: Sun, 15 Dec 2019 08:43:28 -0600 Subject: [PATCH 05/10] add more tests --- test/measures.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/measures.jl b/test/measures.jl index 7e52b2657..ee88c74c8 100644 --- a/test/measures.jl +++ b/test/measures.jl @@ -596,6 +596,8 @@ end @test measure_data(meas2).supports == [0.3, 0.7] meas2 = measure(inf2, par2, 0.5, 0.9, use_existing_supports = true) @test measure_data(meas2).supports == [0.7] + meas2 = measure(inf2, [par2], [0.5], [0.9], use_existing_supports = true) + @test measure_data(meas2).supports == [0.7] meas3 = measure(inf4, num_supports = 5) @test pars1[1] in measure_data(meas3).parameter_ref @@ -612,6 +614,7 @@ end @test_throws ErrorException measure(x) @test_throws ErrorException measure(inf, ParameterRef[]) @test_throws ErrorException measure(inf2) + @test_throws ErrorException measure(inf2, [par, par2]) @test_throws ErrorException measure(inf2, par, 1., 3.) @test_throws ErrorException measure(inf2, par, [0., 1.]) @test_throws ErrorException measure(inf2, par, 0., [1., 1.]) From 626b0cfc939238f66c87106ab4c47437ced37af9 Mon Sep 17 00:00:00 2001 From: Weiqi Date: Sun, 15 Dec 2019 10:37:02 -0600 Subject: [PATCH 06/10] fix lb/ub types and tests --- src/measures.jl | 9 +++++++++ test/measures.jl | 13 ++++++++++--- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/src/measures.jl b/src/measures.jl index b6f4dc6c1..2b25f98ac 100644 --- a/src/measures.jl +++ b/src/measures.jl @@ -514,6 +514,15 @@ function measure(expr::JuMP.AbstractJuMPScalar, "upper bounds.") end + if num_params == 1 + if isa(lb, AbstractArray) + lb = lb[findfirst(x->isa(x, Number), lb)] + end + if isa(ub, AbstractArray) + ub = ub[findfirst(x->isa(x, Number), ub)] + end + end + # Check the input lower bounds and upper bounds are reasonable for i in eachindex(lb) if lb[i] >= ub[i] diff --git a/test/measures.jl b/test/measures.jl index ee88c74c8..a9e2ea837 100644 --- a/test/measures.jl +++ b/test/measures.jl @@ -610,6 +610,15 @@ end meas5 = measure(inf6, use_existing_supports = true) @test measure_data(meas4).supports == measure_data(meas5).supports + add_supports(pars3[1], [0.3, 0.7]) + add_supports(pars3[2], [0.3, 0.7]) + meas6 = measure(inf7, pars3, 0.5, 1.0, use_existing_supports = true) + @test measure_data(meas6).supports == [JuMPC.SparseAxisArray( + Dict([((1,),0.7), ((2,), 0.7)]))] + meas6 = measure(inf7, pars3, [0.5, 0.5], [1.0, 1.0], use_existing_supports = true) + @test measure_data(meas6).supports == [JuMPC.SparseAxisArray( + Dict([((1,),0.7), ((2,), 0.7)]))] + # test errors @test_throws ErrorException measure(x) @test_throws ErrorException measure(inf, ParameterRef[]) @@ -619,7 +628,7 @@ end @test_throws ErrorException measure(inf2, par, [0., 1.]) @test_throws ErrorException measure(inf2, par, 0., [1., 1.]) @test_throws ErrorException measure(inf2, par, 0.5, 0.) - @test_throws ErrorException measure(inf7, use_existing_supports = true) + @test_throws ErrorException measure(inf8, use_existing_supports = true) @test_throws ErrorException measure(meas1) end # test support_sum @@ -628,8 +637,6 @@ end @test measure_data(sum1).parameter_ref == par2 @test measure_data(sum1).supports == [0.3, 0.7] @test name(sum1) == "sum(inf2(par, par2))" - add_supports(pars3[1], [0.3, 0.7]) - add_supports(pars3[2], [0.3, 0.7]) sum2 = support_sum(inf7) @test pars3[1] in measure_data(sum2).parameter_ref @test pars3[2] in measure_data(sum2).parameter_ref From 3af71a7540ee86eff33f224503eda72b9172bee7 Mon Sep 17 00:00:00 2001 From: Weiqi Date: Mon, 16 Dec 2019 14:00:49 -0600 Subject: [PATCH 07/10] extendible measure data generation and truncated distribution --- src/InfiniteOpt.jl | 6 +- .../MeasureEvalMethods.jl} | 4 +- .../methods.jl | 109 +++++++++++------- src/measures.jl | 15 +-- .../methods.jl | 26 ++++- test/runtests.jl | 13 ++- 6 files changed, 106 insertions(+), 67 deletions(-) rename src/{MeasureEvalMethod/MeasureEvalMethod.jl => MeasureEvalMethods/MeasureEvalMethods.jl} (75%) rename src/{MeasureEvalMethod => MeasureEvalMethods}/methods.jl (56%) rename test/{MeasureEvalMethod => MeasureEvalMethods}/methods.jl (82%) diff --git a/src/InfiniteOpt.jl b/src/InfiniteOpt.jl index b6db47ca3..ed96ede75 100644 --- a/src/InfiniteOpt.jl +++ b/src/InfiniteOpt.jl @@ -29,8 +29,8 @@ include("utilities.jl") include("TranscriptionOpt/TranscriptionOpt.jl") using .TranscriptionOpt -include("MeasureEvalMethod/MeasureEvalMethod.jl") -using .MeasureEvalMethod +include("MeasureEvalMethods/MeasureEvalMethods.jl") +using .MeasureEvalMethods # Export model object datatype export InfiniteModel @@ -86,7 +86,7 @@ expand_all_measures!, expect, support_sum # Export measure evaluation functions export generate_measure_data, MC_sampling, Gauss_Legendre, Gauss_Hermite, -Gauss_Laguerre, infinite_transform, support_formatting +Gauss_Laguerre, infinite_transform, support_formatting, measure_dispatch # Export transcription datatypes export TranscriptionData, TranscriptionModel diff --git a/src/MeasureEvalMethod/MeasureEvalMethod.jl b/src/MeasureEvalMethods/MeasureEvalMethods.jl similarity index 75% rename from src/MeasureEvalMethod/MeasureEvalMethod.jl rename to src/MeasureEvalMethods/MeasureEvalMethods.jl index 28fe00a42..24c636308 100644 --- a/src/MeasureEvalMethod/MeasureEvalMethod.jl +++ b/src/MeasureEvalMethods/MeasureEvalMethods.jl @@ -1,4 +1,4 @@ -module MeasureEvalMethod +module MeasureEvalMethods # Import the necessary packages. import FastGaussQuadrature @@ -11,6 +11,6 @@ using ..InfiniteOpt include("methods.jl") # Export functions here export generate_measure_data, MC_sampling, Gauss_Legendre, Gauss_Hermite, -Gauss_Laguerre, infinite_transform, support_formatting +Gauss_Laguerre, infinite_transform, support_formatting, measure_dispatch end # end module diff --git a/src/MeasureEvalMethod/methods.jl b/src/MeasureEvalMethods/methods.jl similarity index 56% rename from src/MeasureEvalMethod/methods.jl rename to src/MeasureEvalMethods/methods.jl index 78cf0337a..1604b485f 100644 --- a/src/MeasureEvalMethod/methods.jl +++ b/src/MeasureEvalMethods/methods.jl @@ -1,44 +1,69 @@ -function generate_measure_data(param::InfiniteOpt.ParameterRef, - lb::Union{Number, Nothing}, - ub::Union{Number, Nothing}, num_supports::Int; +function generate_measure_data(params::Union{InfiniteOpt.ParameterRef, + AbstractArray{<:InfiniteOpt.ParameterRef}}, + num_supports::Int, + lb::Union{Number, JuMPC.SparseAxisArray, Nothing} = nothing, + ub::Union{Number, JuMPC.SparseAxisArray, Nothing} = nothing; method::Function = MC_sampling, name::String = "", - weight_func::Function = InfiniteOpt._w - )::InfiniteOpt.DiscreteMeasureData - set = InfiniteOpt._parameter_set(param) - if isa(set, DistributionSet) - (supports, coeffs) = method(set.distribution, num_supports) + weight_func::Function = InfiniteOpt._w, kwargs... + )::InfiniteOpt.AbstractMeasureData + if isa(params, InfiniteOpt.ParameterRef) + set = InfiniteOpt._parameter_set(params) else - (supports, coeffs) = method(lb, ub, num_supports) + params = convert(JuMPC.SparseAxisArray, params) + set = InfiniteOpt._parameter_set(first(params)) end - return InfiniteOpt.DiscreteMeasureData(param, coeffs, supports, + (supports, coeffs) = measure_dispatch(set, params, num_supports, lb, ub, method; kwargs...) + return InfiniteOpt.DiscreteMeasureData(params, coeffs, supports, name = name, weight_function = weight_func) end -#function generate_measure_data(params::JuMPC.SparseAxisArray{<:InfiniteOpt.ParameterRef}, -function generate_measure_data(params::AbstractArray{<:InfiniteOpt.ParameterRef}, - lb::Union{JuMPC.SparseAxisArray, Nothing}, - ub::Union{JuMPC.SparseAxisArray, Nothing}, - num_supports::Int; - method::Function = MC_sampling, name::String = "", - weight_func::Function = InfiniteOpt._w - )::InfiniteOpt.MultiDiscreteMeasureData - params = convert(JuMPC.SparseAxisArray, params) - set = 1; - for i in params - set = InfiniteOpt._parameter_set(i) - break - end - if isa(set, DistributionSet) - (supports, coeffs) = method(set.distribution, params, num_supports) - else - (supports, coeffs) = method(lb, ub, num_supports) +function measure_dispatch(set::InfiniteOpt.IntervalSet, + params::Union{InfiniteOpt.ParameterRef, + AbstractArray{<:InfiniteOpt.ParameterRef}}, + num_supports::Int, + lb::Union{Number, JuMPC.SparseAxisArray, Nothing}, + ub::Union{Number, JuMPC.SparseAxisArray, Nothing}, + method::Function; kwargs...)::Tuple + return method(lb, ub, num_supports; kwargs...) +end + +function measure_dispatch(set::InfiniteOpt.DistributionSet, + params::Union{InfiniteOpt.ParameterRef, + AbstractArray{<:InfiniteOpt.ParameterRef}}, + num_supports::Int, + lb::Union{Number, JuMPC.SparseAxisArray, Nothing}, + ub::Union{Number, JuMPC.SparseAxisArray, Nothing}, + method::Function; kwargs...)::Tuple + dist = set.distribution + # create truncated distribution if necessary + if !isa(lb, Nothing) || !isa(ub, Nothing) + if isa(dist, Distributions.MultivariateDistribution) + @warn("Truncated distribution for multivariate distribution is " * + "not supported. Lower bounds and upper bounds are ignored.") + else + isa(lb, Number) ? temp_lb = lb : temp_lb = -Inf + isa(ub, Number) ? temp_ub = ub : temp_ub = Inf + dist = Distributions.Truncated(dist, temp_lb, temp_ub) + end end - return InfiniteOpt.DiscreteMeasureData(params, coeffs, supports, - name = name, weight_function = weight_func) + return method(dist, params, num_supports; kwargs...) +end + +# fallback +function measure_dispatch(set::InfiniteOpt.AbstractInfiniteSet, + params::Union{InfiniteOpt.ParameterRef, + AbstractArray{<:InfiniteOpt.ParameterRef}}, + num_supports::Int, + lb::Union{Number, JuMPC.SparseAxisArray, Nothing}, + ub::Union{Number, JuMPC.SparseAxisArray, Nothing}, + method::Function; kwargs...)::Tuple + error("Measure dispatch function is not extended for parameters in sets " * + "of type $(typeof(set)).") + return end # MC sampling from uniform distribution over the interval [lb, ub] -function MC_sampling(lb::Number, ub::Number, num_supports::Int)::Tuple +function MC_sampling(lb::Number, ub::Number, num_supports::Int; kwargs...)::Tuple if lb == -Inf || ub == Inf return infinite_transform(lb, ub, num_supports) else @@ -49,30 +74,30 @@ end # MC sampling - multi-dim version function MC_sampling(lb::JuMPC.SparseAxisArray, ub::JuMPC.SparseAxisArray, - num_supports::Int)::Tuple + num_supports::Int; kwargs...)::Tuple samples_dict = Dict() for i in eachindex(lb) (samples_dict[i], _) = MC_sampling(lb[i], ub[i], num_supports) end - samples = Array{JuMPC.SparseAxisArray, 1}() + samples = Array{JuMPC.SparseAxisArray, 1}(undef, num_supports) for j in 1:num_supports - append!(samples, [JuMP.Containers.SparseAxisArray(Dict(k => samples_dict[k][j] - for k in eachindex(lb)))]) + samples[j] = JuMP.Containers.SparseAxisArray(Dict(k => samples_dict[k][j] + for k in eachindex(lb))) end return (samples, ones(num_supports) / num_supports * prod(ub .- lb)) - end # MC sampling - Distribution Set function MC_sampling(dist::Distributions.UnivariateDistribution, - num_supports::Int)::Tuple + param::InfiniteOpt.ParameterRef, + num_supports::Int; kwargs...)::Tuple samples = rand(dist, num_supports) return (samples, ones(num_supports) / num_supports) end function MC_sampling(dist::Distributions.MultivariateDistribution, params::AbstractArray{<:InfiniteOpt.ParameterRef}, - num_supports::Int)::Tuple + num_supports::Int; kwargs...)::Tuple samples_matrix = rand(dist, num_supports) ordered_pairs = sort(collect(params.data), by=x->x.second.index) samples_dict = Dict() @@ -88,7 +113,7 @@ function MC_sampling(dist::Distributions.MultivariateDistribution, return (samples, ones(num_supports) / num_supports) end # default Gaussian quadrature method for bounded domain -function Gauss_Legendre(lb::Number, ub::Number, num_supports::Int)::Tuple +function Gauss_Legendre(lb::Number, ub::Number, num_supports::Int; kwargs...)::Tuple (supports, coeffs) = FastGaussQuadrature.gausslegendre(num_supports) supports = (ub - lb) / 2 * supports .+ (ub + lb) / 2 coeffs = (ub - lb) / 2 * coeffs @@ -96,7 +121,7 @@ function Gauss_Legendre(lb::Number, ub::Number, num_supports::Int)::Tuple end # default Gaussian quadrature method for infinite domain -function Gauss_Hermite(lb::Number, ub::Number, num_supports::Int)::Tuple +function Gauss_Hermite(lb::Number, ub::Number, num_supports::Int; kwargs...)::Tuple if lb != -Inf || ub != Inf error("Lower/upper bound is not infinity. Use other measure evaluation " * "methods.") @@ -107,7 +132,7 @@ function Gauss_Hermite(lb::Number, ub::Number, num_supports::Int)::Tuple end # default Gaussian quadrature method for semi-infinite domain -function Gauss_Laguerre(lb::Number, ub::Number, num_supports::Int)::Tuple +function Gauss_Laguerre(lb::Number, ub::Number, num_supports::Int; kwargs...)::Tuple if ub == Inf if lb == -Inf error("The range is infinite. Use other measure evaluation methods.") @@ -132,7 +157,7 @@ function infinite_transform(lb::Number, ub::Number, num_supports::Int; transform_x::Function = _default_x, transform_dx::Function = _default_dx, t_lb::Number = -convert(Number, lb == -Inf && ub == Inf), - t_ub::Number = 1.)::Tuple + t_ub::Number = 1., kwargs...)::Tuple if lb != -Inf && ub != Inf error("The range is not (semi-)infinite. Use evaluation methods for " * "bounded domains.") diff --git a/src/measures.jl b/src/measures.jl index 2b25f98ac..cb17e3b5f 100644 --- a/src/measures.jl +++ b/src/measures.jl @@ -464,11 +464,7 @@ function measure(expr::JuMP.AbstractJuMPScalar, if num_params == 0 error("No infinite parameter is provided.") elseif num_params == 1 - temp = 1 - for a in params - temp = a - end - params = temp + params = first(params) else ids = unique(group_id.(params)) if length(ids) > 1 @@ -484,12 +480,7 @@ function measure(expr::JuMP.AbstractJuMPScalar, end end - set = 1; - for a in params - set = _parameter_set(a) - break - end - + set = _parameter_set(first(params)) if num_params > 1 if isa(lb, Number) lb = JuMPC.SparseAxisArray(Dict(k => lb for k in keys(params))) @@ -552,7 +543,7 @@ function measure(expr::JuMP.AbstractJuMPScalar, end # construct AbstractMeasureData as data - data = generate_measure_data(params, lb, ub, num_supports, method = eval_method, + data = generate_measure_data(params, num_supports, lb, ub, method = eval_method, name = name, weight_func = weight_func) # call measure function to construct the measure diff --git a/test/MeasureEvalMethod/methods.jl b/test/MeasureEvalMethods/methods.jl similarity index 82% rename from test/MeasureEvalMethod/methods.jl rename to test/MeasureEvalMethods/methods.jl index 40f6a99c1..6c6228164 100644 --- a/test/MeasureEvalMethod/methods.jl +++ b/test/MeasureEvalMethods/methods.jl @@ -114,20 +114,40 @@ end m = InfiniteModel() @infinite_parameter(m, x in [0., 1.]) @infinite_parameter(m, y[1:2] in [0., 1.]) + dist1 = Normal(0., 1.) + dist2 = MvNormal([0., 0.], [1. 0.;0. 10.]) + @infinite_parameter(m, β in dist1) + @infinite_parameter(m, ω[1:2] in dist2) @testset "generate_measure_data (univariate)" begin - @test_throws ErrorException generate_measure_data(x, 1., 2., 10) - measure_data = generate_measure_data(x, 0.3, 0.7, 10, method = Gauss_Legendre) + @test_throws ErrorException generate_measure_data(x, 10, 1., 2.) + measure_data = generate_measure_data(x, 10, 0.3, 0.7, method = Gauss_Legendre) (expect_supports, expect_weights) = Gauss_Legendre(0.3, 0.7, 10) @test measure_data.supports == expect_supports @test measure_data.coefficients == expect_weights + measure_data = generate_measure_data(β, 10, -1., Inf) + @test length(measure_data.supports) == 10 + @test all([i >= -1. for i in measure_data.supports]) + measure_data = generate_measure_data(β, 10, -Inf, 1.) + @test length(measure_data.supports) == 10 + @test all([i <= 1. for i in measure_data.supports]) end @testset "generate_measure_data (multivariate)" begin lb = convert(JuMPC.SparseAxisArray, [0.3, 0.2]) ub = convert(JuMPC.SparseAxisArray, [0.5, 0.9]) - measure_data = generate_measure_data(y, lb, ub, 10) + measure_data = generate_measure_data(y, 10, lb, ub) @test length(measure_data.supports) == 10 @test length(measure_data.coefficients) == 10 @test sum(measure_data.coefficients) == 0.14 + warn = "Truncated distribution for multivariate distribution is " * + "not supported. Lower bounds and upper bounds are ignored." + @test_logs (:warn, warn) generate_measure_data(ω, 10, + convert(JuMPC.SparseAxisArray, [-1., -10.]), + convert(JuMPC.SparseAxisArray, [1., 10.])) + end + + @testset "measure_dispatch for unextended types" begin + @test_throws ErrorException measure_dispatch(BadSet(), x, 10, + 0., 1., MC_sampling) end end diff --git a/test/runtests.jl b/test/runtests.jl index 0050a3240..d510a2527 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,7 +7,7 @@ const JuMPC = JuMP.Containers const IOTO = InfiniteOpt.TranscriptionOpt using Test using Random -const MEM = InfiniteOpt.MeasureEvalMethod +const MEM = InfiniteOpt.MeasureEvalMethods using FastGaussQuadrature const FGQ = FastGaussQuadrature @@ -17,6 +17,7 @@ include("utilities.jl") println("-----------------------------------------") println("----------------UNIT TESTS---------------") println("-----------------------------------------") +#= @time @testset "Datatypes" begin include("datatypes.jl") end println("") @time @testset "Utilities" begin include("utility_tests.jl") end @@ -38,12 +39,14 @@ println("") println("") @time @testset "Expression Methods" begin include("expressions.jl") end println("") -@time @testset "Measure Methods" begin include("measures.jl") end -println("") +=# @time @testset "Measure Evaluation Methods" begin - include("MeasureEvalMethod/methods.jl") + include("MeasureEvalMethods/methods.jl") end println("") +@time @testset "Measure Methods" begin include("measures.jl") end +println("") +#= @time @testset "Objective Methods" begin include("objective.jl") end println("") @time @testset "Constraint Methods" begin include("constraints.jl") end @@ -68,7 +71,7 @@ println("") # TODO add extension tests # TODO add involved deletion tests - +=# println("-----------------------------------------") println("-------------TESTING COMPLETE------------") println("-----------------------------------------") From cf20f56f78eeb015ddd13bd5a69dd1453f5f4f59 Mon Sep 17 00:00:00 2001 From: Weiqi Date: Mon, 16 Dec 2019 16:36:45 -0600 Subject: [PATCH 08/10] add docstrings --- src/MeasureEvalMethods/methods.jl | 156 ++++++++++++++++++++++++++++++ src/measures.jl | 88 +++++++++++++++-- 2 files changed, 235 insertions(+), 9 deletions(-) diff --git a/src/MeasureEvalMethods/methods.jl b/src/MeasureEvalMethods/methods.jl index 1604b485f..e6aa0bbea 100644 --- a/src/MeasureEvalMethods/methods.jl +++ b/src/MeasureEvalMethods/methods.jl @@ -1,3 +1,28 @@ +""" + generate_measure_data(params::Union{InfiniteOpt.ParameterRef, + AbstractArray{<:InfiniteOpt.ParameterRef}}, + num_supports::Int, + lb::Union{Number, JuMPC.SparseAxisArray, Nothing} = nothing, + ub::Union{Number, JuMPC.SparseAxisArray, Nothing} = nothing; + method::Function = MC_sampling, name::String = "", + weight_func::Function = InfiniteOpt._w, kwargs... + )::InfiniteOpt.AbstractMeasureData +Generate an [`AbstractMeasureData`](@ref) object that automatically generate +supports based on a set of given information. The information is required to +include parameters and number of supports. Other optional information to input +includes lower and upper bounds, measure name, support generation method, and +weight functions. The users could supply extra keyword arguments if necessary +for their custom support generation methods. + +**Example** +```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel()) +julia> @infinite_parameter(m, x in [0., 1.]) +x + +julia> measure_data = generate_measure_data(x, 3, 0.3, 0.7, method = Gauss_Legendre) +DiscreteMeasureData(x, [0.1111111111111111, 0.17777777777777776, 0.1111111111111111], [0.3450806661517033, 0.5, 0.6549193338482967], "", InfiniteOpt._w) +``` +""" function generate_measure_data(params::Union{InfiniteOpt.ParameterRef, AbstractArray{<:InfiniteOpt.ParameterRef}}, num_supports::Int, @@ -62,6 +87,28 @@ function measure_dispatch(set::InfiniteOpt.AbstractInfiniteSet, return end +""" + MC_sampling(lb::Union{JuMPC.SparseAxisArray, Number}, + ub::Union{JuMPC.SparseAxisArray, Number}, + num_supports::Int; kwargs...)::Tuple + +Return a tuple that contains supports and coefficients generated by Monte Carlo +sampling from a uniform distribution between the lower and upper bounds provided. + +**Example** +```jldoctest; setup = :(using InfiniteOpt, Random; Random.seed!(0)) +julia> (supps, coeffs) = MC_sampling(0., 1., 5) +([0.8236475079774124, 0.9103565379264364, 0.16456579813368521, 0.17732884646626457, 0.278880109331201], [0.2, 0.2, 0.2, 0.2, 0.2]) + +julia> supps +5-element Array{Float64,1}: + 0.8236475079774124 + 0.9103565379264364 + 0.16456579813368521 + 0.17732884646626457 + 0.278880109331201 +``` +""" # MC sampling from uniform distribution over the interval [lb, ub] function MC_sampling(lb::Number, ub::Number, num_supports::Int; kwargs...)::Tuple if lb == -Inf || ub == Inf @@ -87,6 +134,26 @@ function MC_sampling(lb::JuMPC.SparseAxisArray, ub::JuMPC.SparseAxisArray, return (samples, ones(num_supports) / num_supports * prod(ub .- lb)) end +""" + MC_sampling(dist::Distributions.NonMatrixDistribution, + param::InfiniteOpt.ParameterRef, + num_supports::Int; kwargs...)::Tuple + +Return a tuple that contains supports and coefficients generated by Monte Carlo +sampling from a given distribution. + +**Example** +```jldoctest; setup = :(using InfiniteOpt, Distributions, Random; m = InfiniteModel(seed = true)) +julia> dist = Normal(0., 1.) +Normal{Float64}(μ=0.0, σ=1.0) + +julia> @infinite_parameter(m, x in dist) +x + +julia> MC_sampling(dist, x, 10) +([0.6791074260357777, 0.8284134829000359, -0.3530074003005963, -0.13485387193052173, 0.5866170746331097, 0.29733585084941616, 0.06494754854834232, -0.10901738508171745, -0.514210390833322, 1.5743302021369892], [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]) +``` +""" # MC sampling - Distribution Set function MC_sampling(dist::Distributions.UnivariateDistribution, param::InfiniteOpt.ParameterRef, @@ -112,6 +179,28 @@ function MC_sampling(dist::Distributions.MultivariateDistribution, end return (samples, ones(num_supports) / num_supports) end + +""" + Gauss_Legendre(lb::Number, ub::Number, num_supports::Int; kwargs...)::Tuple + +Return a tuple that contains supports and coefficients generated using +Gauss-Legendre quadrature method. This is useful for univariate parameter in a +finite interval. + +**Example** +```jldoctest; setup = :(using InfiniteOpt) +julia> (supps, coeffs) = Gauss_Legendre(0., 1., 5) +([0.04691007703066802, 0.23076534494715845, 0.5, 0.7692346550528415, 0.9530899229693319], [0.11846344252809454, 0.23931433524968324, 0.28444444444444444, 0.23931433524968324, 0.11846344252809454]) + +julia> supps +5-element Array{Float64,1}: + 0.04691007703066802 + 0.23076534494715845 + 0.5 + 0.7692346550528415 + 0.9530899229693319 +``` +""" # default Gaussian quadrature method for bounded domain function Gauss_Legendre(lb::Number, ub::Number, num_supports::Int; kwargs...)::Tuple (supports, coeffs) = FastGaussQuadrature.gausslegendre(num_supports) @@ -120,6 +209,27 @@ function Gauss_Legendre(lb::Number, ub::Number, num_supports::Int; kwargs...)::T return (supports, coeffs) end +""" + Gauss_Hermite(lb::Number, ub::Number, num_supports::Int; kwargs...)::Tuple + +Return a tuple that contains supports and coefficients generated using +Gauss-Hermite quadrature method. This is useful for univariate parameter in an +infinite interval. + +**Example** +```jldoctest; setup = :(using InfiniteOpt) +julia> (supps, coeffs) = Gauss_Hermite(-Inf, Inf, 5) +([-2.0201828704560856, -0.9585724646138196, -8.881784197001252e-16, 0.9585724646138196, 2.0201828704560856], [1.1814886255359844, 0.986580996751429, 0.9453087204829428, 0.986580996751429, 1.1814886255359844]) + +julia> supps +5-element Array{Float64,1}: + -2.0201828704560856 + -0.9585724646138196 + -8.881784197001252e-16 + 0.9585724646138196 + 2.0201828704560856 +``` +""" # default Gaussian quadrature method for infinite domain function Gauss_Hermite(lb::Number, ub::Number, num_supports::Int; kwargs...)::Tuple if lb != -Inf || ub != Inf @@ -131,6 +241,27 @@ function Gauss_Hermite(lb::Number, ub::Number, num_supports::Int; kwargs...)::Tu return (supports, coeffs) end +""" + Gauss_Laguerre(lb::Number, ub::Number, num_supports::Int; kwargs...)::Tuple + +Return a tuple that contains supports and coefficients generated using +Gauss-Laguerre quadrature method. This is useful for univariate parameter in a +semi-infinite interval. + +**Example** +```jldoctest; setup = :(using InfiniteOpt) +julia> (supps, coeffs) = Gauss_Laguerre(-Inf, 0., 5) +([-0.2635603197181408, -1.413403059106515, -3.596425771040715, -7.08581000585883, -12.640800844275773], [0.6790940422077494, 1.638487873602747, 2.7694432423708255, 4.3156569009208585, 7.219186354354335]) + +julia> supps +5-element Array{Float64,1}: + -0.2635603197181408 + -1.413403059106515 + -3.596425771040715 + -7.08581000585883 + -12.640800844275773 +``` +""" # default Gaussian quadrature method for semi-infinite domain function Gauss_Laguerre(lb::Number, ub::Number, num_supports::Int; kwargs...)::Tuple if ub == Inf @@ -151,6 +282,31 @@ function Gauss_Laguerre(lb::Number, ub::Number, num_supports::Int; kwargs...)::T end end +""" + infinite_transform(lb::Number, ub::Number, num_supports::Int; + [sub_method::Function = MC_sampling, + transform_x::Function = _default_x, + transform_dx::Function = _default_dx, + t_lb::Number = -convert(Number, lb == -Inf && ub == Inf), + t_ub::Number = 1., kwargs...])::Tuple + +Returns a tuple that contains supports and coefficients generated for a +parameter in an infinite or semi-infinite interval. It works by transforming +the original unbounded interval to a finite interval, on which a support +generation method for finite intervals is applied. Then, the generated supports +are transformed back to the original interval. The user is allowed to specify +the support generation method for finite intevals to use, as well as the +transform function. The default transform function is +``t \\in [-\\infty, \\infty] \\rightarrow x \\in [-1, 1]: t(x) = \\frac{t}{1-t^2}`` +``t \\in [a, \\infty] \\rightarrow x \\in [0, 1]: t(x) = a + \\frac{t}{1-t}`` +``t \\in [-\\infty, a] \\rightarrow x \\in [0, 1]: t(x) = a - \\frac{1-t}{t}`` + +**Example** +```jldoctest; setup = :(using InfiniteOpt) +julia> (supps, coeffs) = infinite_transform(-Inf, Inf, 5, sub_method = Gauss_Legendre) +([-5.06704059565454, -0.7583532171678754, 0.0, 0.7583532171678754, 5.06704059565454], [13.490960583398396, 1.2245949721571516, 0.5688888888888889, 1.2245949721571516, 13.490960583398396]) +``` +""" # transform (semi-)infinite domain to finite domain function infinite_transform(lb::Number, ub::Number, num_supports::Int; sub_method::Function = MC_sampling, diff --git a/src/measures.jl b/src/measures.jl index cb17e3b5f..bc711c650 100644 --- a/src/measures.jl +++ b/src/measures.jl @@ -410,25 +410,45 @@ end """ measure(expr::JuMP.AbstractJuMPScalar, - params::Union{ParameterRef, Vector{ParameterRef}}, - lb::Union{Number, Vector{Number}}, - ub::Union{Number, Vector{Number}}; - eval_method::Function, num_supports::Int, weight_func::Function, - use_existing_supports::Bool = false)::MeasureRef + [params::Union{ParameterRef, AbstractArray{<:ParameterRef}, + Nothing} = nothing, + lb::Union{Number, AbstractArray{<:Number}, Nothing}, + ub::Union{Number, AbstractArray{<:Number}, Nothing}]; + [eval_method::Function = MC_sampling, num_supports::Int = 50, + weight_func::Function = _w, name = "", + use_existing_supports::Bool = false, + call_from_expect::Bool = false])::MeasureRef Returns a measure reference that evaluates `expr` without using an object of [`AbstractMeasureData`](@ref) type. Similar to the main [`measure`](@ref) method, this function aims to implement measures of the form: ``\\int_{p \\in P} expr(p) w(p) dp`` where ``p`` is an infinite parameter (scalar -or vector) and ``w`` is the weight function. This function will serve as a -flexible interface where users only have to provide necessary data about the +or vector) and ``w`` is the weight function. This function serves as a flexible +interface where users only have to provide necessary data about the integration. Instead of taking an [`AbstractMeasureData`](@ref) object as input, this function constructs the [`AbstractMeasureData`](@ref) object using some -default numerical integration schemes. +default numerical integration schemes. By default, the function will generate +points by Monte Carlo sampling from the interval if the parameter is in an +[`IntervalSet`](@ref), or from the underlying distribution if the parameter is +in a [`DistributionSet`](@ref). If the expression involves multiple groups of +parameters, the user needs to specify the parameter. The user can also specify +lower bounds and upper bounds for the parameters, number of supports to +generate, and the function to generate the supports with. The last option makes +the method extendible for more schemes to generate supports. **Example** -```julia +```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(seed = true)) +julia> @infinite_parameter(model, x in [0., 1.]) +x + +julia> @infinite_variable(model, f(x)) +f(x) + +julia> meas = measure(f, num_supports = 5) +(f(x)) +julia> expand(meas) +0.2 f(0.8236475079774124) + 0.2 f(0.9103565379264364) + 0.2 f(0.16456579813368521) + 0.2 f(0.17732884646626457) + 0.2 f(0.278880109331201) ``` """ # Measure function that takes non-AbstractMeasureData types @@ -550,6 +570,32 @@ function measure(expr::JuMP.AbstractJuMPScalar, return measure(expr, data) end +""" + expect(expr::JuMP.AbstractJuMPScalar, + [params::Union{ParameterRef, AbstractArray{<:ParameterRef}, + Nothing} = nothing]; + [num_supports::Int = 50, + use_existing_supports::Bool = false])::MeasureRef + +Creates a measure that represents the expected value of an expression in +a random parameter involved in the expression. Return the [`MeasureRef`](@ref) +of the created measure. + +**Example** +```jldoctest; setup = :(using InfiniteOpt, JuMP, Distributions; model = InfiniteModel(seed = true)) +julia> @infinite_parameter(model, x in Normal(0., 1.)) +x + +julia> @infinite_variable(model, f(x)) +f(x) + +julia> meas = expect(f, num_supports = 2) +expect(f(x)) + +julia> expand(meas) +0.5 f(0.6791074260357777) + 0.5 f(0.8284134829000359) +``` +""" # expectation measure function expect(expr::JuMP.AbstractJuMPScalar, params::Union{ParameterRef, AbstractArray{<:ParameterRef}, @@ -561,6 +607,30 @@ function expect(expr::JuMP.AbstractJuMPScalar, call_from_expect = true) end +""" + support_sum(expr::JuMP.AbstractJuMPScalar, + [params::Union{ParameterRef, AbstractArray{<:ParameterRef}, + Nothing} = nothing])::MeasureRef + +Creates a measure that represents the sum of the expression over a parameter +using its existing supports. Return the [`MeasureRef`](@ref) of the created +measure. + +**Example** +```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel()) +julia> @infinite_parameter(model, x in [0, 1], supports = [0.3, 0.7]) +x + +julia> @infinite_variable(model, f(x)) +f(x) + +julia> meas = support_sum(f) +sum(f(x)) + +julia> expand(meas) +f(0.3) + f(0.7) +``` +""" # sum measure function support_sum(expr::JuMP.AbstractJuMPScalar, params::Union{ParameterRef, AbstractArray{<:ParameterRef}, Nothing} From 40280b37e2014b3e9c6f6fe581dc86421e4ff342 Mon Sep 17 00:00:00 2001 From: Weiqi Date: Tue, 17 Dec 2019 13:39:41 -0600 Subject: [PATCH 09/10] fix tests --- test/runtests.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index d510a2527..cf5c6ba80 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,7 +17,6 @@ include("utilities.jl") println("-----------------------------------------") println("----------------UNIT TESTS---------------") println("-----------------------------------------") -#= @time @testset "Datatypes" begin include("datatypes.jl") end println("") @time @testset "Utilities" begin include("utility_tests.jl") end @@ -39,14 +38,12 @@ println("") println("") @time @testset "Expression Methods" begin include("expressions.jl") end println("") -=# @time @testset "Measure Evaluation Methods" begin include("MeasureEvalMethods/methods.jl") end println("") @time @testset "Measure Methods" begin include("measures.jl") end println("") -#= @time @testset "Objective Methods" begin include("objective.jl") end println("") @time @testset "Constraint Methods" begin include("constraints.jl") end @@ -71,7 +68,6 @@ println("") # TODO add extension tests # TODO add involved deletion tests -=# println("-----------------------------------------") println("-------------TESTING COMPLETE------------") println("-----------------------------------------") From aee5f4de8d8aec0c2d50ed9cdab140b3038a1af3 Mon Sep 17 00:00:00 2001 From: Weiqi Date: Tue, 17 Dec 2019 15:18:22 -0600 Subject: [PATCH 10/10] objective.md fixes --- docs/src/guide/objective.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/guide/objective.md b/docs/src/guide/objective.md index e163a1828..9bd520848 100644 --- a/docs/src/guide/objective.md +++ b/docs/src/guide/objective.md @@ -66,7 +66,7 @@ julia> objective_function(model) 0.5 x[1] + 0.5 x[2] + expect(y(ξ)² - y(ξ)) julia> objective_function_type(model) -GenericQuadExpr{Float64,MeasureFiniteVariableRef} +GenericAffExpr{Float64,MeasureFiniteVariableRef} ``` The objective sense can be one of three possibilities: `MIN_SENSE`, `MAX_SENSE`, or `FEASIBILITY_SENSE`. The later sense applies to models that contain no @@ -94,7 +94,7 @@ macro. The objective function is specified via [`set_objective_function`](@ref J simply updates the expression stored in the objective. For example, let's update out objective to simply be ``0.5x_1 + 0.5x_2``: ```jldoctest obj -julia> set_objective_function(model, 0.5 x[1] + 0.5 x[2]) +julia> set_objective_function(model, 0.5x[1] + 0.5x[2]) julia> objective_function(model) 0.5 x[1] + 0.5 x[2]