client = jcmoptimizer.Client(); 

%Amplitude of a driven harmonic oscillator. 
function a = amplitudes(omegas, F, omega0, gamma)
    a = F./sqrt((2*omegas*omega0*gamma).^2 + (omega0^2 - omegas.^2).^2);
end

% Definition of the search domain to tune the oscillator
design_space = { ...
    struct('name', 'F', 'type', 'continuous', 'domain', [0.1, 40.0]), ... 
    struct('name', 'omega0', 'type', 'continuous', 'domain', [1.0, 4.0]), ...
    struct('name', 'gamma', 'type', 'continuous', 'domain', [0.03, 0.3]) ...
};

%The amplitude is scaned for 10 different driving frequencies
omegas = linspace(1,3,10);

 % Creation of the study object with study_id 'harmonic_oscillator_fit'
study = client.create_study( ...
    'design_space', design_space, ...
    'driver','ActiveLearning',...
    'study_name','Optimization of resonant system based on Gaussian fit',...
    'study_id', 'harmonic_oscillator_fit');

%definition of variables

%The amplitudes are fitted to a Gaussian + linear expression with
%amplitude A, resonance frequency tau, linear gradient B, and constant offset C.
%The output are the fitted values and the mean-squared error (MSE).
fit = struct();
fit.type = "Fit";
fit.name = "fit";
fit.input = "amplitudes";
fit.expression = "A*exp(-0.5*(omega-tau)^2/sigma^2) + B*omega + C";
fit.output_names = {"A", "B", "C", "tau", "sigma", "MSE"};
fit.output_dim = 6;
fit.model_variables = {"omega"};
fit.variable_values = omegas.';
fit.initial_parameters = [1.0, 0.0, 0.0, 1.0, 0.5];
fit.prior_uncertainties = [10.0, 10.0, 10.0, 10.0, 10.0];

%The expression that defines the loss function:
expression = struct();
expression.type = "Expression";
expression.name = "loss";
expression.expression = "(A-20)^2 + (tau - 2)^2 + sigma^2 + 0.001*MSE";
 
%configuration of study
study.configure( ...
    'max_iter', 30, ...
    'surrogates', { ...
        ...A single-output Gaussian process learns the dependence of the amplitudes
        ...for all wavelengths on the design parameters.
        struct('type', "GP", 'name', "amplitudes", 'output_dim', length(omegas)) ...
    }, ...
    'variables', {fit, expression}, ...
    'objectives', { ...
        ...The only objective of the study is to minimize the loss.
        struct("type", "Minimizer", "variable", "loss"), ...
    }, ...
    'acquisition_optimizer', struct( ...
        ...Sample computation based on fits is more expensive. In this case
        ...advanced sample computation is usually too laborious.
        'compute_suggestion_in_advance',  false ...
    ) ...
);

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

    pause(2); % make objective expensive
    observation = study.new_observation();
    %The harmonic-oscillator amplitude is evaluated for all omega values and the
    %observed values are used as training input for the Gaussian process "amplitudes"
    observation.add(amplitudes(omegas, sample.F, sample.omega0, sample.gamma));
    
end  

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


% Print result
best = study.get_state("driver.best_sample");
fprintf("Best design parameters: F=%f, omega0=%f, gamma=%f\n", ...
        best.F, best.omega0, best.gamma);
fprintf(['Objective: (A-20)^2 + (tau - 2)^2 + sigma^2 + 0.001*MSE = ', ...
         '%.3f\n'], study.get_state('driver.min_objective'));
client.shutdown_server();