function runTSEMOptimizationpm(designFunction,initdesignfile,ub,lb,no_outputs,init_dataset_size,maxeval,niter,maxNumberOfProcesses,save2filename,y_refdesign,labelObjective,Userdata)
      
% designFunction - function to process generated designs
% initdesignfile - initial mxa-file (base design)
% ub and lb - lower and upper bounds on design variables
% no_outputs - number of objectives
% init_dataset_size - number of initial designs generated using latin hypercube sampling
% maxeval - total number of generated designs
% niter - number of optimization algorithm iterations
% maxNumberOfProcesses - maximum number of simulations run in parallel by designFunction
% save2filename - save optimization algorithm results to this mat-file
% y_refdesign - (optional) reference design objectives to be displayer on the Pareto plot as a red "x"
% labelObjective - objective names

no_inputs  = length(ub);                                                   % number of decision variables
opt = TSEMO_options;                                                       % call options for solver, see TSEMO_options file to adjust
opt.maxeval = maxeval-init_dataset_size;                                   % number of function evaluations before termination
opt.NoOfBachSequential = ceil((maxeval-init_dataset_size)/niter);          % number of function evaluations per iteration
opt.no_outputs = no_outputs;                                               % number of objectives
% Total number of iterations = opt.maxeval/opt.NoOfBachSequential
ContinueFromIt = 0;                   % by default, continue from generating initial dataset
if exist('save2filename','var') && ~isempty(save2filename) && exist(save2filename,'file')
    button = questdlg({['File ' save2filename]...
        'already exists. Do you want to continue the optimization from where it was interrupted?'...
        'Yes – continue from where it left off'...
        'No – start the optimization from the beginning'}...
        ,'Optimization data file exists','Yes','No','Yes');
    if strcmp(button,'Yes')
        try
            load(save2filename,'X','Y','Ypareto','motorProps','settingsMagnetostatic','Opt','Userdata');
            opt = Opt;
            opt.maxeval = maxeval-init_dataset_size;    % number of function evaluations before termination
            nY = size(Y,1);
        catch
            error([save2filename ' is not correct optimization data file']);
        end
        if nY<init_dataset_size
            ContinueFromIt = 0;     % continue from generating initial dataset
        elseif nY>=maxeval
            disp('Optimization complete');
            return
        else
            ContinueFromIt = ceil((nY-init_dataset_size)/opt.NoOfBachSequential)+1;
        end
    end
end
if ContinueFromIt == 0     % continue from generating initial dataset
    % Generate initial dataset
    X = lhsu(zeros(1,no_inputs),ones(1,no_inputs),init_dataset_size);     % Latin hypercube design
    for k = 1:size(X,1)
        X(k,:) = X(k,:).*(ub-lb)+lb;        % adjustment of bounds
    end
    [motorProps,settingsMagnetostatic] = openMXP(initdesignfile);
    msg0 = ['(Iteration 0 of ' num2str(niter) ':'];
    [Y, ~, Ypareto, Userdata] = designFunction(motorProps,settingsMagnetostatic,initdesignfile,maxNumberOfProcesses,X,msg0,Userdata);

    if exist('save2filename','var') && ~isempty(save2filename)
        Opt = opt;
        save(save2filename);   % save optimization algorithm results to mat-file
    end
    ContinueFromIt = 1;   % start optimization from first iteration
end


% Start algorithm to find Pareto front
[Xpareto,Ypareto,X,Y,XparetoGP,YparetoGP,YparetoGPstd,hypf] = TSEMO_V4(designFunction,X,Y,Ypareto,lb,ub,opt,motorProps,settingsMagnetostatic,initdesignfile,maxNumberOfProcesses,ContinueFromIt,save2filename,Userdata);

% INPUTS
% designFunction denotes the function to be optimized
% X and Y   are the initial datasets to create a surrogate model
% lb and ub are the lower and upper bound of the decision variables
% opt is the option structure of the algorithm

% OUTPUTS
%   Xpareto and Ypareto correspond to the current best Pareto set and Pareto
%   front respectively. 
%   X and Y are the complete dataset of the decision variables and 
%   the objectives respectively. 
%   XparetoGP and YparetoGP represent the Pareto set and Pareto front of the 
%   final Gaussian process model within the algorithm. It is recommended to
%   use these as final result for problems with measurement noise. 
%   YparetoGPstd denotes the standard deviations of the predictions of
%   the GP pareto front YparetoGP  
%   hypf represents the final hyperparameters found for analysis

% For each iteration the current iteration number is displayed, the
% predicted hypervolume improvement and the time taken.

% TS-EMO creates a log file named "TSEMO_log.txt" that contains all relevant information 
% over the entire algorithm run. 
if exist('save2filename','var') && ~isempty(save2filename)
    save(save2filename);   % save optimization algorithm results to mat-file
end

% Visualise results
if ~exist('y_refdesign','var')
    y_refdesign=[];
end
if ~exist('labelObjective','var')
    labelObjective=[];
end
plotPareto(X,Ypareto,init_dataset_size,initdesignfile,motorProps,settingsMagnetostatic,designFunction,y_refdesign,labelObjective,[],[],[],Userdata);





