Merge pull request #279 from marcoq/marcoq-SectionSampler

Section sampler

docs/src/samples.md

@@ -34,6 +34,12 @@ LowDiscrepancySample{T}
 sample(n,lb,ub,S::LowDiscrepancySample)
 ```
 
+* Sample on section
+```@docs
+SectionSample
+sample(n,lb,ub,S::SectionSample)
+```
+
 ## Adding a new sampling method
 
 Adding a new sampling method is a two- step process:

src/NeuralSurrogate.jl

@@ -1,5 +1,6 @@
 using Flux
 using Flux: @epochs
+
 mutable struct NeuralSurrogate{X,Y,M,L,O,P,N,A,U} <: AbstractSurrogate
     x::X
     y::Y

src/Optimization.jl

@@ -697,11 +697,11 @@ function surrogate_optimize(obj::Function,::DYCORS,lb::Number,ub::Number,surr1::
             while new_points[i] < lb || new_points[i] > ub
                 if new_points[i] > ub
                     #reflection
-                    new_points[i] = maximum(surr1.x) - norm(new_points[i] - maximum(surr1.x))
+                    new_points[i] = max(lb, maximum(surr1.x) - norm(new_points[i] - maximum(surr1.x)))
                 end
                 if new_points[i] < lb
                     #reflection
-                    new_points[i] = minimum(surr1.x) + norm(new_points[i]-minimum(surr1.x))
+                    new_points[i] = min(ub, minimum(surr1.x) + norm(new_points[i]-minimum(surr1.x)))
                 end
             end
         end
@@ -831,10 +831,10 @@ function surrogate_optimize(obj::Function,::DYCORS,lb,ub,surrn::AbstractSurrogat
             for j = 1:d
                 while new_points[i,j] < lb[j] || new_points[i,j] > ub[j]
                     if new_points[i,j] > ub[j]
-                        new_points[i,j] = maximum(surrn.x)[j] - norm(new_points[i,j] - maximum(surrn.x)[j])
+                        new_points[i,j] = max(lb[j], maximum(surrn.x)[j] - norm(new_points[i,j] - maximum(surrn.x)[j]))
                     end
                     if new_points[i,j] < lb[j]
-                        new_points[i,j] = minimum(surrn.x)[j] + norm(new_points[i]-minimum(surrn.x)[j])
+                        new_points[i,j] = min(ub[j], minimum(surrn.x)[j] + norm(new_points[i]-minimum(surrn.x)[j]))
                     end
                 end
             end
@@ -1711,3 +1711,82 @@ end
      end
      return pareto_set,pareto_front
  end
+
+
+function surrogate_optimize(obj::Function,::EI,lb,ub,krig,sample_type::SectionSample;maxiters=100,num_new_samples=100)
+    dtol = 1e-3*norm(ub-lb)
+    eps = 0.01
+    for i = 1:maxiters
+        d = length(krig.x)
+        new_sample = sample(num_new_samples,lb,ub,sample_type)
+        f_max = maximum(krig.y)
+        evaluations = zeros(eltype(krig.x[1]),num_new_samples)
+        point_found = false
+        new_x_max = zero(eltype(krig.x[1]))
+        new_y_max = zero(eltype(krig.x[1]))
+        diff_x = zeros(eltype(krig.x[1]),d)
+        while point_found == false
+            for j = 1:length(new_sample)
+                std = std_error_at_point(krig,new_sample[j])
+                u = krig(new_sample[j])
+                if abs(std) > 1e-6
+                    z = (u - f_max - eps)/std
+                else
+                    z = 0
+                end
+                evaluations[j] = (u-f_max-eps)*cdf(Normal(),z) + std*pdf(Normal(),z)
+            end
+            index_max = argmax(evaluations)
+            x_new = new_sample[index_max]
+            y_new = maximum(evaluations)
+            for l = 1:d
+                diff_x[l] = norm(krig.x[l] .- x_new)
+            end
+            bit_x = diff_x .> dtol
+            #new_min_x has to have some distance from krig.x
+            if false in bit_x
+                #The new_point is not actually that new, discard it!
+                deleteat!(evaluations,index_max)
+                deleteat!(new_sample,index_max)
+                if length(new_sample) == 0
+                    println("Out of sampling points")
+                    return section_sampler_returner(sample_type, krig.x, krig.y, lb, ub, krig)
+                end
+             else
+                point_found = true
+                new_x_max = x_new
+                new_y_max = y_new
+            end
+        end
+        if new_y_max < 1e-6*norm(maximum(krig.y)-minimum(krig.y))
+            return section_sampler_returner(sample_type, krig.x, krig.y, lb, ub, krig)
+        end
+        add_point!(krig,Tuple(new_x_max),obj(new_x_max))
+    end
+    println("Completed maximum number of iterations")
+end
+
+function section_sampler_returner(
+        sample_type::SectionSample, surrn_x, surrn_y, 
+        lb, ub, surrn)
+    d_fixed = Surrogates.fixed_dimensions(sample_type)
+    @assert length(surrn_y) == size(surrn_x)[1]
+    surrn_xy = [(surrn_x[y], surrn_y[y]) for y in 1:length(surrn_y)]
+    section_surr1_xy = filter(
+        xyz->xyz[1][d_fixed]==Tuple(sample_type.x0[d_fixed]),
+        surrn_xy)
+    section_surr1_x = [xy[1] for xy in section_surr1_xy]
+    section_surr1_y = [xy[2] for xy in section_surr1_xy]
+    if length(section_surr1_xy) == 0
+        @debug "No new point added - surrogate locally stable"
+        N_NEW_POINTS = 100
+        section_surr1_x = sample(N_NEW_POINTS, lb, ub, sample_type)
+        section_surr1_y = zeros(N_NEW_POINTS)
+        for i in 1:size(section_surr1_x, 1)
+            xi = Tuple([section_surr1_x[i, :]...])[1]
+            section_surr1_y[i] = surrn(xi)
+        end
+    end
+    index = argmin(section_surr1_y)
+    return (section_surr1_x[index, :][1], section_surr1_y[index])
+end

src/Sampling.jl

@@ -35,6 +35,10 @@ end
 
 struct GoldenSample <: SamplingAlgorithm end
 
+struct SectionSample{T} <: SamplingAlgorithm
+    x0::Vector{T}
+    sa::SamplingAlgorithm
+end
 
 """
 sample(n,lb,ub,S::GridSample)
@@ -236,3 +240,41 @@ function sample(n,lb,ub,G::GoldenSample)
         return Tuple.(y)
     end
 end
+
+fixed_dimensions(
+        section_sampler::SectionSample)::Vector{Int64} = findall(
+    x->x == false, isnan.(section_sampler.x0))
+
+free_dimensions(
+        section_sampler::SectionSample)::Vector{Int64} = findall(
+    x->x == true, isnan.(section_sampler.x0))
+
+"""
+sample(n,d,K::SectionSample)
+
+Returns Tuples constrained to a section.
+
+In surrogate-based identification and control, optimization can alternate between unconstrained sampling in the full-dimensional parameter space, and sampling constrained on specific sections (e.g. a planes in a 3D volume),
+
+A SectionSampler allows sampling and optimizing on a subset of 'free' dimensions while keeping 'fixed' ones constrained.
+The sampler is defined as in e.g. 
+
+section_sampler_y_is_10 = SectionSample([NaN64, NaN64, 10.0, 10.0], Surrogates.UniformSample())
+
+where the first argument is a Vector{T} in which numbers are fixed coordinates and `NaN`s correspond to free dimensions, and the second argument is a SamplingAlgorithm which is used to sample in the free dimensions.
+"""
+function sample(n,lb,ub,section_sampler::SectionSample)
+    if lb isa Number
+        return rand(Uniform(lb,ub),n)
+    else
+        d_free = Surrogates.free_dimensions(section_sampler)
+        new_samples = sample(n, lb[d_free], ub[d_free], section_sampler.sa)
+        out_as_vec = repeat(section_sampler.x0', n, 1)
+        for y in 1:size(out_as_vec,1)
+            for xi in 1:length(d_free)
+                out_as_vec[y,xi] = new_samples[y][xi]
+            end
+        end
+        [Tuple(out_as_vec[y,:]) for y in 1:size(out_as_vec,1)]
+    end
+end

src/Surrogates.jl

@@ -95,7 +95,7 @@ export AbstractSurrogate, SamplingAlgorithm
 export Kriging, RadialBasis, add_point!, current_estimate, std_error_at_point
 export linearRadial,cubicRadial,multiquadricRadial,thinplateRadial
 export sample, GridSample, UniformSample, SobolSample, LatinHypercubeSample, LowDiscrepancySample
-export RandomSample, KroneckerSample, GoldenSample
+export RandomSample, KroneckerSample, GoldenSample, SectionSample
 export SRBF,LCBS,EI,DYCORS,SOP,EGO,RTEA,SMB,surrogate_optimize
 export LobachevskySurrogate, lobachevsky_integral, lobachevsky_integrate_dimension
 export LinearSurrogate

test/SectionSampleTests.jl

@@ -0,0 +1,53 @@
+"""
+    Sobel-sample x+y in [0,10]x[0,10],
+    then minimize it on Section([NaN,10.0]),
+    and verify that the minimum is on x,y=(0,10)
+    rather than in (0,0)
+"""
+
+using Surrogates
+using Test
+
+lb = [ 0.0, 0.0, 0.0]
+ub = [10.0,10.0,10.0]
+x = Surrogates.sample(10,lb,ub,LatinHypercubeSample())
+f = x -> x[1]+x[2]+x[3]
+y = f.(x)
+f([0,0,0]) == 0
+
+f_hat = Kriging(x,y,lb,ub)
+
+f_hat([0,0,0])
+
+isapprox(f([0,0,0]), f_hat([0,0,0]))
+
+""" The global minimum is at (0,0) """
+
+(xy_min, f_hat_min) = surrogate_optimize(
+    f,
+    DYCORS(), lb, ub,
+    f_hat,
+    SobolSample())
+
+isapprox(xy_min[1], 0.0, atol=1e-3)
+
+""" The minimum on the (0,10) section is around (0,10) """
+
+section_sampler_z_is_10 = SectionSample(
+    [NaN64, NaN64, 10.0],
+    Surrogates.UniformSample())
+
+@test [3] == Surrogates.fixed_dimensions(section_sampler_z_is_10)
+@test [1,2] == Surrogates.free_dimensions(section_sampler_z_is_10)
+
+Surrogates.sample(5, lb, ub, section_sampler_z_is_10)
+
+(xy_min, f_hat_min) = surrogate_optimize(
+    f,
+    EI(), lb, ub,
+    f_hat,
+    section_sampler_z_is_10, maxiters=1000)
+
+isapprox(xy_min[1],  0.0, atol=0.1)
+isapprox(xy_min[2],  0.0, atol=0.1)
+isapprox(xy_min[3], 10.0, atol=0.1)

test/runtests.jl

@@ -19,3 +19,4 @@ using Surrogates
 @testset "VariableFidelity" begin include("VariableFidelity.jl") end
 @testset "Earth" begin include("earth.jl") end
 @testset "Gradient Enhanced Kriging" begin include("GEK.jl") end
+@testset "Section Samplers" begin include("SectionSampleTests.jl") end