client = jcmoptimizer.Client(); 

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

 % Creation of the study object with study_id 'multi_objective_optimization'
study = client.create_study( ...
    'design_space', design_space, ...
    'driver','ActiveLearning',...
    'study_name','Multi-objective optimization',...
    'study_id', 'multi_objective_optimization');
%Lower and upper reference point for hypervolume definition (see figure)
lower_ref=[0, 0];
upper_ref=[5, 50];
 
study.configure( ...
    'max_iter', 50, ...
    'surrogates', { ...
        ...A multi-output Gaussian process that learns the dependence of
        ...the model on the design parameters.
        struct('type', "GP", 'name', "model_vector", 'output_dim', 2, ...
             'correlate_outputs', false)
    }, ...
    'variables', { ...
        ...Selectors for the values of f1 and f2
        struct('type', "SingleSelector", 'name', "f1", ...
             'input', "model_vector", 'select_by_name', "model_vector0"), ...
        struct('type', "SingleSelector", 'name', "f2", ...
              'input', "model_vector", 'select_by_name', "model_vector1") ...
    }, ...
    'objectives', { ...
        ... Multi-minimization objective for f1, f2
        struct('type', "MultiMinimizer", 'variables', ["f1", "f2"], 'name', "objective", ...
             'lower_reference', lower_ref, 'upper_reference', upper_ref) ...
    } ...
);

% Evaluation of the black-box function for specified design parameters
function observation = evaluate(study, sample)

    pause(2); % make objective expensive
    observation = study.new_observation();
    %determine list of the two objective values
    g = 20 + (sample.x2^2 - 10*cos(2*pi*sample.x2));
    observation.add([sample.x1, g/sample.x1]);
    
end  

% Run the minimization
study.set_evaluator(@evaluate);
study.run(); 


fig = figure('Position', [0, 0, 1000, 500]);

%Plot all observations of (f1, f2)
f1 = study.driver.get_observed_values('object_type', "variable", 'name', "f1");
f2 = study.driver.get_observed_values('object_type', "variable", 'name', "f2");
plot(f1.means, f2.means, "o", 'DisplayName', "Donimated observations");
hold on;
 
%Plot the estimate of the Paretro front, i.e. nondominated observations
pf = study.driver.get_state("pareto_front");
nondom_f1 = pf.f1;
nondom_f2 = pf.f2;
plot(nondom_f1, nondom_f2, "o", 'DisplayName', "Pareto front estimate");
 
%Plot analytic Pareto front
X = linspace(1, 4, 100);
plot(X, 10./X, "k--", 'DisplayName', "Analytic Pareto front");
 
%Plot nondominated hypervolume
plot(lower_ref(1), lower_ref(2), "rx", 'DisplayName', "Lower reference");
plot(upper_ref(1), upper_ref(2), "bx", 'DisplayName', "Upper reference");
X = [lower_ref(1), repelem(nondom_f1, 2), upper_ref(1)];
nondom_lower = repelem(lower_ref(2), length(X));
nondom_f2_aug = [upper_ref(2), nondom_f2];
nondom_upper = repelem(nondom_f2_aug, 2);
patch([X, fliplr(X)], [nondom_lower, fliplr(nondom_upper)], 'b', ...
             'FaceAlpha', 0.2, 'EdgeAlpha', 0.2, ...
             'DisplayName', "Nondominated hypervolume");

xlabel("f1");
ylabel("f2");
legend('Location', "north");
grid();
saveas(fig, "multi_objective_minimization.png");

client.shutdown_server();