client = jcmoptimizer.Client(); 

% Definition of the search domain
design_space = {...
   struct('name', 'x1', 'type', 'continuous', 'domain', [-1.5,1.5]), ...
   struct('name', 'x2', 'type', 'continuous', 'domain', [-1.5,1.5]), ...
};

% Definition of fixed environment parameter
environment = {...     
   struct('name', 'radius', 'type', 'fixed', 'domain', 1.5) ...
};

% Definition of a constraint on the search domain
constraints = {...
   struct('name', 'circle', 'expression', 'sqrt(x1^2 + x2^2) <= radius')...
};

%Creation of studies to benchmark against each other
studies = struct();
drivers = {"BayesianOptimization", "ActiveLearning", "CMAES", ...
           "DifferentialEvolution", "ScipyMinimizer"};
for i = 1:length(drivers)
   studies.(drivers{i}) = client.create_study( ...
        'design_space', design_space, ...
        'environment', environment, ...
        'constraints', constraints, ...
        'driver', drivers{i}, ...
        'study_name', drivers{i}, ...
        'study_id', sprintf("benchmark_%s", drivers{i}), ...
        'open_browser' , false ...
   );
end

%Configuration of all studies
config = {'num_parallel', 1, 'max_iter', 250};
min_val = 1e-3; %Stop study when this value was observed
for i = 1:length(drivers)
    driver = drivers{i};
    study = studies.(driver);
    if strcmp(driver, "ActiveLearning")
        study.configure( ...
           'surrogates', {struct('type', "NN")}, ...
           'objectives', {struct('type', "Minimizer", 'min_val', min_val)}, ...
           config{:});
    elseif strcmp(driver, "ScipyMinimizer")
        %For a more global search, 3 initial gradient-free Nelder-Mead are started
        study.configure( ...
           'method', "Nelder-Mead", 'num_initial', 3, ...
           'min_val', min_val, config{:});
    else
        study.configure('min_val', min_val, config{:});
    end
end
% Evaluation of the black-box function for specified design parameters
function observation = evaluate(study, sample)

    pause(2); % make objective expensive
    observation = study.new_observation();
    x1 = sample.x1;
    x2 = sample.x2;
    observation.add(10*2 ...
                    + (x1.^2-10*cos(2*pi*x1)) ...
                    + (x2.^2-10*cos(2*pi*x2)) ...
                   );
end  

% Creation of a benchmark with 6 repetitions and add all 4 studies
benchmark = client.create_benchmark('num_average', 6);
for i = 1:length(drivers)
   benchmark.add_study(studies.(drivers{i}));
end

% Run the benchmark - this will take a while
benchmark.set_evaluator(@evaluate);
benchmark.run();

% Plot cummin convergence w.r.t. number of evaluations and time
fig = figure('Position', [0, 0, 1000, 500]);
subplot(1, 2, 1) 
data = benchmark.get_data('x_type', "num_evaluations");

for i = 1:length(data.names)
  std_error = data.sdev(i,:) ./ sqrt(6);
  p = plot(data.X(i,:), data.Y(i,:), 'LineWidth', 2.0);
  hold on
  patch([data.X(i,:), fliplr(data.X(i,:))], ...
        [data.Y(i,:) - std_error, fliplr(data.Y(i,:) + std_error)], ...
        p(1).Color, 'FaceAlpha', 0.2, 'EdgeAlpha', 0.2);
end
grid()
xlabel("Number of Evaluations", 'FontSize', 12)
ylabel("Average Cummulative Minimum", 'FontSize', 12)
ylim([-0.4, 10])

subplot(1, 2, 2)
data = benchmark.get_data('x_type', "time");
plots=[];
for i = 1:length(data.names)
  std_error = data.sdev(i,:) ./ sqrt(6);
  p = plot(data.X(i,:), data.Y(i,:), 'LineWidth', 2.0, 'DisplayName', data.names{i});
  plots(i) = p;
  hold on
  patch([data.X(i,:), fliplr(data.X(i,:))], ...
        [data.Y(i,:) - std_error, fliplr(data.Y(i,:) + std_error)], ...
        p(1).Color, 'FaceAlpha', 0.2, 'EdgeAlpha', 0.2);
end
legend(plots, data.names{:});
grid()
xlabel("Time (sec)", 'FontSize', 12)
ylim([-0.4, 10])
saveas(fig, "benchmark.png")


client.shutdown_server();