Skip to content

Commit

Permalink
prep to test integration
Browse files Browse the repository at this point in the history
Former-commit-id: cd8246a
  • Loading branch information
chriscoey committed Mar 27, 2019
1 parent 530477b commit b3291b1
Show file tree
Hide file tree
Showing 2 changed files with 52 additions and 23 deletions.
28 changes: 14 additions & 14 deletions Manifest.toml
Original file line number Diff line number Diff line change
Expand Up @@ -65,9 +65,9 @@ version = "0.2.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "195a3ffcb8b0762684b6821de18f83a16455c6ea"
git-tree-sha1 = "84aa74986c5b9b898b0d1acaf3258741ee64754f"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "2.0.0"
version = "2.1.0"

[[Conda]]
deps = ["Compat", "JSON", "VersionParsing"]
Expand All @@ -77,9 +77,9 @@ version = "1.2.0"

[[Crayons]]
deps = ["Test"]
git-tree-sha1 = "416737eea5c50ee5a08c588ea73d77d5eebc94e7"
git-tree-sha1 = "f621b8ef51fd2004c7cf157ea47f027fdeac5523"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "3.0.0"
version = "4.0.0"

[[DataFrames]]
deps = ["CategoricalArrays", "CodecZlib", "Compat", "DataStreams", "Dates", "InteractiveUtils", "IteratorInterfaceExtensions", "LinearAlgebra", "Missings", "Printf", "Random", "Reexport", "SortingAlgorithms", "Statistics", "StatsBase", "TableTraits", "Tables", "Test", "TranscodingStreams", "Unicode", "WeakRefStrings"]
Expand Down Expand Up @@ -125,9 +125,9 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Distributions]]
deps = ["Distributed", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "c24e9b6500c037673f0241a2783472b8c3d080c7"
git-tree-sha1 = "dec0ebacfbc3a2126c614ab5e903c9ef063688d0"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.16.4"
version = "0.17.0"

[[DynamicPolynomials]]
deps = ["LinearAlgebra", "MultivariatePolynomials", "Pkg", "Reexport", "Test"]
Expand All @@ -154,8 +154,10 @@ uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
[[GSL]]
deps = ["BinaryProvider", "Compat", "Libdl", "Test"]
git-tree-sha1 = "b001d22c3d3efb302f74ee6b946a19f8e59a86c1"
repo-rev = "v0.4.0"
repo-url = "https://github.com/JuliaMath/GSL.jl.git"
uuid = "92c85e6c-cbff-5e0c-80f7-495c94daaecd"
version = "0.4.0"
version = "0.5.0+"

[[InteractiveUtils]]
deps = ["Markdown"]
Expand Down Expand Up @@ -214,10 +216,8 @@ uuid = "a63ad114-7e13-5084-954f-fe012c677804"
[[MultivariateMoments]]
deps = ["LinearAlgebra", "MultivariatePolynomials", "RowEchelon", "SemialgebraicSets"]
git-tree-sha1 = "d85604b2068aa8806d14bf4a2c08df5a732cad57"
repo-rev = "master"
repo-url = "https://github.com/JuliaAlgebra/MultivariateMoments.jl.git"
uuid = "f4abf1af-0426-5881-a0da-e2f168889b5e"
version = "0.2.0+"
version = "0.2.0"

[[MultivariatePolynomials]]
deps = ["LinearAlgebra", "Test"]
Expand Down Expand Up @@ -245,9 +245,9 @@ version = "0.9.6"

[[Parsers]]
deps = ["Dates", "Mmap", "Test", "WeakRefStrings"]
git-tree-sha1 = "dc0469396c9b221223861558ad89141be38635b4"
git-tree-sha1 = "d4e634c9cab597f0ec8f5a1d62c0787e98602c55"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "0.2.20"
version = "0.2.22"

[[Pkg]]
deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
Expand Down Expand Up @@ -402,9 +402,9 @@ version = "0.5.0"

[[TranscodingStreams]]
deps = ["Pkg", "Random", "Test"]
git-tree-sha1 = "d4718bd4db6b7850af3c3833556c6aedbb8e9904"
git-tree-sha1 = "f42956022d8084539f1d7219f632542b0ea686ce"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.2"
version = "0.9.3"

[[URIParser]]
deps = ["Test", "Unicode"]
Expand Down
47 changes: 38 additions & 9 deletions examples/integration/jump.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
#=
Copyright 2018, Chris Coey and contributors
Copyright 2019, Chris Coey, Lea Kapelevich and contributors
approximate integration (upper bound) of a polynomial over a basic semialgebraic set
adapted from "Approximate volume and integration for basic semialgebraic sets"
Expand All @@ -15,30 +15,59 @@ const MU = HYP.ModelUtilities
import MathOptInterface
const MOI = MathOptInterface
import JuMP
using LinearAlgebra
import Random
import MultivariatePolynomials
import DynamicPolynomials
import SemialgebraicSets
# import Random
using Test

function build_quadrature(
d::Int,
p, # polynomial parameter (in objective of moment problem; must be positive on K)
K_dom::MU.Domain,
B_dom::MU.Domain; # canonical set (for which we have quadrature weights) containing K_dom
B_dom::MU.Domain,
ppar,
)
# generate interpolation for K
(U, pts, P0, PWts, _) = MU.interpolate(K_dom, d, sample = true, calc_w = false)

# get quadrature weights for B
y2 = MU.get_weights(B_dom, pts)
y2 = # TODO

# build JuMP model
model = JuMP.Model(JuMP.with_optimizer(HYP.Optimizer, verbose = true))
JuMP.@variable(model, y1[1:U]) # moments of μ1
JuMP.@objective(model, Max, sum(y1[i] * p(pts[i, :]) for i in 1:U))
JuMP.@objective(model, Max, sum(y1[i] * ppar(pts[i, :]) for i in 1:U))
JuMP.@constraint(model, y1 in HYP.WSOSPolyInterpCone(U, [P0, PWts...], true))
JuMP.@constraint(model, y2 - y1 in HYP.WSOSPolyInterpCone(U, [P0], true))
JuMP.@constraint(model, y2 .- y1 in HYP.WSOSPolyInterpCone(U, [P0], true))

return (model, y1, pts)
end


function integrate_poly(
p, # poly to integrate
d::Int, # degree of moment problem
K_dom::MU.Domain, # domain on which quadrature weights are valid
B_dom::MU.Domain, # canonical set (for which we have quadrature weights) containing K_dom
ppar, # polynomial as the parameter in objective of moment problem
)
# optimize to get quadrature weights
(model, y1, pts) = build_quadrature(d, K_dom, B_dom, ppar)
JuMP.optimize!(model)

w = JuMP.value.(y1)
integral = sum(w[i] * p(pts[i, :]) for i in eachindex(w))

println(w)
println(integral)
end


n = 1
DynamicPolynomials.@polyvar x
d = 10
p = 1.0 + 0.0x
K_dom = MU.Box([0.0], [0.5]) # TODO SemialgebraicSets.@set(x * (1/2 - x) >= 0)
B_dom = MU.Box([-1.0], [1.0])
ppar = 1.0 + 0.0x

integrate_poly(p, d, K_dom, B_dom, ppar)

0 comments on commit b3291b1

Please sign in to comment.