WIP: add example for integral/quadrature-weights on BSS

Manifest.toml

@@ -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"]
@@ -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"]
@@ -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"]
@@ -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"]
@@ -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"]
@@ -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"]
@@ -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"]

examples/integration/jump.jl

@@ -0,0 +1,73 @@
+#=
+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"
+by Henrion, Lasserre, & Savorgnan (2009)
+https://pdfs.semanticscholar.org/893b/e70a990901d7b6b2f052cbcb5883a043b5d9.pdf
+=#
+
+import Hypatia
+const HYP = Hypatia
+const CO = HYP.Cones
+const MU = HYP.ModelUtilities
+
+import MathOptInterface
+const MOI = MathOptInterface
+import JuMP
+import MultivariatePolynomials
+import DynamicPolynomials
+import SemialgebraicSets
+# import Random
+using Test
+
+function build_quadrature(
+    d::Int,
+    K_dom::MU.Domain,
+    B_dom::MU.Domain,
+    ppar,
+    )
+    # generate interpolation for B
+    (U, pts, P0, PWts, w) = MU.interpolate(B_dom, d, sample = true, calc_w = true)
+
+    # get weights for B
+    y2 = # TODO
+
+    # build JuMP model
+    model = JuMP.Model(JuMP.with_optimizer(HYP.Optimizer, verbose = true))
+    JuMP.@variable(model, y[1:U]) # moments of μ1
+    JuMP.@objective(model, Max, sum(y[i] * ppar(pts[i, :]) for i in 1:U))
+    JuMP.@constraint(model, y in HYP.WSOSPolyInterpCone(U, [P0, PWts...], true))
+    JuMP.@constraint(model, w - y in HYP.WSOSPolyInterpCone(U, [P0], true))
+
+    return (model, y, 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, y, pts) = build_quadrature(d, K_dom, B_dom, ppar)
+    JuMP.optimize!(model)
+
+    w = JuMP.value.(y)
+    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)