function [time, valuesALL, labels] = simulationLoop(model, iterations, maxTime, maxResets, epsilon, filter, seed_option, plot_option)
%SIMULATIONLOOP Iteratively simulates single trajectories and plots them all
%
% example usage: simulationLoop('mymodel', 10, 20, 5, 0.1, {}, 1, 1)
%
% seed_option:
% 1 = always shuffle the random seed in each iteration
% 2 = make results reproducible
%
% plot_option:
% 0 = one plot
% 1 = subplots for each variable
%
% ------------------------------------------------------------------------------
% authors: Christian Schilling, Taylor T. Johnson
% ------------------------------------------------------------------------------
    %global time values labels ScopeData;

    % add classes and functions subfolders to path
    addpath(genpath('functions'));
    addpath(genpath('classes'));
    
    % open model
%     open_system(['./output/' model]);
    load_system(['./output/' model])
    
    % plot_option enum cases
    plot_phase = 2; % show phase portrait
    
    opt_plot = 1; % 0 = do not plot, 1 = do plot
    
    opt_config_change = 0;
    opt_time = 1;
    
    if opt_time
        tic;
    end
    
    % holds names of model parameters
    printer = com.verivital.hyst.printers.StateflowSpPrinter();
    
    % set standard values
    if ((nargin < 8) || (plot_option == -1))
        plot_option = 1;
    end
    if ((nargin < 7) || (seed_option == -1))
        seed_option = 1;
    end
    if (nargin < 6)
        filter = cell(1, 0);
    end
    if ((nargin < 5) || (epsilon == -1))
        epsilon = 0.0001;
    end
    if ((nargin < 4) || (maxResets == -1))
        maxResets = 3;
    end
    if ((nargin < 3) || (maxTime == -1))
        % try to load maximum simulation time from a model constant
        rt = sfroot;
        modelInstance = rt.find('-isa', 'Simulink.BlockDiagram');
        chart = modelInstance.find('-isa', 'Stateflow.Chart', ...
            '-and', 'Name', ['SF_', model]);
        for i = 1 : 4
            switch i
                % TODO: we should handle this max time detection in Hyst,
                % as this shows up as a common problem in translation
                % between tool formats
                %
                % Solution idea: search like as done here
                % Set a hyst parameter somewhere that corresponds with the
                % maximum time (e.g., as a field)
                % Perhaps this should be done in the configuration options.
                
                case 1
                    x = chart.find('Name', [char(printer.PREFIX_VARIABLE), 'Tmax']);
                case 2
                    x = chart.find('Name', [char(printer.PREFIX_VARIABLE), 'maxT']);
                case 3
                    x = chart.find('Name', [char(printer.PREFIX_VARIABLE), 'tmax']);
                case 4
                    x = chart.find('Name', [char(printer.PREFIX_VARIABLE), 'maxt']);
            end
            if (~ isempty(x))
                break;
            end
        end
        if (~ isempty(x))
            maxTime = str2double(x.Props.InitialValue);
        else
            fprintf('Error: not able to infer a maximum simulation time');
            return;
        end
    end
    
    % deactivate zero crossing detection
%     set_param(model, 'MaxConsecutiveZCsMsg', 'none');
    
    % initialize plot related data
    %figure();
    hold on; % activate several plots in one figure
    
    % set colors
    %colors = ['rgbcmyk'];
    %colors = ['gbcmk'];
    colors = 'k';
    
    % set random seed to default for reproducing results
    if (seed_option)
        rng('default');
    end
    
    % set model parameters
    assignin('base', char(printer.V_maxT.name), maxTime); % maximum simulation time
    assignin('base', char(printer.V_eps.name), epsilon); % epsilon for equality-interval
    assignin('base', char(printer.V_maxResets.name), max(maxResets, 1)); % # backtrackings
    
    % additional configuration
    %simulation_config(model,maxTime,1);
    
    for i = 1 : iterations
        if opt_time
            fprintf(['start iteration: ', num2str(i), ' at ', num2str(toc)]);
        end
        
        % set model parameter: fresh random seed
        assignin('base', char(printer.V_seed.name), randi(1e6));
        
         if opt_config_change
%             configSet_mdl = configSet.copy;
%             configSetNames = get_param(configSet_mdl, 'ObjectParameters');
% 
%             cfg_name = get_param(configSet_mdl, 'Name')
%             try
%                 cfg_max_timestep = get_param(configSet_mdl, 'MaxStep')
%                 set_param(configSet_mdl, 'MaxStep', num2str(rand(1,1) * maxTime / 1000));
%             end
%             set_param(configSet_mdl, 'Name', ['rand', cfg_name, num2str(i)]);
%             cfg_name = get_param(configSet_mdl, 'Name')
%             try
%                 cfg_max_timestep = get_param(configSet_mdl, 'MaxStep');
%             end
%             attachConfigSet(model, configSet_mdl, false);
%             setActiveConfigSet(model, cfg_name);
            
             simulation_config(model,maxTime,rand(1,1));
        end
        
        
        % simulate model
        [~] = sim(model);
        % problem with scopedata if we change parameters here; maybe set
        % active config set elsewhere before simulation loop? doesn't seem
        % to be a problem to set this in general (e.g., manually in the
        % model), just something strange with this particular command
        %[~] = sim(model, 'MaxStep', '1e-4');
        %[~] = sim(model, configSet);
        
        % read data from scope block
        [time, values, labels] = resultsToMatrix(ScopeData, filter, 1);
        
        valuesALL(i).values = values;
        
        plotSimulationLoop;
    end
    
    % deactivate several plots again
    hold off;
end
