Comparison of GEKPLS on Welded Beam Test with Kriging, RBF and PolyChaos for JuliaCon 2022 Presentation

comparison_3_dimensions_juliacon.jl

using Surrogates
using Zygote
using SurrogatesPolyChaos
using SurrogatesAbstractGPs
using AbstractGPs

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

## welded beam tests
function welded_beam(x)
    h = x[1]
    l = x[2]
    t = x[3]
    a = 6000 / (sqrt(2) * h * l)
    b = (6000 * (14 + 0.5 * l) * sqrt(0.25 * (l^2 + (h + t)^2))) /
        (2 * (0.707 * h * l * (l^2 / 12 + 0.25 * (h + t)^2)))
    return (sqrt(a^2 + b^2 + l * a * b)) / (sqrt(0.25 * (l^2 + (h + t)^2)))
end

n_train = 50 #number of training points
lb = [0.125, 5.0, 5.0]
ub = [1.0, 10.0, 10.0]
x_train = sample(n_train,lb,ub,SobolSample())
X_train = vector_of_tuples_to_matrix(x_train)
grads = vector_of_tuples_to_matrix2(gradient.(welded_beam, x_train))
y_train_gekpls = reshape(welded_beam.(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 = welded_beam.(x_test)
n_comp = 2
delta_x = 0.01
extra_points = 5
initial_theta = 0.1
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 = welded_beam.(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: 908.7384065223049
println("rmse_krig: ", rmse(y_krig, y_true)) #rmse_krig: 1070.0688052722955
println("rmse_poly: ", rmse(y_poly, y_true)) #rmse_poly: 469.3247547003671
println("rmse_gekpls: ", rmse(y_gekpls, y_true)) #rmse_gekpls: 65.28877678066098