comparison of GEKPLS with RBF and Kriging

gekpls_benchmarked.jl

using Surrogates
using Zygote
using SurrogatesPolyChaos

function vector_of_tuples_to_matrix(v)
    #helper function to convert training data generated by surrogate sampling into a matrix suitable for GEKPLS
    num_rows = length(v)
    num_cols = length(first(v))
    K = zeros(num_rows, num_cols)
    for row in 1:num_rows
        for col in 1:num_cols
            K[row, col]=v[row][col]
        end
    end
    return K
end

function vector_of_tuples_to_matrix2(v)
    #helper function to convert gradients into matrix form for GEKPLS
    num_rows = length(v)
    num_cols = length(first(first(v)))
    K = zeros(num_rows, num_cols)
    for row in 1:num_rows
        for col in 1:num_cols
            K[row, col] = v[row][1][col]
        end
    end
    return K
end


function rmse(a,b)
    #to calculate root mean squared error
    a = vec(a)
    b = vec(b)
    if(size(a)!=size(b))
        println("error in inputs")
        return
    end
    n = size(a,1)
    return sqrt(sum((a - b).^2)/n)
end

function water_flow(x)
    r_w = x[1]
    r = x[2]
    T_u = x[3]
    H_u = x[4]
    T_l = x[5]
    H_l = x[6]
    L = x[7]
    K_w = x[8]
    log_val = log(r/r_w)
    return (2*pi*T_u*(H_u - H_l))/ ( log_val*(1 + (2*L*T_u/(log_val*r_w^2*K_w)) + T_u/T_l))
end

n_train = 1000 #number of training points
lb = [0.05,100,63070,990,63.1,700,1120,9855]
ub = [0.15,50000,115600,1110,116,820,1680,12045]
x_train = sample(n_train,lb,ub,SobolSample())
X_train = vector_of_tuples_to_matrix(x_train)
grads = vector_of_tuples_to_matrix2(gradient.(water_flow, x_train))
y_train_gekpls = reshape(water_flow.(x_train),(size(x_train,1),1)) #gekpls needs input to be in the form of a matrix
xlimits = hcat(lb, ub)

n_test = 100
x_test = sample(n_test,lb,ub,GoldenSample())
X_test = vector_of_tuples_to_matrix(x_test)
y_true = water_flow.(x_test)
n_comp = 3
delta_x = 0.0001
extra_points = 4
initial_theta = 0.01
g = GEKPLS(X_train, y_train_gekpls, grads, n_comp, delta_x, xlimits, extra_points, initial_theta)
y_gekpls = g(X_test)


### models with kriging, radials and polychaos surrogates
y_train = water_flow.(x_train)
my_rad = RadialBasis(x_train, y_train, lb, ub)
y_rad = my_rad.(x_test)
my_krig = Kriging(x_train, y_train, lb, ub)
y_krig = my_krig.(x_test)
my_poly = SurrogatesPolyChaos.PolynomialChaosSurrogate(x_train, y_train, lb, ub)
y_poly = my_poly.(x_test)


### results
println("rmse_rad: ", rmse(y_rad, y_true)) #rmse_rad: 50.38857855638598
println("rmse_krig: ", rmse(y_krig, y_true)) #rmse_krig: 46.08551997741417
println("rmse_gekpls: ", rmse(y_gekpls, y_true)) #rmse_gekpls: 0.01853324953941199
println("rmse_poly: ", rmse(y_poly, y_true)) #rmse_poly: 1.5181687125431087