scm_example_script.m

function scm_example_script(step)
dbstop if error
%scm_example_script(step)
%
%% Step description:
% Attention: In some of the steps precomputed scm_offline_data can be loaded!
%
% step 1: generates scm_offline_data (stored in reduced_data) and saves them. Shows the typical approach.
% step 2: loads scm_offline_data generated with step 1 and does some tests (scm_lower_bound =< real constant =< scm_upper_bound for random mus and
%         scm_lower_bound = real constant = scm_upper_bound for mu in C)
% step 3: test for scm_offline_data computet with step 1. Compares error estimators computed in the rb_simulation using scm and not using it.
% step 4: test for the small thermal_block (B1 = 2, B2 = 1). One shall see, that the scm for alpha and for beta produce the same upper and lower bounds.
%
%% Notes:
%
% - the fields that can be changed are always marked with 'choose', 'specify' or 'insert' 
%
% - regarding step 1: See 'Additional information on the functionality of the SCM' on which models are supported yet.
%                     The standard model used in this step (thermalblock_model_old) doesn't call gen_reduced_data
%                     in the basis generation (while gen_detailed_data). Other models might do so, which would result in the yet
%                     unsolved problem of generating scm_offline_data in each basis generation step.
% - regarding step 2: Due to the limited amounts of supportet models of step 1 the exact constants in this step are computed rather direct and
%                     expensive (huge 'lin_evol' problems may cause RAM issues here)
% - regarding step 3: If this step shall work for a 'lin_evol' model in the current version (24.09 2012) one has to choos model.error_norm  = 'energy' and
%                     model.rb_simulation = @lin_evol_rb_simulation_primal_dual.
% - regarding step 4: In the end is a display showing the maximum error of any of the computet lower bounds to the respective exact constant. By
%                     this one can see how exact the scm is.
%                     This step can easily be expanded to support other models.
%
%% Additional information on the functionality of the SCM:
% - the SCM for coercive 'lin_stat' problems works fine
% - the SCM for coercive 'lin_evol' problems hasn't been tested yet
% - the SCM for small inf-sup 'lin_stat' problems works fine (see scm_demo.m and step 4)
% - the SCM for small inf-sup 'lin_evol' problems hasn't been tested yet
%
% Huge (in the sense of Q and therefore Q^2 being quite large) inf-sup
% problems (3x3 thermalblock and european_option_pricing_model) caused the
% following problem:
% the bounding box allowed many components to be negativ, which should be
% fine, but in the upcoming online-phase it seemed impossible to guarantee 
% a positiv lower bound. This problem could not be solved in time and will
% have to be tackled in the future. One posibility might be to severely
% increase the size of C (500 and more) via model.scm_size_C.

%% todo: die standard scm values müssen doch ins default model, geht nicht in scm offline -> da dann checks, falls n feld fehlt warning + keyboard oder so

switch step
    
    case 1
        
        model = thermalblock_model_old; % insert the desired 'lin_stat' (or 'lin_evol') model here
        model.verbose = 21; % choose verbosity level (there will be no output in case of model.verbose =< 20)
        % specifying SCM parameters (optional - see scm_offline.m for a detailed description)
        model.scm_M_alpha = 70;
        model.scm_M_plus = 300;
        model.scm_eps_tol = 1e-5;
        model.scm_size_C = 200;
        % necessary to choose either coercive or inf-sup stable SCM
        model.scm_desired_constant = 1; %1 = alpha (coercive), 2 = beta (inf-sup stable)
        % don't change the rest
        model.use_scm = 1;     
        model_data = gen_model_data(model);
        detailed_data = gen_detailed_data(model, model_data);
        reduced_data = gen_reduced_data(model, detailed_data); % scm_offline.m is included here
        save('step1_data.mat')
        
    case 2
        
        load('step1_data.mat') % computet with step 1. Contains model, model_data, detailed_data and reduced_data(including reduced_data.scm_offline_data)
        number_of_test_mus = 50; % choose number of random mus for the test (but 50 should be fine) - nothing else to choose in this step
        K = model.get_inner_product_matrix(model_data);
        if isfield(model_data,'df_info') % in FEM-cases cut out dirgid-values, since they can bias the upcoming eigenvalueproblems
            K(:,model_data.df_info.dirichlet_gids) = [];
            K(model_data.df_info.dirichlet_gids,:) = [];
        end
        % generate a set M_test with the desired number of random mus
        build1 = zeros(size(model.mu_ranges,2),1);
        build2 = zeros(size(model.mu_ranges,2),1);
        for i = 1:size(model.mu_ranges,2)
            build1(i,1) = model.mu_ranges{i}(1);
            build2(i,1) = model.mu_ranges{i}(2);
        end  
        M_test = unifrnd(repmat(build1,1,number_of_test_mus),repmat(build2,1,number_of_test_mus));
        %% first test: scm_lower_bound =< real constant =< scm_upper_bound for mu in M_test
        disp('first test')
        gather_LB_random_mu = zeros(number_of_test_mus,1);
        gather_UB_random_mu = zeros(number_of_test_mus,1);
        gather_constant_random_mu = zeros(number_of_test_mus,1);
        % compute the components (for the computation of the exact constant)
        model.decomp_mode = 1;
        if strcmp(model.rb_problem_type, 'lin_stat') == 1
            A_comp = model.operators(model, model_data);
            M = 1; % only 1 timestep;
            H_total = size(K,1)*M; % total size of the system
        elseif strcmp(model.rb_problem_type, 'lin_evol') == 1
            [L_I_comp,L_E_comp,~] = model.operators_ptr(model, model_data);
            M = model.nt + 1; % number of timesteps
            H_total = size(K,1)*M;  % total size of the system
        end
        warningstate_error = warning('error','MATLAB:eigs:NoEigsConverged'); % sets the state of this warning to error. This way you can tackle it with a try catch block and adjust the problem.
        for i = 1:number_of_test_mus
            fprintf('.');
            model = set_mu(model, M_test(:,i));
            % compute the lower and upper bound
            [gather_LB_random_mu(i),gather_UB_random_mu(i)] = scm_lower_bound(model, reduced_data);
            % compute the exact constant
            model.decomp_mode = 2;
            if strcmp(model.rb_problem_type, 'lin_stat') == 1
                A_coeff = model.operators(model, model_data);
                A = lincomb_sequence(A_comp, A_coeff);
                if isfield(model_data,'df_info') % again cut out dirgid values in FEM case
                    A(:,model_data.df_info.dirichlet_gids) = [];
                    A(model_data.df_info.dirichlet_gids,:) = [];
                end
            elseif strcmp(model.rb_problem_type, 'lin_evol') == 1
                [L_I_coeff,L_E_coeff,~] = model.operators_ptr(model, model_data);
                L_I = lincomb_sequence(L_I_comp, L_I_coeff);
                L_E = lincomb_sequence(L_E_comp, L_E_coeff);
                % build Petrov-Galerkin-BLF-matrix
                kron_1 = sparse(diag([1;zeros(model.nt,1)]));
                kron_2 = sparse(diag([0;ones(model.nt,1)]));
                kron_3 = sparse(diag(-ones(model.nt,1),-1));
                A = kron(kron_1, K) + kron(kron_2, L_I'*K) + kron(kron_3, L_E'*K);
            end
            if model.scm_desired_constant == 1
                try % solve the EWP in a try catch block that prevents covergence issues and other problems with eigs    
                  LM_alpha = eigs(0.5*(A+A'),kron(speye(M),K),1,'LM'); %first try to solve the EWP with default options
                catch err 
                  if (strcmp(err.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err.identifier,'MATLAB:eigs:NoEigsConverged'))...
                    || (strcmp(err.message,'Error with ARPACK routine dnaupd: info = -8'))
                    % in these cases (convergence problems or heavy clustering of eigenvalues) change the options and try again. Can be timeconsuming.
                    if H_total < 600 % opts.p has to be smaller than H_total
                        opts.p = ceil(H_total/2);
                    else
                        opts.p = 300;
                    end
                    opts.maxit = 1500;
                    opts.tol = 1e-14;
                    no_solution = 1;
                    while no_solution
                        disp('There was no solution to an eigenvalueproblem. I changed the options.');
                        try
                            LM_alpha = eigs(0.5*(A+A'),kron(speye(M),K),1,'LM',opts);
                            no_solution  = 0; %end when there is no error
                        catch err2 %when it still fails change the options and try again!
                            no_solution = 1;
                            if (strcmp(err2.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err2.identifier,'MATLAB:eigs:NoEigsConverged'))...
                                || (strcmp(err2.message,'Error with ARPACK routine dnaupd: info = -8'))
                                opts.tol = 100 * opts.tol;
                                if H_total < 600 % opts.p has to be smaller than H_total
                                    opts.p = opts.p + ceil(H_total/50);
                                else
                                    opts.p = opts.p +50;
                                end
                                opts.maxit = opts.maxit + 100;
                            end
                            if opts.tol > 1e-3 %if the options have changed several times and there is still no solution produce an error.
                                error(['I changed the options of an eigenvalueproblem up to opts.tol = ' mat2str(opts.tol) ', opts.p = ' mat2str(opts.p) ' and opts.maxit = ' mat2str(opts.maxit) ' but I still could not find a solution. Reconsider your problem']);
                            end
                        end
                    end
                  else %if there is another error rethrow
                        rethrow(err)
                  end
                end
                
                if LM_alpha >= 0
                    try % solve the EWP in a try catch block that prevents covergence issues and other problems with eigs
                      gather_constant_random_mu(i) = eigs(0.5*(A+A') - (LM_alpha+1e-10)*kron(speye(M),K),kron(speye(M),K),1,'LM') + (LM_alpha+1e-10); %first try to solve the EWP with default options
                    catch err 
                      if (strcmp(err.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err.identifier,'MATLAB:eigs:NoEigsConverged'))...
                        || (strcmp(err.message,'Error with ARPACK routine dnaupd: info = -8'))
                        % in these cases (convergence problems or heavy clustering of eigenvalues) change the options and try again. Can be timeconsuming.
                        if H_total < 600 % opts.p has to be smaller than H_total
                            opts.p = ceil(H_total/2);
                        else
                            opts.p = 300;
                        end
                        opts.maxit = 1500;
                        opts.tol = 1e-14;
                        no_solution = 1;
                        while no_solution
                            disp('There was no solution to an eigenvalueproblem. I changed the options.');
                            try
                                gather_constant_random_mu(i) = eigs(0.5*(A+A') - (LM_alpha+1e-10)*kron(speye(M),K),kron(speye(M),K),1,'LM',opts) + (LM_alpha+1e-10);
                                no_solution  = 0; %end when there is no error
                            catch err2 %when it still fails change the options and try again!
                                no_solution = 1;
                                if (strcmp(err2.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err2.identifier,'MATLAB:eigs:NoEigsConverged'))...
                                    || (strcmp(err2.message,'Error with ARPACK routine dnaupd: info = -8'))
                                    opts.tol = 100 * opts.tol;
                                    if H_total < 600 % opts.p has to be smaller than H_total
                                        opts.p = opts.p + ceil(H_total/50);
                                    else
                                        opts.p = opts.p +50;
                                    end
                                    opts.maxit = opts.maxit + 100;
                                end
                                if opts.tol > 1e-3 %if the options have changed several times and there is still no solution produce an error.
                                    error(['I changed the options of an eigenvalueproblem up to opts.tol = ' mat2str(opts.tol) ', opts.p = ' mat2str(opts.p) ' and opts.maxit = ' mat2str(opts.maxit) ' but I still could not find a solution. Reconsider your problem']);
                                end
                            end
                        end
                      else %if there is another error rethrow
                        rethrow(err)
                      end
                    end
                else
                    gather_constant_random_mu(i) = LM_alpha;
                end
            elseif model.scm_desired_constant == 2
                try % solve the EWP in a try catch block that prevents covergence issues and other problems with eigs
                  LM_beta = sqrt(eigs(A*kron(speye(M),inv(K))*A',kron(speye(M),K),1,'LM')); %first try to solve the EWP with default options
                catch err 
                  if (strcmp(err.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err.identifier,'MATLAB:eigs:NoEigsConverged'))...
                    || (strcmp(err.message,'Error with ARPACK routine dnaupd: info = -8'))
                    % in these cases (convergence problems or heavy clustering of eigenvalues) change the options and try again. Can be timeconsuming.
                    if H_total < 600 % opts.p has to be smaller than H_total
                        opts.p = ceil(H_total/2);
                    else
                        opts.p = 300;
                    end
                    opts.maxit = 1500;
                    opts.tol = 1e-14;
                    no_solution = 1;
                    while no_solution
                        disp('There was no solution to an eigenvalueproblem. I changed the options.');
                        try
                            LM_beta = sqrt(eigs(A*kron(speye(M),inv(K))*A',kron(speye(M),K),1,'LM',opts));
                            no_solution  = 0; %end when there is no error
                        catch err2 %when it still fails change the options and try again!
                            no_solution = 1;
                            if (strcmp(err2.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err2.identifier,'MATLAB:eigs:NoEigsConverged'))...
                                || (strcmp(err2.message,'Error with ARPACK routine dnaupd: info = -8'))
                                opts.tol = 100 * opts.tol;
                                if H_total < 600 % opts.p has to be smaller than H_total
                                    opts.p = opts.p + ceil(H_total/50);
                                else
                                    opts.p = opts.p +50;
                                end
                                opts.maxit = opts.maxit + 100;
                            end
                            if opts.tol > 1e-3 %if the options have changed several times and there is still no solution produce an error.
                                error(['I changed the options of an eigenvalueproblem up to opts.tol = ' mat2str(opts.tol) ', opts.p = ' mat2str(opts.p) ' and opts.maxit = ' mat2str(opts.maxit) ' but I still could not find a solution. Reconsider your problem']);
                            end
                        end
                    end
                  else %if there is another error rethrow
                        rethrow(err)
                  end
                end
                if LM_beta >= 0
                    try % solve the EWP in a try catch block that prevents covergence issues and other problems with eigs 
                      gather_constant_random_mu(i) = sqrt(eigs(A*kron(speye(M),inv(K))*A' - (LM_beta+1e-10)*kron(speye(M),K),kron(speye(M),K),1,'LM') + (LM_beta+1e-10)); %first try to solve the EWP with default options
                    catch err 
                      if (strcmp(err.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err.identifier,'MATLAB:eigs:NoEigsConverged'))...
                        || (strcmp(err.message,'Error with ARPACK routine dnaupd: info = -8'))
                        % in these cases (convergence problems or heavy clustering of eigenvalues) change the options and try again. Can be timeconsuming.
                        if H_total < 600 % opts.p has to be smaller than H_total
                            opts.p = ceil(H_total/2);
                        else
                            opts.p = 300;
                        end
                        opts.maxit = 1500;
                        opts.tol = 1e-14;
                        no_solution = 1;
                        while no_solution
                            disp('There was no solution to an eigenvalueproblem. I changed the options.');
                            try
                                gather_constant_random_mu(i) = sqrt(eigs(A*kron(speye(M),inv(K))*A' - (LM_beta+1e-10)*kron(speye(M),K),kron(speye(M),K),1,'LM',opts) + (LM_beta+1e-10));
                                no_solution  = 0; %end when there is no error
                            catch err2 %when it still fails change the options and try again!
                                no_solution = 1;
                                if (strcmp(err2.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err2.identifier,'MATLAB:eigs:NoEigsConverged'))...
                                    || (strcmp(err2.message,'Error with ARPACK routine dnaupd: info = -8'))
                                    opts.tol = 100 * opts.tol;
                                    if H_total < 600 % opts.p has to be smaller than H_total
                                        opts.p = opts.p + ceil(H_total/50);
                                    else
                                        opts.p = opts.p +50;
                                    end
                                    opts.maxit = opts.maxit + 100;
                                end
                                if opts.tol > 1e-3 %if the options have changed several times and there is still no solution produce an error.
                                    error(['I changed the options of an eigenvalueproblem up to opts.tol = ' mat2str(opts.tol) ', opts.p = ' mat2str(opts.p) ' and opts.maxit = ' mat2str(opts.maxit) ' but I still could not find a solution. Reconsider your problem']);
                                end
                            end
                        end
                      else %if there is another error rethrow
                        rethrow(err)
                      end
                    end
                else
                    gather_constant_random_mu(i) = LM_beta;
                end
            end
        end
        warning(warningstate_error) % return the state of the warning to 'on' instead of 'error'
        fprintf('\n');
        
        %% second test: scm_lower_bound = real constant = scm_upper_bound for mu in C
        disp('second test')
        size_of_C = size(reduced_data.scm_offline_data.C,2);
        gather_LB_mu_in_C = zeros(size_of_C,1);
        gather_UB_mu_in_C = zeros(size_of_C,1);
        gather_constant_mu_in_C = reduced_data.scm_offline_data.alpha_C; % they were already computed in the offline-phase
        % compute the lower and upper bounds
        for j = 1:size_of_C
            fprintf('.');
            model = set_mu(model, reduced_data.scm_offline_data.C(:,j));
            [gather_LB_mu_in_C(j),gather_UB_mu_in_C(j)] = scm_lower_bound(model, reduced_data);
        end
        fprintf('\n');
        
        %% plots and output
        disp('plotting...')
        f = figure;
        subplot(1,2,1)
        plot(1:number_of_test_mus, gather_LB_random_mu, 'r:', 1:number_of_test_mus, gather_constant_random_mu, 'b', 1:number_of_test_mus, gather_UB_random_mu, 'g--')
        title('results of the first test, i.e. using the set of random \mu')
        xlabel('number of test \mu')
        axis tight
        leg = legend('scm lower bound','exact constant','scm upper bound');
        set(leg,'Location','Best')
        subplot(1,2,2)
        plot(1:size_of_C, gather_LB_mu_in_C, 'r:', 1:size_of_C, gather_constant_mu_in_C, 'b', 1:size_of_C, gather_UB_mu_in_C, 'g--')
        title('results of the second tests, i.e. using the set C')
        xlabel('number of \mu in C')
        axis tight
        leg = legend('scm lower bound','exact constant','scm upper bound');
        set(leg,'Location','Best')
        set(gcf,'Color','white')
        % additional output information
        disp(['maximum distance between scm lower bound and exact constant for random mu is ' mat2str(max(abs(gather_constant_random_mu - gather_LB_random_mu)))])
        disp(['maximum distance between scm upper bound and exact constant for random mu is ' mat2str(max(abs(gather_constant_random_mu - gather_UB_random_mu)))])
        disp(['maximum distance between scm lower bound and exact constant for mu in C is ' mat2str(max(abs(gather_constant_mu_in_C - gather_LB_mu_in_C)))])
        disp(['maximum distance between scm upper bound and exact constant for mu in C is ' mat2str(max(abs(gather_constant_mu_in_C - gather_UB_mu_in_C)))])
        
    case 3
        
        load('step1_data.mat') % load offline_data computet with step 1 - includes model, model_data, detailed_data and reduced_data
        number_of_tests = 100; % choose the size of the random mu set
        gather_scm_estimators = zeros(number_of_tests,1);
        gather_non_scm_estimators = zeros(number_of_tests,1);
        % build the random mu set
        build1 = zeros(size(model.mu_ranges,2),1);
        build2 = zeros(size(model.mu_ranges,2),1);
        for i = 1:size(model.mu_ranges,2)
            build1(i,1) = model.mu_ranges{i}(1);
            build2(i,1) = model.mu_ranges{i}(2);
        end  
        M_test = unifrnd(repmat(build1,1,number_of_tests),repmat(build2,1,number_of_tests));
        for i = 1:number_of_tests
            fprintf('.');
            model = set_mu(model, M_test(:,i));
            % first do rb_simulation (error estimator computation) without scm
            model.use_scm = 0;
            rb_sim_data = rb_simulation(model, reduced_data);
            gather_non_scm_estimators(i) = rb_sim_data.Delta;
            % now with scm
            model.use_scm = 1;
            rb_sim_data = rb_simulation(model, reduced_data);
            gather_scm_estimators(i) = rb_sim_data.Delta;
        end
        fprintf('\n');
        % compare results with plot
        disp('plotting...')
        f = figure;
        plot(1:number_of_tests, gather_non_scm_estimators, 'r:', 1:number_of_tests, gather_scm_estimators, 'g--')
        xlabel('number of test \mu')
        leg = legend('estimator without scm','estimator with scm');
        set(leg,'Location','Best')
        set(gcf,'Color','white')
        % additional output information
        disp(['maximum error bestween an estimator with and without scm is ' mat2str(max(abs(gather_non_scm_estimators - gather_scm_estimators)))])

    case 4
        
        mode = 2; % choose to either load precomputed scm_offline_data of the small thermal_block (params.B1 = 2, params.B2 = 1) for the coercive and inf-sup stable case (mode = 1),
                  % or produce the data here (mode = 2)
        if mode == 1
          load('step1_data_for_small_coercive_TB.mat')
          model_1 = model;
          %reduced_data1 = reduced_data;
          reduced_data1 = [];
          reduced_data1.scm_offline_data = scm_offline_data;
          load('step1_data_for_small_infup_stable_TB.mat')
          model_2 = model;
          %reduced_data2 = reduced_data;
          reduced_data2 = [];
          reduced_data2.scm_offline_data = scm_offline_data;
        elseif mode == 2
          params.B1 = 2;
          params.B2 = 1;
          model_1 = thermalblock_model_old(params);
          model_2 = thermalblock_model_old(params);
          % choose to have output (=< 20) or to have output (>20) from the upcoming offline phase
          model_1.verbose = 21; 
          model_2.verbose = 21;
          % make sure to use the scm
          model_1.use_scm = 1;
          model_2.use_scm = 1;
          % configuration fields of the scm
          model_1.scm_M_alpha = 70;
          model_2.scm_M_alpha = 70;
          model_1.scm_M_plus = 200;
          model_2.scm_M_plus = 200;
          model_1.scm_eps_tol = 1e-14;
          model_2.scm_eps_tol = 1e-14;
          model_1.scm_size_C = 50;
          model_2.scm_size_C = 50;
          model_1.scm_desired_constant = 1; % for the coercive scm
          model_2.scm_desired_constant = 2; % for the inf-sup stable scm
          % for the other optional scm-parameters the standard values are sufficient
          % model and detailed_data are the sam in both cases
          model_data = gen_model_data(model_1);
          model_1.RB_mu_list = {[0.1,1],[1,0.1]}; % required so gen_detailed_data works
          detailed_data = gen_detailed_data(model_1, model_data);
          reduced_data1 = gen_reduced_data(model_1, detailed_data);
          reduced_data2 = gen_reduced_data(model_2, detailed_data);
        end
        
        number_of_tests = 100; % choose the size of the random mu set
        gather_lower_bounds_1 = zeros(size(number_of_tests,2));
        gather_upper_bounds_1 = zeros(size(number_of_tests,2));
        gather_lower_bounds_2 = zeros(size(number_of_tests,2));
        gather_upper_bounds_2 = zeros(size(number_of_tests,2));
        gather_constant = zeros(size(number_of_tests,2));
        % build the random mu set
        build1 = zeros(size(model_1.mu_ranges,2),1);
        build2 = zeros(size(model_1.mu_ranges,2),1);
        for i = 1:size(model_1.mu_ranges,2)
            build1(i,1) = model_1.mu_ranges{i}(1);
            build2(i,1) = model_1.mu_ranges{i}(2);
        end  
        M_test = unifrnd(repmat(build1,1,number_of_tests),repmat(build2,1,number_of_tests));
        % collect the exact constant and the various upper and lower bounds
        % for every mu in the set
        for  i = 1:number_of_tests
            fprintf('.');
            model_1 = set_mu(model_1, M_test(:,i));
            [gather_lower_bounds_1(i), gather_upper_bounds_1(i)] = scm_lower_bound(model_1, reduced_data1);
            model_2 = set_mu(model_2, M_test(:,i));
            [gather_lower_bounds_2(i), gather_upper_bounds_2(i)] = scm_lower_bound(model_2, reduced_data2);
            gather_constant(i) = min(M_test(:,i)); % for the thermal block this is known to always be the exact constant
        end
        fprintf('\n');
        % plot the upper and lower bounds
        disp('plotting...')
        f = figure;
        subplot(1,2,1)
        plot(1:number_of_tests, gather_lower_bounds_1, 'r:', 1:number_of_tests, gather_lower_bounds_2, 'g--')
        xlabel('number of test \mu')
        leg = legend('\alpha_{scm,LB}','\beta_{scm,LB}');
        set(leg,'Location','Best')
        subplot(1,2,2)
        plot(1:number_of_tests, gather_upper_bounds_1, 'r:', 1:number_of_tests, gather_upper_bounds_2, 'g--')
        xlabel('number of test \mu')
        leg = legend('\alpha_{scm,UB}','\beta_{scm,UB}');
        set(leg,'Location','Best')
        set(gcf,'Color','white') 
        % error result (how exact is the lower bound to the exact constant)
        disp(['maximal error of a lower bound to the exact constant is ' mat2str(max(abs(gather_lower_bounds_1 - gather_constant)))])
        
    otherwise
        
        error('step number unknown!')
        
end