function traj = priv_simulateConstrainedRandom(sys,params,options)
% priv_simulateConstrainedRandom - performs several random simulation of the system
%    so that the simulations stay within a given reachable set; this 
%    reachable set is typically a backwards minmax reachable set. These 
%    reachable sets assume that a control input exists such that the 
%    solutions stay within the reachable set. However, the inputs are yet
%    to be determined by simulateConstrained of the individual system 
%    dynamics classes. It can be set how many simulations should be 
%    performed and what percentage of initial states should start at 
%    vertices of the initial set.
%
% Syntax:
%    res = priv_simulateConstrainedRandom(sys,params,options)
%
% Inputs:
%    sys - contDynamics object
%    params - model parameters
%    options - settings for random simulation
%
% Outputs:
%    traj - object of class trajectory storing time and states of the 
%          simulated trajectories.

% Authors:       Matthias Althoff
% Written:       05-January-2023
% Last update:   29-June-2024 (TL, bug fix for empty R0)
%                08-August-2025 (LL, create trajectory object)
% Last revision: ---

% ------------------------------ BEGIN CODE -------------------------------

% generate random initial points
nrExtreme = ceil(options.points*options.fracVert);
nrStandard = options.points - nrExtreme;
X0 = [];
if nrExtreme > 0
	X0 = [X0, randPoint(params.R0,nrExtreme,'extreme')]; 
end
if nrStandard > 0
	X0 = [X0, randPoint(params.R0,nrStandard,'standard')];
end

% check number of generated points (might be different for e.g. empty sets)
nrPoints = size(X0,2);

% initialize time and state
y_r = [];

% output equation only for linearSys and linearSysDT currently
comp_y = (isa(sys,'linearSys') || isa(sys,'linearSysDT')) && ~isempty(sys.C);

% loop over all starting points in X0
traj(nrPoints,1) = trajectory();
for r = 1:nrPoints
    
    % start of trajectory
    params.x0 = X0(:,r);

    % simulate dynamical system
    if comp_y
        % Is output desired?
        [t_r,x_r,~,y_r] = simulateConstrained(sys,params,options);
    else
        [t_r,x_r] = simulateConstrained(sys,params,options);
    end

    % construct trajectory object storing the simulation results
    traj(r) = trajectory([],x_r,y_r,t_r);
end

end

% ------------------------------ END OF CODE ------------------------------
