From 62d69251cd431ac6bb9a8a3c978511776337fa80 Mon Sep 17 00:00:00 2001 From: LudovicoBessi <21291898+ludoro@users.noreply.github.com> Date: Mon, 6 Jul 2020 17:49:15 +0200 Subject: [PATCH 1/3] Update README.md --- README.md | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 63860bfc7..adeadc4a6 100644 --- a/README.md +++ b/README.md @@ -29,7 +29,8 @@ The construction of a surrogate model can be seen as a three-step process: ## ALL the currently available surrogate models: - Kriging -- Radial Basis Function +- Radial Basis +- Wendland - Linear - Second Order Polynomial - Support Vector Machines (SVM) @@ -37,6 +38,10 @@ The construction of a surrogate model can be seen as a three-step process: - Random Forests - Lobachesky - Inverse-distance +- Polynomial expansions +- Variable fidelity +- Mixture of experts + ## ALL the currently available optimization methods: From 9af10dab62acdbf9a66a9be8f9a6244417f73fbb Mon Sep 17 00:00:00 2001 From: LudovicoBessi <21291898+ludoro@users.noreply.github.com> Date: Mon, 6 Jul 2020 18:06:34 +0200 Subject: [PATCH 2/3] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index d7108a1d4..39ccb713a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Surrogates" uuid = "6fc51010-71bc-11e9-0e15-a3fcc6593c49" authors = ["ludoro "] -version = "1.1.2" +version = "1.2" [deps] Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5" From df2a98b621bf12362d7a02cfdf9fbeb562fef851 Mon Sep 17 00:00:00 2001 From: Andrea Cognolato Date: Tue, 7 Jul 2020 14:05:47 +0200 Subject: [PATCH 3/3] Test implementation of multivariable polynomial basis (#188) * Test implementation of multivariable polynomial basis * Implement multivariate polynomial basis * Make multivar_poly_basis zygote-compatible * Add Zygote as build dependency --- Project.toml | 1 + src/Radials.jl | 96 +++++++++++++++++++++++++++++++++++--------------- 2 files changed, 68 insertions(+), 29 deletions(-) diff --git a/Project.toml b/Project.toml index 39ccb713a..a644f0638 100644 --- a/Project.toml +++ b/Project.toml @@ -19,6 +19,7 @@ Sobol = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Stheno = "8188c328-b5d6-583d-959b-9690869a5511" XGBoost = "009559a3-9522-5dbb-924b-0b6ed2b22bb9" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [compat] Clustering = "0.14" diff --git a/src/Radials.jl b/src/Radials.jl index d50e4dd53..ffed21930 100644 --- a/src/Radials.jl +++ b/src/Radials.jl @@ -14,13 +14,15 @@ mutable struct RadialBasis{F,Q,X,Y,L,U,C,S,D} <: AbstractSurrogate end mutable struct RadialFunction{Q,P} - q::Q + q::Q # degree of polynomial phi::P end linearRadial = RadialFunction(0,z->norm(z)) + cubicRadial = RadialFunction(1,z->norm(z)^3) multiquadricRadial = RadialFunction(1,z->sqrt(norm(z)^2+1)) + thinplateRadial = RadialFunction(2, z->begin result = norm(z)^2 * log(norm(z)) ifelse(iszero(z), zero(result), result) @@ -32,7 +34,7 @@ end) Constructor for RadialBasis surrogate. """ function RadialBasis(x, y, lb::Number, ub::Number; rad::RadialFunction=linearRadial, scale_factor::Real=1.0, sparse = false) - q = rad.q + 1 + q = rad.q phi = rad.phi coeff = _calc_coeffs(x, y, lb, ub, phi, q,scale_factor, sparse) return RadialBasis(phi, q, x, y, lb, ub, coeff,scale_factor,sparse) @@ -44,22 +46,29 @@ RadialBasis(x,y,lb,ub,rad::RadialFunction, scale_factor::Float = 1.0) Constructor for RadialBasis surrogate """ function RadialBasis(x, y, lb, ub; rad::RadialFunction = linearRadial, scale_factor::Real=1.0, sparse = false) - q = rad.q + 1 + q = rad.q phi = rad.phi coeff = _calc_coeffs(x, y, lb, ub, phi, q, scale_factor, sparse) return RadialBasis(phi, q, x, y, lb, ub, coeff,scale_factor, sparse) end function _calc_coeffs(x, y, lb, ub, phi, q, scale_factor, sparse) + nd = length(first(x)) + num_poly_terms = binomial(q + nd, q) + D = _construct_rbf_interp_matrix(x, first(x), lb, ub, phi, q, scale_factor, sparse) - Y = _construct_rbf_y_matrix(y, first(y), length(y) + q) + Y = _construct_rbf_y_matrix(y, first(y), length(y) + num_poly_terms) + coeff = D \ Y return coeff end function _construct_rbf_interp_matrix(x, x_el::Number, lb, ub, phi, q, scale_factor, sparse) n = length(x) - m = n + q + + num_poly_terms = binomial(q + 1, q) + m = n + num_poly_terms + if sparse D = ExtendableSparseMatrix{eltype(x_el),Int}(m,m) else @@ -70,7 +79,7 @@ function _construct_rbf_interp_matrix(x, x_el::Number, lb, ub, phi, q, scale_fac D[i,j] = phi( (x[i] .- x[j]) ./ scale_factor ) end if i <= n - for k = 1:q + for k = 1:num_poly_terms D[i,n+k] = _scaled_chebyshev(x[i], k-1, lb, ub) end end @@ -82,10 +91,10 @@ end function _construct_rbf_interp_matrix(x, x_el, lb, ub, phi, q, scale_factor,sparse) n = length(x) nd = length(x_el) - central_point = _center_bounds(x_el, lb, ub) - sum_half_diameter = sum((ub[k]-lb[k])/2 for k = 1:nd) - mean_half_diameter = sum_half_diameter / nd - m = n+q + + num_poly_terms = binomial(q + nd, q) + m = n + num_poly_terms + if sparse D = ExtendableSparseMatrix{eltype(x_el),Int}(m,m) else @@ -96,8 +105,8 @@ function _construct_rbf_interp_matrix(x, x_el, lb, ub, phi, q, scale_factor,spar D[i,j] = phi( (x[i] .- x[j]) ./ scale_factor) end if i < n + 1 - for k = 1:q - D[i,n+k] = centralized_monomial(x[i], k-1, mean_half_diameter, central_point) + for k = 1:num_poly_terms + D[i,n+k] = multivar_poly_basis(x[i], k-1, nd, q) end end end @@ -108,22 +117,49 @@ end _construct_rbf_y_matrix(y, y_el::Number, m) = [i <= length(y) ? y[i] : zero(y_el) for i = 1:m] _construct_rbf_y_matrix(y, y_el, m) = [i <= length(y) ? y[i][j] : zero(first(y_el)) for i=1:m, j=1:length(y_el)] -""" - centralized_monomial(vect,alpha,mean_half_diameter,central_point) - -Returns the value at point vect[] of the alpha degree monomial centralized. +using Zygote: @nograd + +function _make_combination(n, d, ix) + exponents_combinations = [ + e + for e + in collect( + Iterators.product( + Iterators.repeated(0:n, d)... + ) + )[:] + if sum(e) <= n + ] + + return exponents_combinations[ix + 1] +end +# TODO: Is this correct? Do we ever want to differentiate w.r.t n, d, or ix? +# By using @nograd we force the gradient to be 1 for n, d, ix +@nograd _make_combination -#Arguments: --'vect': vector of points i.e [x,y,...,w] --'alpha': degree --'mean_half_diameter': half diameter of the domain --'central_point': central point in the domain """ -function centralized_monomial(vect,alpha,mean_half_diameter,central_point) - if iszero(alpha) return one(eltype(vect)) end - centralized_product = prod(vect .- central_point) - return (centralized_product / mean_half_diameter)^alpha -end + multivar_poly_basis(x, ix, d, n) + +Evaluates in `x` the `ix`-th element of the multivariate polynomial basis of maximum +degree `n` and `d` dimensions. + +Time complexity: `(n+1)^d.` + +# Example +For n=2, d=2 the multivariate polynomial basis is +```` +1, +x,y +x^2,y^2,xy +```` +Therefore the 3rd (ix=3) element is `y` . +Therefore when x=(13,43) and ix=3 this function will return 43. +""" +multivar_poly_basis(x, ix, d, n) = prod( + a^d + for (a, d) + in zip(x, _make_combination(n, d, ix)) +) """ Calculates current estimate of value 'val' with respect to the RadialBasis object. @@ -136,13 +172,14 @@ end function _approx_rbf(val::Number, rad) n = length(rad.x) q = rad.dim_poly + num_poly_terms = binomial(q + 1, q) lb = rad.lb ub = rad.ub approx = zero(rad.coeff[1, :]) for i = 1:n approx += rad.coeff[i, :] * rad.phi( (val .- rad.x[i]) / rad.scale_factor) end - for k = 1:q + for k = 1:num_poly_terms approx += rad.coeff[n+k, :] * _scaled_chebyshev(val, k-1, lb, ub) end return approx @@ -151,6 +188,7 @@ function _approx_rbf(val, rad) n = length(rad.x) d = length(rad.x[1]) q = rad.dim_poly + num_poly_terms = binomial(q + d, q) lb = rad.lb ub = rad.ub sum_half_diameter = sum((ub[k]-lb[k])/2 for k = 1:d) @@ -159,8 +197,8 @@ function _approx_rbf(val, rad) approx = zero(rad.coeff[1, :]) @views approx += sum( rad.coeff[i, :] * rad.phi( (val .- rad.x[i]) ./rad.scale_factor) for i = 1:n) - for k = 1:q - @views approx += rad.coeff[n+k, :] .* centralized_monomial(val, k-1, mean_half_diameter, central_point) + for k = 1:num_poly_terms + @views approx += rad.coeff[n+k, :] .* multivar_poly_basis(val, k-1, d, q) end return approx end