convection_diffusion_fv_output_opt_script.m

function convection_diffusion_fv_output_opt_script(step, expar)
%function convection_diffusion_fv_output_opt_script(step)
%
%

% Markus Dihlmann 12.09.2011


if nargin<2
    expar = [];
end


fdir = fullfile(rbmatlabhome,'rbasis','lin_evol_opt','convdiff');
fdir_temp = fullfile(rbmatlabtemp,'lin_evol_opt_data');

switch step
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5    
% Basis generation    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 1
        
        if nargin >1
            if isfield(expar,'coarse_factor')
                params.coarse_factor = expar.coarse_factor;
            else
                params.coarse_factor = 16;
            end
        else
            params.coarse_factor = 16;
        end
        
        model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
        model.compute_derivative_info=0;
        
        if nargin>1
            if isfield(expar,'k')
                model.k = expar.k;
            else
                model.k = 0.25;
            end
        end
        
        
        model_data = gen_model_data(model);
        
        model.RB_stop_Nmax = 50;
        
        
        detailed_data = gen_detailed_data(model, model_data);
        
        
        save(fullfile(fdir,['DiffGreedybase_coarse',num2str(model.coarse_factor),'_Nmax_',num2str(model.RB_stop_Nmax),'_k_',num2str(model.k),'.mat']),...
            'model', 'detailed_data');
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate separate reduced basis with constant diffusion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 2
            
        
        if nargin >1
            if isfield(expar,'coarse_factor')
                params.coarse_factor = expar.coarse_factor;
            else
                params.coarse_factor = 16;
            end
        else
            params.coarse_factor = 16;
        end
        params.separate_bases = 1;
        params.cone_const = 1;
        
        model = convection_diffusion_fv_output_opt_model_rect_normed(params);
        model.compute_derivative_info=1;
        model.RB_error_indicator = 'estimator_parallel';
        
        
        if nargin>1
            if isfield(expar,'k')
                model.k = expar.k;
            else
                model.k = 0.25;
            end
            if isfield(expar,'gen_detailed_data')
                model.gen_detailed_data = expar.gen_detailed_data;
                model.RB_extension_algorithm = @RB_extension_PCA_fixspace_without_filecache;
            end
        else
            model.k = 0.25;
        end
        
        
        model_data = gen_model_data(model);
        
        model.RB_stop_Nmax = 10;
        
        
        detailed_data = gen_detailed_data(model, model_data);
        
        
        save(fullfile(fdir,['DiffGreedybase_sep_coarse',num2str(model.coarse_factor),'_Nmax_',num2str(model.RB_stop_Nmax),'_k_',num2str(model.k),'.mat']),...
            'model', 'detailed_data');
    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% detailed gradient optimization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 3
        params.coarse_factor = 4;
        params.cone_const = 1;
        
        model = convection_diffusion_fv_output_opt_model_rect_normed(params);
        model_data = gen_model_data(model);
        disp(['Diffusion coefficient = ',num2str(model.k)]);
        
        model.optimization.init_params = [0.3,0.3];
        model.compute_derivative_info=1;
        tic
        opt_data = optimize(model, model_data);
        t_opt_d = toc;
        opt_data.t_opt_d = t_opt_d;
        
        disp(['Time for optimization : ',num2str(t_opt_d)]);
        
        keyboard;
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% grid optimization detailed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 4
        disp('grid optimization of diffusion model')
        
        params.coarse_factor = 4;
        params.cone_const = 1;
        
        model = convection_diffusion_fv_output_opt_model_rect_normed(params);
        model.k=0.25;
        model_data = gen_model_data(model);
        disp(['Diffusion coefficient = ',num2str(model.k)]);
        model.mu_ranges = {[0.788,0.86],[0.6720,0.704]};
        
        model.optimization.optimizer = @detailed_grid_search;
        model.opt_params.grid_density = 6;
        
        
        opt_data = optimize(model, model_data);
        
        
        keyboard;
        
        
        

        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% reduced optimization (simplified) using a specialized reduced basis
% and a cone_const conv_diff_model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 5
        
        load('DiffGreedyBase_ptpart10_2p_coarse4_Nmax30_k_0.25_selected.mat');
        
%         params.coarse_factor = 4;
%         params.separate_bases = 1;
%         params.cone_const = 1;
%         
%         model = convection_diffusion_fv_output_opt_model_rect_normed(params);
%         
        model.optimization.init_params = [0.3,0.3];
        model.optimization.opt_mode = 'reduced';
        model.compute_derivative_info = 1;
        model.optimization.lower_bound = [0,0];
        model.optimization.upper_bound = [1,1];
        model.lipschitz_type = 'norm';
        
        reduced_data = gen_reduced_data(model, detailed_data);
         
        tic;
        opt_data = optimize(model,reduced_data);
        t_opt_red = toc;
        opt_data.t_opt_red = t_opt_red;
        disp(['Time for reduced optimization: ',num2str(t_opt_red)]);
        plot_PM(opt_data);
        
        keyboard;
    
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% reduced optimization coarse 4 using different sizes of bases
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 501
        
        load('DiffGreedyBase_ptpart10_2p_coarse4_Nmax30_k_0.25_selected.mat');
        load('opt_data_graddet_0303_q10_tol0.01.mat');
        opt_data_det = opt_data;
                
        part_used = [34,45,67,68,68,69,69,69];
        nr_tparts = length(detailed_data.parts_detailed_data{part_used(1)}.t_part_detailed_data);
%         params.coarse_factor = 4;
%         params.separate_bases = 1;
%         params.cone_const = 1;
%         
%         model = convection_diffusion_fv_output_opt_model_rect_normed(params);
%         
        model.optimization.init_params = [0.3,0.3];
        model.optimization.opt_mode = 'reduced';
        model.compute_derivative_info = 1;
        model.optimization.lower_bound = [0,0];
        model.optimization.upper_bound = [1,1];
        model.lipschitz_type = 'norm';
        
        %add reduced data subset fields
        model.reduced_data_subset = @p_part_reduced_data_subset;
        model.base_model.base_model.reduced_data_subset = @lin_evol_opt_reduced_data_subset_separate_bases;
        
        
        reduced_data = gen_reduced_data(model, detailed_data);
         
        
        N_vec = 55:1:60;
        N_der_vec = 25:1:30;
        
        %quantities to note
        n_N = length(N_vec);
        n_Nder = length(N_der_vec);
        
        N_av = zeros(1,n_N*n_Nder);
        N_der_av = zeros(2,n_N*n_Nder);
        t_opt = zeros(1,n_N*n_Nder);
        opt_err = zeros(1,n_N*n_Nder);
        opt_delta = zeros(1,n_N*n_Nder);
        nbr_steps = zeros(1,n_N*n_Nder);
        
        for i=1:length(N_vec)
            for j=1:length(N_der_vec);
                model.N = N_vec(i);
                model.N_der = N_der_vec(j).*[1,1];
                
                red_data_sub = model.reduced_data_subset(model, reduced_data);
                
                %optimization
                tic;
                opt_data = optimize(model,red_data_sub);
                t_opt_red = toc;
                opt_data.t_opt_red = t_opt_red;
                disp(['Time for reduced optimization: ',num2str(t_opt_red)]);
                
                %save quantities
                t_opt(n_Nder*(i-1)+j)=t_opt_red;
                opt_err(n_Nder*(i-1)+j) = norm(opt_data.max_output_paramsets - opt_data_det.max_output_paramsets);
                opt_delta(n_Nder*(i-1)+j) = opt_data.Delta_mu(end);
                nbr_steps(n_Nder*(i-1)+j) = opt_data.nr_opti_iter;
                
                sumN=[0,0,0];
                for p=1:length(part_used)
                    rd=red_data_sub.parts_reduced_data{part_used(p)};
                    for k=1:length(rd.t_part_reduced_data)
                        sumN(1) = sumN(1)+length(rd.t_part_reduced_data{k}.a0{1});
                        sumN(2) = sumN(2)+length(rd.t_part_reduced_data{k}.c0{1}{1});
                        sumN(3) = sumN(3)+length(rd.t_part_reduced_data{k}.c0{2}{1});
                    end
                end
        
                N=size(part_used,2)*nr_tparts;
                sumN = sumN./N;
                N_av(n_Nder*(i-1)+j) = sumN(1);
                N_der_av(:,n_Nder*(i-1)+j) = sumN([2,3])';
                
            end
        end
        
        disp('Results:')
        disp('---------------------');
        disp('N_av , N_av_der, t_opt, opt_error, opt est err, numbr steps')
        for i=1:length(N_av)
            disp([' ',num2str(N_av(i)),' | ',num2str(N_der_av(1,i)),' | ',num2str(N_der_av(2,i)),...
               ' || ',num2str(t_opt(i)),' | ',num2str(opt_err(i)),' | ',num2str(opt_delta(i))...
               ' | ',num2str(nbr_steps(i)),'||']);
        end
        
        
        keyboard;
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% detailed gradient optimization (high dimensionsl coarse=1 oder 2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 6
        params.coarse_factor = 2;
        params.cone_const = 1;
        
        model = convection_diffusion_fv_output_opt_model_rect_normed(params);
        model_data = gen_model_data(model);
        disp(['Diffusion coefficient = ',num2str(model.k)]);
        
        model.optimization.init_params = [0.3, 0.3];
        model.compute_derivative_info=1;
        tic
        opt_data = optimize(model, model_data);
        t_opt_d = toc;
        opt_data.t_opt_d = t_opt_d;
        
        disp(['Time for optimization : ',num2str(t_opt_d)]);
        opt_data_det = opt_data;
        
        save(['opt_data_det_gradient_step_error_init',num2str(model.optimization.init_params),'.mat'],'opt_data_det');
        
        keyboard;
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% reduced optimization (gradient step error estimator)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 7
        
        load(fullfile(fdir,'coarse2','reduced_data_complete_N200.mat'));
        
        
%         params.coarse_factor = 4;
%         params.separate_bases = 1;
%         params.cone_const = 1;
%         
%         model = convection_diffusion_fv_output_opt_model_rect_normed(params);
%       

        %load(fullfile(fdir,'coarse2','basisNmax200','reduced_data_complete_2ndfunctional.mat'));
        
        
        model.optimization.init_params = [0.3,0.3];
        model.optimization.opt_mode = 'reduced';
        model.base_model.optimization.init_params = [0.3,0.3];
        model.base_model.optimization.opt_mode = 'reduced';
        model.base_model.base_model.optimization.init_params = [0.3,0.3];
        model.base_model.base_model.optimization.opt_mode = 'reduced';
        
         model.optimization.tol = 0.001;
        model.base_model.optimization.tol = 0.001;
        model.base_model.base_model.optimization.tol =0.001;
         
        %___________________________________________________________________
        %functional 2
        
        
        
       
        %set 2nd functional fields
%         model.base_model.save_time_indices = 0:model.base_model.nt;
%         model.base_model.base_model.save_time_indices = 0:model.base_model.nt;
%         model.verbose = 4;
%         model.base_model.verbose = 4;
%         model.name_output_functional = 'box_mean_time';
%         model.get_output = @lin_evol_opt_output_time_integrated;
%         model.output_function_ptr = @output_function_box_mean;
%         model.sbox_xmin = 0.5;
%         model.sbox_xmax = 1;
%         model.sbox_ymin = 0.5;
%         model.sbox_ymax = 1;
%         model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
%         model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
%         model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
%         model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
%         model.base_model.name_output_functional = 'box_mean_time';
%         model.base_model.get_output = @lin_evol_opt_output_time_integrated;
%         model.base_model.output_function_ptr = @output_function_box_mean;
%         model.base_model.sbox_xmin = 0.5;
%         model.base_model.sbox_xmax = 1;
%         model.base_model.sbox_ymin = 0.5;
%         model.base_model.sbox_ymax = 1;
%         model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
%         model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
%         model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
%         model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
%         model.base_model.base_model.name_output_functional = 'box_mean_time';
%         model.base_model.base_model.get_output = @lin_evol_opt_output_time_integrated;
%         model.base_model.base_model.output_function_ptr = @output_function_box_mean;
%         model.base_model.base_model.sbox_xmin = 0.5;
%         model.base_model.base_model.sbox_xmax = 1;
%         model.base_model.base_model.sbox_ymin = 0.5;
%         model.base_model.base_model.sbox_ymax = 1;
%         model.base_model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
%         model.base_model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
%         model.base_model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
%         model.base_model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        %set optimzation parameters
       
        model.base_model.optimization.optimizer = @gradient_opt;
        model.base_model.base_model.optimization.optimizer = @gradient_opt;
        model.optimization.stepsize_rule = 'quotient';        
        model.base_model.optimization.stepsize_rule = 'quotient';
        model.base_model.base_model.optimization.stepsize_rule = 'quotient';
        model.optimization.initial_stepsize = 1;
        model.base_model.optimization.initial_stepsize = 1;
        model.base_model.base_model.optimization.initial_stepsize = 1;
        model.stepsize_coefficient = 3;
        model.base_model.stepsize_coefficient = 3;
        model.base_model.base_model.stepsize_coefficient = 3;
        
%         model.optimization.stepsize_rule = 'armijo';
%         model.optimization.armijo_t_min = 0.8;
%         model.optimization.stepsize_factor = 0.75;
%         model.base_model.optimization.stepsize_rule = 'armijo';
%         model.base_model.optimization.armijo_t_min = 0.8;
%         model.base_model.optimization.stepsize_factor = 0.75;
%         model.base_model.base_model.optimization.stepsize_rule = 'armijo';
%         model.base_model.base_model.optimization.armijo_t_min = 0.8;
%         model.base_model.base_model.optimization.stepsize_factor = 0.75;
%         
        
        tic;
        opt_data = optimize(model,reduced_data);
        t_opt_red = toc;
        opt_data.t_opt_red = t_opt_red;
        disp(['Time for reduced optimization: ',num2str(t_opt_red)]);
        plot_PM(opt_data);
        
        keyboard;
   
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conducting a reduced simulation on symphony
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 8
        
        cd /home/dihlmann/matlab/rbmatlab-core/rbasis/lin_evol_opt/convdiff/coarse2
        
        load('reduced_data_complete.mat');
        
        model.optimization.init_params = [0.3,0.3];
        model.optimization.opt_mode = 'reduced';
        model.compute_derivative_info = 1;
        model.optimization.lower_bound = [0,0];
        model.optimization.upper_bound = [1,1];
        model.lipschitz_type = 'norm';
        
        tic;
        opt_data = optimize(model,reduced_data);
        t_opt_red = toc;
        opt_data.t_opt_red = t_opt_red;
        disp(['Time for reduced optimization: ',num2str(t_opt_red)]);
        save(['opt_data_',num2str(model.optimization.init_params),'.dat'],'opt_data');
        
        disp('opt_data gespeichert')
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% doing several reduced optimization for randomly chosen initial values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 9
        
        %load('reduced_data_complete.mat');
        load(fullfile(fdir_temp,'coarse2','reduced_data_complete.mat'));
        model.optimization.opt_mode = 'reduced';
        model.base_model.optimization.opt_mode = 'reduced';
        model.base_model.base_model.optimization.opt_mode = 'reduced';
        
        n_sample = 30;
        rand_mu=rand_uniform(n_sample,model.mu_ranges);
        stepsizes = 5:10;
        
        opt_data_ar = cell(n_sample, length(stepsizes));
        
        for i=1:n_sample
            model.optimization.init_params = rand_mu(:,i);
            for j=1:length(stepsizes)
                model.stepsize_coefficient = stepsizes(j);
                tic;
                opt_data_ar{i}{j} = optimize(model, reduced_data);
                t_opt = toc;
                opt_data_ar{i}{j}.t_opt = t_opt;
            end
        end
        
        save('opt_data_arrays25.mat','opt_data_ar','rand_mu','stepsizes');

        %find out best step size
        stepsize_v = zeros(n_sample,1);
        Dmu_min_v = zeros(n_sample,1);
        Dmu_sum = 0;
        for i=1:n_sample
            Dmu_min=opt_data_ar{i}{1}.Delta_mu(end);
            stepsize_v(i) = stepsizes(1);
            for j=2:length(stepsizes)
                if (opt_data_ar{i}{j}.Delta_mu(end)<Dmu_min)
                    Dmu_min = opt_data_ar{i}{j}.Delta_mu(end);
                    stepsize_v(i) = stepsizes(j);
                end
            end
            Dmu_min_v(i) = Dmu_min;
            Dmu_sum = Dmu_sum+Dmu_min;
        end
        Dmu_av = Dmu_sum/n_sample;
        
        save('opt_data_arrays_eval.mat','opt_data_ar','rand_mu','stepsizes',...
            'stepsize_v','Dmu_min_v','Dmu_av');
        
        keyboard
        
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% detailed optimization for optimization with new error estimator
% (this is for comparison purposes)
% data is saved in opt_data_det.mat
% !!!exchange model to outflow model as soon as new model is available!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 10
        
        
        params.coarse_factor = 2;
        params.cone_const =1;
        
%         model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
%         model.verbose = 5;
%         model.optimization.optimizer = @gradient_opt_non_iter_err;
%         model.optimization.tol = 0.0005;
%         model.optimization.init_params = [0.3, 0.3];
%         %model.optimization.stepsize_rule = 'quotient';
%         %model.stepsize_coefficient = 7;
%         %model.optimization.stepsize_rule = 'fixed_stepsize';
%         %model.optimization.stepsize_fix = 1;
%         model.optimization.stepsize_rule = 'armijo';
%         model.optimization.sigma = 0.8;
%         model.optimization.stepsize_factor = 0.75;
%         
%         model_data = gen_model_data(model);
%         
%         
%         tic
%         opt_data_det = optimize(model, model_data);
%         t_opt = toc;
%         
%         opt_data_det.t_opt = t_opt;
%         
%         %keyboard;
%         
%         name = ['opt_data_det_',num2str(model.optimization.init_params(1)),'_',num2str(model.optimization.init_params(2)),...
%             'coarse',num2str(model.coarse_factor),'.mat'];
%         
%         save(name,'opt_data_det');
        
%         keyboard;
        
        % zweites fucktional
        disp('Optimizing 2nd functional')
        disp('------------------------')
        params.functional = 'box_mean_time';
        
        model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
        model.verbose = 9;
        
        model_data = gen_model_data(model);
        
        %set optimization paramters
        model.optimization.tol = 0.00001;
        model.optimization.init_params = [0.8, 0.8];
        model.optimization.stepsize_rule = 'armijo';
        model.optimization.sigma = 0.7;
        model.optimization.stepsize_factor = 0.75;
        model.optimization.initial_stepsize = 50;
        
        tic
        opt_data_det = optimize(model, model_data);
        t_opt = toc;
        
        opt_data_det.t_opt = t_opt;
        
        name = ['opt_data_det_func2_',num2str(model.optimization.init_params(1)),'_',num2str(model.optimization.init_params(2)),...
            'coarse',num2str(model.coarse_factor),'.mat'];
        
        save(name,'opt_data_det');
        
        
        %keyboard;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reduced optimization using the new non-iterative error estimator. We
% verify the error estimator using the result from step 10. The constants
% are assumed to be known at the optimum: replace later by a more general
% offline computation... see step 120 and 121.
% An important quantity for the quality of the error estimtor is the
% breaking condition: model.optimization.tol
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
    case 11
        
        %load detailed optimization data
        %load('opt_data_det_0.3_0.3coarse2.mat');
        
        
        %load the reduced data and model (later: outflow model)
        load(fullfile(fdir,'coarse2','reduced_data_complete_N200.mat'));
        %load(fullfile(fdir,'coarse16','reduced_data_complete_coarse16.mat'));
        
        model.verbose = 4;
        %set optimization parameters
        model.optimization.tol = 0.0005;
        model.base_model.optimization.tol = 0.0005;
        model.base_model.base_model.optimization.tol =0.0005;
        model.optimization.init_params = [0.3,0.3];
        model.base_model.optimization.init_params = [0.3,0.3];
        model.base_model.base_model.optimization.init_params = [0.3,0.3];
        
        %set new fields in model
        model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.base_model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        
        %set fixed stepsize
%         model.optimization.stepsize_rule = 'fixed_stepsize';
%         model.base_model.optimization.stepsize_rule = 'fixed_stepsize';
%         model.base_model.base_model.optimization.stepsize_rule = 'fixed_stepsize';
%         model.optimization.stepsize_fix = 1;
%         model.base_model.optimization.stepsize_fix = 1;
%         model.base_model.base_model.optimization.stepsize_fix = 1;
        %model.optimization.stepsize_rule = 'quotient';
        %model.base_model.optimization.stepsize_rule = 'quotient';
        %model.base_model.base_model.optimization.stepsize_rule = 'quotient';
        model.optimization.stepsize_rule = 'armijo';
        model.optimization.sigma = 0.8;
        model.optimization.stepsize_factor = 0.75;
        model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.optimization.sigma = 0.8;
        model.base_model.optimization.stepsize_factor = 0.75;
        
        %set constants(replace later!!)
        model.gamma_H = @my_get_gamma_H_constant;
        model.base_model.gamma_H = @my_get_gamma_H_constant;
        model.base_model.base_model.gamma_H = @my_get_gamma_H_constant;
        model.L_H = @my_get_L_H_constant;
        model.base_model.L_H = @my_get_L_H_constant;
        model.base_model.base_model.L_H = @my_get_L_H_constant;
        
        %optimization
        tic
        opt_data = optimize(model, reduced_data);
        t_opt = toc;
        opt_data.t_opt = t_opt;
        
         save(fullfile(fdir,'optimization_case11.mat'),'opt_data')
        
        %verify and compare
%         disp(['reduced optimum: ',num2str(opt_data.optimal_params)])
%         disp(['Delta_mu       : ', num2str(opt_data.Delta_mu)])
%         diff = opt_data.optimal_params-opt_data_det.optimal_params;
%         disp(['true error     : ',num2str(norm(diff))])
%         disp(['Efficiency     : ',num2str(opt_data.Delta_mu/norm(diff))])
%         disp('-----')
%         disp(['Time reduced optimization: ', num2str(opt_data.t_opt)])
%         disp(['Time detailed optimization: ', num2str(opt_data_det.t_opt)]);
        
        
        
        %keyboard;
       
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Optimizae with a second functional (but NON TIME DEPENDENT!!!)
%
% the functional fields have to be added manually, as they are not included
% in the large model.
% reduced data has to be inhanced by projected output vectors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 12
        
        %load detailed optimization data
        %load('opt_data_det_func2_0.3_0.3coarse2.mat');
        
        
        %load the reduced data and model (later: outflow model)
        load(fullfile(fdir,'coarse2_neumann','reduced_data_complete_N200.mat'));
        
        %set 2nd functional fields
        model.name_output_functional = 'box_mean';
        model.get_output = @lin_evol_opt_output_time_independent;
        model.output_function_ptr = @output_function_box_integral;
        model.sbox_xmin = 0.75;
        model.sbox_xmax = 1.75;
        model.sbox_ymin = 0.25;
        model.sbox_ymax = 1.25;
        model.s_l2norm = sqrt(2);
        model.base_model.name_output_functional = 'box_mean';
        model.base_model.get_output = @lin_evol_opt_output_time_independent;
        model.base_model.output_function_ptr = @output_function_box_integral;
        model.base_model.sbox_xmin = 0.75;
        model.base_model.sbox_xmax = 1.75;
        model.base_model.sbox_ymin = 0.25;
        model.base_model.sbox_ymax = 1.25;
        model.base_model.s_l2norm = sqrt(2);
        model.base_model.base_model.name_output_functional = 'box_mean';
        model.base_model.base_model.get_output = @lin_evol_opt_output_time_independent;
        model.base_model.base_model.output_function_ptr = @output_function_box_integral;
        model.base_model.base_model.sbox_xmin = 0.75;
        model.base_model.base_model.sbox_xmax = 1.75;
        model.base_model.base_model.sbox_ymin = 0.25;
        model.base_model.base_model.sbox_ymax = 1.25;
        model.base_model.base_model.s_l2norm = sqrt(2);
        
        
        
        %set tolerance level
        model.optimization.tol = 0.0001;
        model.base_model.optimization.tol = 0.0001;
        model.base_model.base_model.optimization.tol =0.0001;
        
        %set new fields in model
        model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.base_model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        
        %set fixed stepsize
        model.optimization.stepsize_rule = 'fixed_stepsize';
        model.base_model.optimization.stepsize_rule = 'fixed_stepsize';
        model.base_model.base_model.optimization.stepsize_rule = 'fixed_stepsize';
        model.optimization.stepsize_fix = 10;
        model.base_model.optimization.stepsize_fix = 10;
        model.base_model.base_model.optimization.stepsize_fix = 10;
        
        %set constants(replace later!!)
%         model.gamma_H = @(model)[2,model];
%         model.base_model.gamma_H = @(model)[2,model];
%         model.base_model.base_model.gamma_H = @(model)[2,model];
%         model.L_H = @(model)[0.08,model];
%         model.base_model.L_H = @(model)[0.08,model];
%         model.base_model.base_model.L_H = @(model)[0.08,model];
%         
        %optimization
        tic
        opt_data = optimize(model, reduced_data);
        t_opt = toc;
        opt_data.t_opt = t_opt;
        
        %verify and compare
        disp('Results for 2nd Functional reduced optimization');
        disp('----------------------------------------------');
        disp(['reduced optimum: ',num2str(opt_data.optimal_params)])
        disp(['Delta_mu       : ', num2str(opt_data.Delta_mu)])
        diff = opt_data.optimal_params-opt_data_det.optimal_params;
        disp(['true error     : ',num2str(norm(diff))])
        disp(['Efficiency     : ',num2str(opt_data.Delta_mu/norm(diff))])
        disp('-----')
        disp(['Time reduced optimization: ', num2str(opt_data.t_opt)])
        disp(['Time detailed optimization: ', num2str(opt_data_det.t_opt)]);
        
        
        
        keyboard;
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Optimize with a second functional
%
% important: the functional fields have to be added manually, as they were
% not included while painstakingly creating the reduced model on Symphony.
% The optimization tolerance has to be set very low (1e-5)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 13
        
        disp('Optimizing 2nd functional by gradient method with non-iterative error estimator')
        
        %load detailed optimization data
        %load('opt_data_det_func2_0.8_0.8coarse2.mat');
        
        
        %load the reduced data and model (later: outflow model)
        load(fullfile(fdir,'coarse2','basisNmax200','reduced_data_complete_2ndfunctional.mat'));
        %load(fullfile(fdir,'coarse16','reduced_data_complete_2ndfunctional.mat'));
        
        %set 2nd functional fields
        model.verbose = 4;
        model.base_model.verbose = 4;
        model.name_output_functional = 'box_mean_time';
        model.get_output = @lin_evol_opt_output_time_integrated;
        model.output_function_ptr = @output_function_box_mean;
        model.sbox_xmin = 0.5;
        model.sbox_xmax = 1;
        model.sbox_ymin = 0.5;
        model.sbox_ymax = 1;
        model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.sbox_xmin = 0.5;
        model.base_model.sbox_xmax = 1;
        model.base_model.sbox_ymin = 0.5;
        model.base_model.sbox_ymax = 1;
        model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.base_model.sbox_xmin = 0.5;
        model.base_model.base_model.sbox_xmax = 1;
        model.base_model.base_model.sbox_ymin = 0.5;
        model.base_model.base_model.sbox_ymax = 1;
        model.base_model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.base_model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.base_model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        %set optimzation parameters
        model.optimization.tol = 0.0005;
        model.base_model.optimization.tol = 0.0005;
        model.base_model.base_model.optimization.tol =0.0005;
        
        model.optimization.init_params =[0.75,0.57];
        model.base_model.optimization.init_params =[0.75,0.57];
        model.base_model.base_model.optimization.init_params =[0.75,0.57];
        
        model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.base_model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        
        model.optimization.stepsize_rule = 'armijo';
        model.optimization.sigma = 0.7;
        model.optimization.stepsize_factor = 0.75;
        model.optimization.initial_stepsize = 10;
        model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.optimization.sigma = 0.7;
        model.base_model.optimization.stepsize_factor = 0.75;
        model.base_model.optimization.initial_stepsize = 10;
        model.base_model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.base_model.optimization.sigma = 0.7;
        model.base_model.base_model.optimization.stepsize_factor = 0.75;
        model.base_model.base_model.optimization.initial_stepsize = 10;
        
        model.L_H = @my_get_L_H_constant;
        model.base_model.L_H = @my_get_L_H_constant;
        model.base_model.base_model.L_H = @my_get_L_H_constant;
        
        model.gamma_H = @my_get_gamma_H_constant;
        model.base_model.gamma_H = @my_get_gamma_H_constant;
        model.base_model.base_model.gamma_H = @my_get_gamma_H_constant;
        
        model.optimization.allow_constant_calculation = 1;
        model.base_model.optimization.allow_constant_calculation = 1;
        model.base_model.base_model.optimization.allow_constant_calculation = 1;
        
        %optimization
        disp('starting reduced optimization')
        tic
        opt_data = optimize(model, reduced_data);
        t_opt = toc;
        opt_data.t_opt = t_opt;
        
        save(fullfile(fdir,'optimization_case13.mat'),'opt_data');
        
        %detailed_optimization
        disp('generating model data...')
        model_data = gen_model_data(model);
        disp('starting detailed_optimization...')
        tic
        opt_data_det = optimize(model, model_data);
        t_opt = toc;
        opt_data_det.t_opt = t_opt;
        

        
        %verify and compare
        disp('Results for 2nd Functional reduced optimization');
        disp('----------------------------------------------');
        disp(['reduced optimum: ',num2str(opt_data.optimal_params)])
        disp(['Delta_mu       : ', num2str(opt_data.Delta_mu)])
        diff = opt_data.optimal_params-opt_data_det.optimal_params;
        disp(['true error     : ',num2str(norm(diff))])
        disp(['Efficiency     : ',num2str(opt_data.Delta_mu/norm(diff))])
        disp('-----')
        disp(['Time reduced optimization: ', num2str(opt_data.t_opt)])
        disp(['Time detailed optimization: ', num2str(opt_data_det.t_opt)]);
        
        save(fullfile(fdir,['opt_data_2ndfunctional_initsteps_',...
            num2str(model.optimization.initial_stepsize),...
            'coarse2.mat']),'opt_data','opt_data_det');
        
        %keyboard;
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% optimize with a basis subset (varying basis size)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 14
        
        N_choice = [76,76,70,65,60,50];
        N_der_choice = [80,80;75,75;70,70;65,65;60,60;50,50];
        
        %load the reduced basis
        load(fullfile(fdir,'coarse2','reduced_data_complete_N200.mat'));
        
        
        %set tolerance level
        model.optimization.tol = 0.0005;
        model.base_model.optimization.tol = 0.0005;
        model.base_model.base_model.optimization.tol =0.0005;
        model.optimization.init_params = [0.3,0.3];
        model.base_model.optimization.init_params = [0.3,0.3];
        model.base_model.base_model.optimization.init_params = [0.3,0.3];
        model.optimization.stepsize_rule = 'armijo';
        model.optimization.sigma = 0.8;
        model.optimization.stepsize_factor = 0.75;
        model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.optimization.sigma = 0.8;
        model.base_model.optimization.stepsize_factor = 0.75;
        
        %set new fields in model
        model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.base_model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        
         %set constants(replace later!!)
        model.gamma_H = @my_get_gamma_H_constant;
        model.base_model.gamma_H = @my_get_gamma_H_constant;
        model.base_model.base_model.gamma_H = @my_get_gamma_H_constant;
        model.L_H = @my_get_L_H_constant;
        model.base_model.L_H = @my_get_L_H_constant;
        model.base_model.base_model.L_H = @my_get_L_H_constant;
        
        %fix reduced data subset function
        model.reduced_data_subset = @p_part_reduced_data_subset;
        
        %normal optimzation full rb
        tic;
        opt_data_ar{1} = optimize(model, reduced_data);
        t_opt = toc;
        opt_data_ar{1}.t_opt = t_opt;
        
        for i=1:length(N_choice)
            model.N=N_choice(i);
            model.N_der = N_der_choice(i,:);
            reduced_data_subset = model.reduced_data_subset(model,reduced_data);
            
            tic;
            opt_data_ar{i+1} = optimize(model, reduced_data_subset);
            t_opt = toc;
            opt_data_ar{i+1}.t_opt = t_opt;
            clear reduced_data_subset
            %keyboard
        end
        
        save(fullfile(fdir,'case14_sizecomparison.mat'),'opt_data_ar','N_choice','N_der_choice')
        %keyboard;
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
% performing a detailled and a reduced optimization 
% for various initial prameters (randomly chosen)
% the evaluating the results
% functional 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 15
        
        %number of samples
        if (nargin>1)&&(isfield(expar,'n_sample'))
            n_sample = expar.n_sample;
        else
            n_sample = 5;
        end
        
        %introduce sample arrays
        opt_data_det_ar = cell(1,n_sample);
        opt_data_red_ar = cell(1,n_sample);
        
        
        %load the reduced data and model
        load(fullfile(fdir,'coarse2','reduced_data_complete_N200.mat'));
        %load(fullfile(fdir,'coarse16','reduced_data_complete_coarse16.mat'));
        
        model.base_model.save_time_indices = 0:model.base_model.nt;
        model.base_model.base_model.save_time_indices = 0:model.base_model.nt;
        model.verbose = 4;
        model.base_model.verbose = 4;
        %set optimization parameters
        model.optimization.tol = 0.0005;
        model.base_model.optimization.tol = 0.0005;
        model.base_model.base_model.optimization.tol =0.0005;
        %set new fields in model
        model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.base_model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.optimization.stepsize_rule = 'armijo';
        model.optimization.sigma = 0.8;
        model.optimization.stepsize_factor = 0.75;
        model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.optimization.sigma = 0.8;
        model.base_model.optimization.stepsize_factor = 0.75;
        
        %set constants(replace later!!)
        model.gamma_H = @my_get_gamma_H_constant;
        model.base_model.gamma_H = @my_get_gamma_H_constant;
        model.base_model.base_model.gamma_H = @my_get_gamma_H_constant;
        model.L_H = @my_get_L_H_constant;
        model.base_model.L_H = @my_get_L_H_constant;
        model.base_model.base_model.L_H = @my_get_L_H_constant;
        
        %generate model_data
        model_data = gen_model_data(model);
        
        %generating the sample set
        M = rand_uniform(n_sample,model.mu_ranges);
        
       
        for i=1:n_sample
            disp('===================')
            disp(['Testpoint ',num2str(i),' of ',num2str(n_sample)])
            disp(['initial parameters: ',num2str(M(:,i)')])
            disp('----------------')
            %set new initial parameters
            buf_model = model;
            buf_model.optimization.init_params = M(:,i);
            buf_model.base_model.optimization.init_params = M(:,i);
            buf_model.base_model.base_model.optimization.init_params = M(:,i);
            
            
            %optimization detailed
            disp('detailed optimization...')
            tic;
            opt_data_det = optimize(buf_model, model_data);
            t_opt_det = toc;
            opt_data_det.t_opt = t_opt_det;
            %optimization reduced
            disp('reduced optimization')
            tic
            opt_data_red = optimize(buf_model, reduced_data);
            t_opt = toc;
            opt_data_red.t_opt = t_opt;
            
            opt_data_det_ar{i} = opt_data_det;
            opt_data_red_ar{i} = opt_data_red;  
        end
        
        
        %evaluate results
        nRB_samples = [];
        nRB_der_samples = [];
        t_opt_red_samples = zeros(n_sample,1);
        t_opt_det_samples = zeros(n_sample,1);
        Delta_samples = zeros(n_sample,1);
        err_samples = zeros(n_sample,1);
        effectivity_samples = zeros(n_sample,1);
        opti_steps_samples = zeros(n_sample,1);
        nr_fct_calls_samples = zeros(n_sample,1);
        nr_grad_calls_samples = zeros(n_sample,1);
        for i = 1:n_sample
            nRB_samples = [nRB_samples, opt_data_red_ar{i}.nRB];
            nRB_der_samples = [nRB_der_samples, opt_data_red_ar{i}.nRB_der];
            t_opt_red_samples(i) = opt_data_red_ar{i}.t_opt;
            t_opt_det_samples(i) = opt_data_det_ar{i}.t_opt;
            Delta_samples(i) = opt_data_red_ar{i}.Delta_mu(end);
            err_samples(i) = norm(opt_data_det_ar{i}.optimal_params - opt_data_red_ar{i}.optimal_params);
            effectivity_samples(i) = Delta_samples(i) / err_samples(i);
            opti_steps_samples(i) = opt_data_red_ar{i}.nr_opti_iter;
            nr_fct_calls_samples(i) = opt_data_red_ar{i}.nr_fct_calls;
            nr_grad_calls_samples(i) = opt_data_red_ar{i}.nr_grad_calc;
        end
        
        results.M = M;
        results.nRB_max = max(nRB_samples);
        results.nRB_min = min(nRB_samples);
        results.nRB_av = mean(nRB_samples);
        results.nRB_der_max = max(nRB_der_samples');
        results.nRB_der_min = min(nRB_der_samples');
        results.nRB_der_av = mean(nRB_der_samples');
        results.t_opt_det_max = max(t_opt_det_samples);
        results.t_opt_det_min = min(t_opt_det_samples);
        results.t_opt_det_av = mean(t_opt_det_samples);
        results.t_opt_red_max = max(t_opt_red_samples);
        results.t_opt_red_min = min(t_opt_red_samples);
        results.t_opt_red_av = mean(t_opt_red_samples);
        results.Delta_max = max(Delta_samples);
        results.Delta_min = min(Delta_samples);
        results.Delta_av = mean(Delta_samples);
        results.err_max = max(err_samples);
        results.err_min = min(err_samples);
        results.err_av = mean(err_samples);
        results.effectivity_max = max(effectivity_samples);
        results.effectivity_min = min(effectivity_samples);
        results.effectivity_av = mean(effectivity_samples);
        results.opti_steps_max = max(opti_steps_samples);
        results.opti_steps_min = min(opti_steps_samples);
        results.opti_steps_av = mean(opti_steps_samples);
        results.nr_fct_calls_max = max(nr_fct_calls_samples);
        results.nr_fct_calls_min = min(nr_fct_calls_samples);
        results.nr_fct_calls_av = mean(nr_fct_calls_samples);
        results.nr_grad_calls_max = max(nr_grad_calls_samples);
        results.nr_grad_calls_min = min(nr_grad_calls_samples);
        results.nr_grad_calls_av = mean(nr_grad_calls_samples);
        
        name = ['results_case_15_coarse',num2str(model.base_model.coarse_factor),...
            '_nsample_',num2str(n_sample),'.mat'];
        save(fullfile(fdir,name),'results','opt_data_red_ar','opt_data_det_ar');
        
        %email senden
        setpref('Internet','SMTP_Server','wap04.mathematik.uni-stuttgart.de')
        sendmail('dihlmann','case15 fertig',['Initial value changes with n_sample=',num2str(n_sample),' for functional 1 finished. Good work Spock!']);
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% weitere auswertungen der in case 15 erhaltenen Ergebnisse
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

    case 16
        
        load('reusults_case_15_coarse16_nsample_3.mat')
        
        opt_params=[];
        
        for i=1:length(opt_data_det_ar)
            opt_params = [opt_params;opt_data_det_ar{i}.optimal_params];
        end
        
        scatter(opt_params(:,1),opt_params(:,2));
        
        keyboard;
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
% performing a detailled and a reduced optimization 
% for various initial prameters (randomly chosen)
% the evaluating the results
% functional 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 17
        
        %number of samples
        if (nargin>1)&&(isfield(expar,'n_sample'))
            n_sample = expar.n_sample;
        else
            n_sample = 3;
        end
        
        %introduce sample arrays
        opt_data_det_ar = cell(1,n_sample);
        opt_data_red_ar = cell(1,n_sample);
                
        %set optimization parameters
        load(fullfile(fdir,'coarse2','basisNmax200','reduced_data_complete_2ndfunctional.mat'));
       
        %set 2nd functional fields
        model.base_model.save_time_indices = 0:model.base_model.nt;
        model.base_model.base_model.save_time_indices = 0:model.base_model.nt;
        model.verbose = 4;
        model.base_model.verbose = 4;
        model.name_output_functional = 'box_mean_time';
        model.get_output = @lin_evol_opt_output_time_integrated;
        model.output_function_ptr = @output_function_box_mean;
        model.sbox_xmin = 0.5;
        model.sbox_xmax = 1;
        model.sbox_ymin = 0.5;
        model.sbox_ymax = 1;
        model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.sbox_xmin = 0.5;
        model.base_model.sbox_xmax = 1;
        model.base_model.sbox_ymin = 0.5;
        model.base_model.sbox_ymax = 1;
        model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.base_model.sbox_xmin = 0.5;
        model.base_model.base_model.sbox_xmax = 1;
        model.base_model.base_model.sbox_ymin = 0.5;
        model.base_model.base_model.sbox_ymax = 1;
        model.base_model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.base_model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.base_model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        %set optimzation parameters
        model.optimization.tol = 0.0005;
        model.base_model.optimization.tol = 0.0005;
        model.base_model.base_model.optimization.tol =0.0005;
        
        model.optimization.init_params =[0.3,0.3];
        model.base_model.optimization.init_params =[0.3,0.3];
        model.base_model.base_model.optimization.init_params =[0.3,0.3];
        
        model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.base_model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        
        model.optimization.stepsize_rule = 'armijo';
        model.optimization.sigma = 0.8;
        model.optimization.stepsize_factor = 0.75;
        model.optimization.initial_stepsize = 10;
        model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.optimization.sigma = 0.8;
        model.base_model.optimization.stepsize_factor = 0.75;
        model.base_model.optimization.initial_stepsize = 10;
        model.base_model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.base_model.optimization.sigma = 0.8;
        model.base_model.base_model.optimization.stepsize_factor = 0.75;
        model.base_model.base_model.optimization.initial_stepsize = 10;
        
        model.gamma_H = @(model)dummy_constant(model,9.4347);
        model.base_model.gamma_H = @(model)dummy_constant(model,9.4347);
        model.base_model.base_model.gamma_H =@(model)dummy_constant(model,9.4347);
        model.L_H = @(model)dummy_constant(model,3.442);
        model.base_model.L_H =@(model)dummy_constant(model,3.442);
        model.base_model.base_model.L_H = @(model)dummy_constant(model,3.442);
       
        %generate model_data
        model_data = gen_model_data(model);
        
        %generating the sample set
        M = rand_uniform(n_sample,model.mu_ranges);
        
        for i=1:n_sample
            disp('===================')
            disp(['Testpoint ',num2str(i),' of ',num2str(n_sample)])
            disp('----------------')
            %set new initial parameters
            buf_model = model;
            buf_model.optimization.init_params = M(:,i);
            buf_model.base_model.optimization.init_params = M(:,i);
            buf_model.base_model.base_model.optimization.init_params = M(:,i);
            
            
            %optimization detailed
            tic;
            opt_data_det = optimize(buf_model, model_data);
            t_opt_det = toc;
            opt_data_det.t_opt = t_opt_det;
            %optimization reduced
            tic
            opt_data_red = optimize(buf_model, reduced_data);
            t_opt = toc;
            opt_data_red.t_opt = t_opt;
            
            opt_data_det_ar{i} = opt_data_det;
            opt_data_red_ar{i} = opt_data_red;  
        end
        
        %evaluate results
        nRB_samples = [];
        nRB_der_samples = [];
        t_opt_red_samples = zeros(n_sample,1);
        t_opt_det_samples = zeros(n_sample,1);
        Delta_samples = zeros(n_sample,1);
        err_samples = zeros(n_sample,1);
        effectivity_samples = zeros(n_sample,1);
        opti_steps_samples = zeros(n_sample,1);
        nr_fct_calls_samples = zeros(n_sample,1);
        nr_grad_calls_samples = zeros(n_sample,1);
        for i = 1:n_sample
            nRB_samples = [nRB_samples, opt_data_red_ar{i}.nRB];
            nRB_der_samples = [nRB_der_samples, opt_data_red_ar{i}.nRB_der];
            t_opt_red_samples(i) = opt_data_red_ar{i}.t_opt;
            t_opt_det_samples(i) = opt_data_det_ar{i}.t_opt;
            Delta_samples(i) = opt_data_red_ar{i}.Delta_mu(end);
            err_samples(i) = norm(opt_data_det_ar{i}.optimal_params - opt_data_red_ar{i}.optimal_params);
            effectivity_samples(i) = Delta_samples(i) / err_samples(i);
            opti_steps_samples(i) = opt_data_red_ar{i}.nr_opti_iter;
            nr_fct_calls_samples(i) = opt_data_red_ar{i}.nr_fct_calls;
            nr_grad_calls_samples(i) = opt_data_red_ar{i}.nr_grad_calc;
        end
        
        results.M = M;
        results.nRB_max = max(nRB_samples);
        results.nRB_min = min(nRB_samples);
        results.nRB_av = mean(nRB_samples);
        results.nRB_der_max = max(nRB_der_samples');
        results.nRB_der_min = min(nRB_der_samples');
        results.nRB_der_av = mean(nRB_der_samples');
        results.t_opt_det_max = max(t_opt_det_samples);
        results.t_opt_det_min = min(t_opt_det_samples);
        results.t_opt_det_av = mean(t_opt_det_samples);
        results.t_opt_red_max = max(t_opt_red_samples);
        results.t_opt_red_min = min(t_opt_red_samples);
        results.t_opt_red_av = mean(t_opt_red_samples);
        results.Delta_max = max(Delta_samples);
        results.Delta_min = min(Delta_samples);
        results.Delta_av = mean(Delta_samples);
        results.err_max = max(err_samples);
        results.err_min = min(err_samples);
        results.err_av = mean(err_samples);
        results.effectivity_max = max(effectivity_samples);
        results.effectivity_min = min(effectivity_samples);
        results.effectivity_av = mean(effectivity_samples);
        results.opti_steps_max = max(opti_steps_samples);
        results.opti_steps_min = min(opti_steps_samples);
        results.opti_steps_av = mean(opti_steps_samples);
        results.nr_fct_calls_max = max(nr_fct_calls_samples);
        results.nr_fct_calls_min = min(nr_fct_calls_samples);
        results.nr_fct_calls_av = mean(nr_fct_calls_samples);
        results.nr_grad_calls_max = max(nr_grad_calls_samples);
        results.nr_grad_calls_min = min(nr_grad_calls_samples);
        results.nr_grad_calls_av = mean(nr_grad_calls_samples);
        
        name = ['results_case_17_coarse',num2str(model.base_model.coarse_factor),...
            '_nsample_',num2str(n_sample),'.mat'];
        save(fullfile(fdir,name),'results','opt_data_red_ar','opt_data_det_ar');
        
        
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%
%  Optimierung mit dem Simplex (vorprogrammiert in Matlab)
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 18
        
        %load the reduced data and model (later: outflow model)
        %load(fullfile(fdir,'coarse2','reduced_data_complete_N200.mat'));
        load(fullfile(fdir,'coarse16','reduced_data_complete_coarse16.mat'));
        model.verbose = 3;
        
        %set optimization parameters
        model.optimization.tol = 0.0005;
        model.base_model.optimization.tol = 0.0005;
        model.base_model.base_model.optimization.tol =0.0005;
        model.optimization.init_params = [0.3,0.3];
        model.base_model.optimization.init_params = [0.3,0.3];
        model.base_model.base_model.optimization.init_params = [0.3,0.3];
        
        %set new fields in model
        model.base_model.optimization.optimizer = @simplex_nonlinear;
        model.base_model.base_model.optimization.optimizer = @simplex_nonlinear;
        
        
        %set constants(replace later!!)
        model.gamma_H = @my_get_gamma_H_constant;
        model.base_model.gamma_H = @my_get_gamma_H_constant;
        model.base_model.base_model.gamma_H = @my_get_gamma_H_constant;
        model.L_H = @my_get_L_H_constant;
        model.base_model.L_H = @my_get_L_H_constant;
        model.base_model.base_model.L_H = @my_get_L_H_constant;
        
        model_data = gen_model_data(model);
        
        %optimization
        tic
        opt_data = optimize(model, reduced_data);
        t_opt = toc;
        opt_data.t_opt = t_opt;
        
        keyboard;
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lade die ergebnisse von case 15 und führe mit dem Simplex algorithmus für
% die gleichen Testpunkte auch eine Optimierung durch
% Functional 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 19
        
        
         if (nargin>1)&&(isfield(expar,'n_sample'))
            n_sample = expar.n_sample;
        else
            n_sample = 30;
        end
        
        
        % lade reduced data functional 1
        load(fullfile(fdir,'coarse2','reduced_data_complete_N200.mat'));
        %load(fullfile(fdir,'coarse16','reduced_data_complete_coarse16.mat'));
        
        % load results from case 15
%         name = ['reusults_case_15_coarse',num2str(model.base_model.coarse_factor),...
%             '_nsample_',num2str(n_sample),'.mat'];
%         load(name);
%         
        %initialize variables
        M = rand_uniform(n_sample,model.mu_ranges);
        %M = results.M;
        opt_data_simpl_det_ar = cell(1,n_sample);
        opt_data_simpl_red_ar = cell(1,n_sample);
        
        
        %set fields for functional 1 and simplex algorithm
        model.verbose = 3;
        %set optimization parameters
        model.optimization.tol = 0.0005;
        model.base_model.optimization.tol = 0.0005;
        model.base_model.base_model.optimization.tol =0.0005;
        model.base_model.optimization.optimizer = @simplex_nonlinear;
        model.base_model.base_model.optimization.optimizer = @simplex_nonlinear;
        
        %set constants(replace later!!)
        model.gamma_H = @my_get_gamma_H_constant;
        model.base_model.gamma_H = @my_get_gamma_H_constant;
        model.base_model.base_model.gamma_H = @my_get_gamma_H_constant;
        model.L_H = @my_get_L_H_constant;
        model.base_model.L_H = @my_get_L_H_constant;
        model.base_model.base_model.L_H = @my_get_L_H_constant;
        
        %generate model_data
        model_data = gen_model_data(model);
        
        
        %optimize for all test points
        for i=1:n_sample
            disp('===================')
            disp(['Testpoint ',num2str(i),' of ',num2str(n_sample)])
            disp('----------------')
            %set new initial parameters
            model.optimization.init_params = M(:,i);
            model.base_model.optimization.init_params = M(:,i);
            model.base_model.base_model.optimization.init_params = M(:,i);
            
            
            %optimization detailed
            disp('Detailed simplex optimization...')
            tic;
            opt_data_simpl_det = optimize(model, model_data);
            t_opt_det = toc;
            opt_data_simpl_det.t_opt = t_opt_det;
            %optimization reduced
            disp('reduced simplex optimization...')
            tic
            opt_data_simpl_red = optimize(model, reduced_data);
            t_opt = toc;
            opt_data_simpl_red.t_opt = t_opt;
            
            opt_data_simpl_det_ar{i} = opt_data_simpl_det;
            opt_data_simpl_red_ar{i} = opt_data_simpl_red;
        end
        
        
        
        
        
        %evaluate the resulte
        
        t_opt_red_samples = zeros(n_sample,1);
        t_opt_det_samples = zeros(n_sample,1);
        Delta_samples = zeros(n_sample,1);
        err_samples = zeros(n_sample,1);
        effectivity_samples = zeros(n_sample,1);
        opti_steps_samples = zeros(n_sample,1);
        nr_fct_calls_samples = zeros(n_sample,1);
        for i = 1:n_sample
            t_opt_red_samples(i) = opt_data_simpl_red_ar{i}.t_opt;
            t_opt_det_samples(i) = opt_data_simpl_det_ar{i}.t_opt;
            Delta_samples(i) = opt_data_simpl_red_ar{i}.Delta_mu(end);
            err_samples(i) = norm(opt_data_simpl_det_ar{i}.optimal_params - opt_data_simpl_red_ar{i}.optimal_params);
            effectivity_samples(i) = Delta_samples(i) / err_samples(i);
            opti_steps_samples(i) = opt_data_simpl_red_ar{i}.nr_iterations;
            nr_fct_calls_samples(i) = opt_data_simpl_red_ar{i}.nr_fct_calls;
            
        end
        
        results_simpl.M = M;
        results_simpl.t_opt_det_max = max(t_opt_det_samples);
        results_simpl.t_opt_det_min = min(t_opt_det_samples);
        results_simpl.t_opt_det_av = mean(t_opt_det_samples);
        results_simpl.t_opt_red_max = max(t_opt_red_samples);
        results_simpl.t_opt_red_min = min(t_opt_red_samples);
        results_simpl.t_opt_red_av = mean(t_opt_red_samples);
        results_simpl.Delta_max = max(Delta_samples);
        results_simpl.Delta_min = min(Delta_samples);
        results_simpl.Delta_av = mean(Delta_samples);
        results_simpl.err_max = max(err_samples);
        results_simpl.err_min = min(err_samples);
        results_simpl.err_av = mean(err_samples);
        results_simpl.effectivity_max = max(effectivity_samples);
        results_simpl.effectivity_min = min(effectivity_samples);
        results_simpl.effectivity_av = mean(effectivity_samples);
        results_simpl.opti_steps_max = max(opti_steps_samples);
        results_simpl.opti_steps_min = min(opti_steps_samples);
        results_simpl.opti_steps_av = mean(opti_steps_samples);
        results_simpl.nr_fct_calls_max = max(nr_fct_calls_samples);
        results_simpl.nr_fct_calls_min = min(nr_fct_calls_samples);
        results_simpl.nr_fct_calls_av = mean(nr_fct_calls_samples);
        
        
        %save the results
        name = ['results_case_19_simplex_coarse',num2str(model.base_model.coarse_factor),...
            '_nsample_',num2str(n_sample),'.mat'];
        save(fullfile(fdir,name),'results_simpl','opt_data_simpl_red_ar','opt_data_simpl_det_ar');
        
        %email senden
        setpref('Internet','SMTP_Server','wap04.mathematik.uni-stuttgart.de')
        sendmail('dihlmann','case19 fertig',['Simplex optimization with ',num2str(n_sample),' samples for functional 1 finished. Good work Spock!']);
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lade die ergebnisse von case 17 und führe mit dem Simplex algorithmus für
% die gleichen Testpunkte auch eine Optimierung durch
% Functional 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 20
        
        n_sample = 30;
        
        % lade reduced data functional 2
        load(fullfile(fdir,'coarse2','basisNmax200','reduced_data_complete_2ndfunctional.mat'));
        %load(fullfile(fdir,'coarse16','reduced_data_complete_coarse16_2ndfunctional.mat'));
        
        % load results from case 17
%         name = ['reusults_case_17_coarse',num2str(model.base_model.coarse_factor),...
%             '_nsample_',num2str(n_sample),'.mat'];
%         load(name);
        
        %initialize variables
        M = rand_uniform(n_sample,model.mu_ranges);
        %M = results.M;
        opt_data_simpl_det_ar = cell(1,n_sample);
        opt_data_simpl_red_ar = cell(1,n_sample);
        
        
        %set fields for functional 1 and simplex algorithm
         model.base_model.save_time_indices = 0:model.base_model.nt;
        model.base_model.base_model.save_time_indices = 0:model.base_model.nt;
 model.verbose = 4;
        model.base_model.verbose = 4;
        model.name_output_functional = 'box_mean_time';
        model.get_output = @lin_evol_opt_output_time_integrated;
        model.output_function_ptr = @output_function_box_mean;
        model.sbox_xmin = 0.5;
        model.sbox_xmax = 1;
        model.sbox_ymin = 0.5;
        model.sbox_ymax = 1;
        model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.sbox_xmin = 0.5;
        model.base_model.sbox_xmax = 1;
        model.base_model.sbox_ymin = 0.5;
        model.base_model.sbox_ymax = 1;
        model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.base_model.sbox_xmin = 0.5;
        model.base_model.base_model.sbox_xmax = 1;
        model.base_model.base_model.sbox_ymin = 0.5;
        model.base_model.base_model.sbox_ymax = 1;
        model.base_model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.base_model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.base_model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        %set optimzation parameters
        model.optimization.tol = 0.0005;
        model.base_model.optimization.tol = 0.0005;
        model.base_model.base_model.optimization.tol =0.0005;
        
        model.gamma_H = @(model)dummy_constant(model,9.4347);
        model.base_model.gamma_H = @(model)dummy_constant(model,9.4347);
        model.base_model.base_model.gamma_H =@(model)dummy_constant(model,9.4347);
        model.L_H = @(model)dummy_constant(model,3.442);
        model.base_model.L_H =@(model)dummy_constant(model,3.442);
        model.base_model.base_model.L_H = @(model)dummy_constant(model,3.442);
       
        model.base_model.optimization.optimizer = @simplex_nonlinear;
        model.base_model.base_model.optimization.optimizer = @simplex_nonlinear;
                
        
        %generate model_data
        model_data = gen_model_data(model);
        
        
        %optimize for all test points
        for i=1:n_sample
            disp('===================')
            disp(['Testpoint ',num2str(i),' of ',num2str(n_sample)])
            disp('----------------')
            %set new initial parameters
            model.optimization.init_params = M(:,i);
            model.base_model.optimization.init_params = M(:,i);
            model.base_model.base_model.optimization.init_params = M(:,i);
            
            
            %optimization detailed
            tic;
            opt_data_simpl_det = optimize(model, model_data);
            t_opt_det = toc;
            opt_data_simpl_det.t_opt = t_opt_det;
            %optimization reduced
            tic
            opt_data_simpl_red = optimize(model, reduced_data);
            t_opt = toc;
            opt_data_simpl_red.t_opt = t_opt;
            
            opt_data_simpl_det_ar{i} = opt_data_simpl_det;
            opt_data_simpl_red_ar{i} = opt_data_simpl_red;
        end
        
        
        %evaluate the resulte
        t_opt_red_samples = zeros(n_sample,1);
        t_opt_det_samples = zeros(n_sample,1);
        Delta_samples = zeros(n_sample,1);
        err_samples = zeros(n_sample,1);
        effectivity_samples = zeros(n_sample,1);
        opti_steps_samples = zeros(n_sample,1);
        nr_fct_calls_samples = zeros(n_sample,1);
        for i = 1:n_sample
            t_opt_red_samples(i) = opt_data_simpl_red_ar{i}.t_opt;
            t_opt_det_samples(i) = opt_data_simpl_det_ar{i}.t_opt;
            Delta_samples(i) = opt_data_simpl_red_ar{i}.Delta_mu(end);
            err_samples(i) = norm(opt_data_simpl_det_ar{i}.optimal_params - opt_data_simpl_red_ar{i}.optimal_params);
            effectivity_samples(i) = Delta_samples(i) / err_samples(i);
            opti_steps_samples(i) = opt_data_simpl_red_ar{i}.nr_iterations;
            nr_fct_calls_samples(i) = opt_data_simpl_red_ar{i}.nr_fct_calls;
            
        end
        
        results_simpl.M = M;
        results_simpl.t_opt_det_max = max(t_opt_det_samples);
        results_simpl.t_opt_det_min = min(t_opt_det_samples);
        results_simpl.t_opt_det_av = mean(t_opt_det_samples);
        results_simpl.t_opt_red_max = max(t_opt_red_samples);
        results_simpl.t_opt_red_min = min(t_opt_red_samples);
        results_simpl.t_opt_red_av = mean(t_opt_red_samples);
        results_simpl.Delta_max = max(Delta_samples);
        results_simpl.Delta_min = min(Delta_samples);
        results_simpl.Delta_av = mean(Delta_samples);
        results_simpl.err_max = max(err_samples);
        results_simpl.err_min = min(err_samples);
        results_simpl.err_av = mean(err_samples);
        results_simpl.effectivity_max = max(effectivity_samples);
        results_simpl.effectivity_min = min(effectivity_samples);
        results_simpl.effectivity_av = mean(effectivity_samples);
        results_simpl.opti_steps_max = max(opti_steps_samples);
        results_simpl.opti_steps_min = min(opti_steps_samples);
        results_simpl.opti_steps_av = mean(opti_steps_samples);
        results_simpl.nr_fct_calls_max = max(nr_fct_calls_samples);
        results_simpl.nr_fct_calls_min = min(nr_fct_calls_samples);
        results_simpl.nr_fct_calls_av = mean(nr_fct_calls_samples);
        
        
        %save the results
        name = ['results_case_20_simplex_coarse',num2str(model.base_model.coarse_factor),...
            '_nsample_',num2str(n_sample),'.mat'];
        save(fullfile(fdir,name),'results_simpl','opt_data_simpl_red_ar','opt_data_simpl_det_ar');
        

        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Experiment changing the breaking tolerance
% we chacnge the breaking tolerance for the gradient optimization to stop
% Observe: When is the breaking condition for the error estimator
% fullfilled, how does the error estimator change with the changing
% breaking condition, how does the effectivity change?
%-------
% 2nd functional
%---------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
    case 21
        
        %epsilon breaking values
        tol_vec = [0.005,0.002, 0.001, 0.0008,0.0005];
        n_sample = length(tol_vec);
        
        %introduce sample arrays
        opt_data_det_ar = cell(1,n_sample);
        opt_data_red_ar = cell(1,n_sample);
        
        %_________________________________________________________________________
        %functional 1
        %===========================
        %         %set optimization parameters
        %         load(fullfile(fdir,'coarse2','reduced_data_complete_N200.mat'));
        %         %load(fullfile(fdir,'coarse16','reduced_data_complete_coarse16.mat'));
        %
        %         model.verbose = 4;
        %         model.base_model.verbose = 4;
        %         %set optimization parameters
        %         model.optimization.tol = 0.0005;
        %         model.base_model.optimization.tol = 0.0005;
        %         model.base_model.base_model.optimization.tol =0.0005;
        %         %set new fields in model
        %         model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        %         model.base_model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        %         model.optimization.stepsize_rule = 'armijo';
        %         model.optimization.sigma = 0.8;
        %         model.optimization.stepsize_factor = 0.75;
        %         model.base_model.optimization.stepsize_rule = 'armijo';
        %         model.base_model.optimization.sigma = 0.8;
        %         model.base_model.optimization.stepsize_factor = 0.75;
        %
        %         %set constants(replace later!!)
        %         model.gamma_H = @my_get_gamma_H_constant;
        %         model.base_model.gamma_H = @my_get_gamma_H_constant;
        %         model.base_model.base_model.gamma_H = @my_get_gamma_H_constant;
        %         model.L_H = @my_get_L_H_constant;
        %         model.base_model.L_H = @my_get_L_H_constant;
        %         model.base_model.base_model.L_H = @my_get_L_H_constant;
        %
        
        %________________________________________________________________________________
        % functional 2
        %=====================
        
        %set optimization parameters
        load(fullfile(fdir,'coarse2','basisNmax200','reduced_data_complete_2ndfunctional.mat'));
        
        %set 2nd functional fields
        model.base_model.save_time_indices = 0:model.base_model.nt;
        model.base_model.base_model.save_time_indices = 0:model.base_model.nt;
        model.verbose = 4;
        model.base_model.verbose = 4;
        model.name_output_functional = 'box_mean_time';
        model.get_output = @lin_evol_opt_output_time_integrated;
        model.output_function_ptr = @output_function_box_mean;
        model.sbox_xmin = 0.5;
        model.sbox_xmax = 1;
        model.sbox_ymin = 0.5;
        model.sbox_ymax = 1;
        model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.sbox_xmin = 0.5;
        model.base_model.sbox_xmax = 1;
        model.base_model.sbox_ymin = 0.5;
        model.base_model.sbox_ymax = 1;
        model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.base_model.sbox_xmin = 0.5;
        model.base_model.base_model.sbox_xmax = 1;
        model.base_model.base_model.sbox_ymin = 0.5;
        model.base_model.base_model.sbox_ymax = 1;
        model.base_model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.base_model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.base_model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        %set optimzation parameters
        
        model.optimization.init_params =[0.35,0.35];
        model.base_model.optimization.init_params =[0.35,0.35];
        model.base_model.base_model.optimization.init_params =[0.35,0.35];
        
        model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.base_model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        
        model.optimization.stepsize_rule = 'armijo';
        model.optimization.sigma = 0.8;
        model.optimization.stepsize_factor = 0.75;
        model.optimization.initial_stepsize = 10;
        model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.optimization.sigma = 0.8;
        model.base_model.optimization.stepsize_factor = 0.75;
        model.base_model.optimization.initial_stepsize = 10;
        model.base_model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.base_model.optimization.sigma = 0.8;
        model.base_model.base_model.optimization.stepsize_factor = 0.75;
        model.base_model.base_model.optimization.initial_stepsize = 10;
        
        model.gamma_H = @(model)dummy_constant(model,100);
        model.base_model.gamma_H = @(model)dummy_constant(model,100);
        model.base_model.base_model.gamma_H =@(model)dummy_constant(model,100);
        model.L_H = @(model)dummy_constant(model,100);
        model.base_model.L_H =@(model)dummy_constant(model,100);
        model.base_model.base_model.L_H = @(model)dummy_constant(model,100);
        
        model.optimization.allow_constant_calculation = 1;
        model.base_model.optimization.allow_constant_calculation = 1;
        model.base_model.base_model.optimization.allow_constant_calculation = 1;
        
        
        %generate model_data
        model_data = gen_model_data(model);
        
        
        for i=1:n_sample
            
            disp('===================')
            disp(['Testpoint ',num2str(i),' of ',num2str(n_sample)])
            disp('----------------')
            %set new tolerance
            model.optimization.tol = tol_vec(i);
            model.base_model.optimization.tol = tol_vec(i);
            model.base_model.base_model.optimization.tol =tol_vec(i);
            
            %optimization detailed
            disp('detailed optimization ')
            disp('---------------------')
            tic;
            opt_data_det = optimize(model, model_data);
            t_opt_det = toc;
            opt_data_det.t_opt = t_opt_det;
            %optimization reduced
            disp('reduced optimization ')
            disp('---------------------')
            tic
            opt_data_red = optimize(model, reduced_data);
            t_opt = toc;
            opt_data_red.t_opt = t_opt;
            
            opt_data_det_ar{i} = opt_data_det;
            opt_data_red_ar{i} = opt_data_red;
        end
        
        %evaluate results
        nRB_samples = [];
        nRB_der_samples = [];
        t_opt_red_samples = zeros(n_sample,1);
        t_opt_det_samples = zeros(n_sample,1);
        Delta_samples = zeros(n_sample,1);
        err_samples = zeros(n_sample,1);
        effectivity_samples = zeros(n_sample,1);
        opti_steps_samples = zeros(n_sample,1);
        nr_fct_calls_samples = zeros(n_sample,1);
        nr_grad_calls_samples = zeros(n_sample,1);
        break_value = zeros(1,n_sample);
        for i = 1:n_sample
            nRB_samples = [nRB_samples, opt_data_red_ar{i}.nRB];
            nRB_der_samples = [nRB_der_samples, opt_data_red_ar{i}.nRB_der];
            t_opt_red_samples(i) = opt_data_red_ar{i}.t_opt;
            t_opt_det_samples(i) = opt_data_det_ar{i}.t_opt;
            Delta_samples(i) = opt_data_red_ar{i}.Delta_mu(end);
            err_samples(i) = norm(opt_data_det_ar{i}.optimal_params - opt_data_red_ar{i}.optimal_params);
            effectivity_samples(i) = Delta_samples(i) / err_samples(i);
            opti_steps_samples(i) = opt_data_red_ar{i}.nr_opti_iter;
            nr_fct_calls_samples(i) = opt_data_red_ar{i}.nr_fct_calls;
            nr_grad_calls_samples(i) = opt_data_red_ar{i}.nr_grad_calc;
            break_value(i) = opt_data_red_ar{i}.break_value;
        end
        
        
        results.nRB_max = max(nRB_samples);
        results.nRB_min = min(nRB_samples);
        results.nRB_av = mean(nRB_samples);
        results.nRB_der_max = max(nRB_der_samples');
        results.nRB_der_min = min(nRB_der_samples');
        results.nRB_der_av = mean(nRB_der_samples');
        results.t_opt_det_max = max(t_opt_det_samples);
        results.t_opt_det_min = min(t_opt_det_samples);
        results.t_opt_det_av = mean(t_opt_det_samples);
        results.t_opt_red_max = max(t_opt_red_samples);
        results.t_opt_red_min = min(t_opt_red_samples);
        results.t_opt_red_av = mean(t_opt_red_samples);
        results.Delta_max = max(Delta_samples);
        results.Delta_min = min(Delta_samples);
        results.Delta_av = mean(Delta_samples);
        results.err_max = max(err_samples);
        results.err_min = min(err_samples);
        results.err_av = mean(err_samples);
        results.effectivity_max = max(effectivity_samples);
        results.effectivity_min = min(effectivity_samples);
        results.effectivity_av = mean(effectivity_samples);
        results.opti_steps_max = max(opti_steps_samples);
        results.opti_steps_min = min(opti_steps_samples);
        results.opti_steps_av = mean(opti_steps_samples);
        results.nr_fct_calls_max = max(nr_fct_calls_samples);
        results.nr_fct_calls_min = min(nr_fct_calls_samples);
        results.nr_fct_calls_av = mean(nr_fct_calls_samples);
        results.nr_grad_calls_max = max(nr_grad_calls_samples);
        results.nr_grad_calls_min = min(nr_grad_calls_samples);
        results.nr_grad_calls_av = mean(nr_grad_calls_samples);
        results.break_value = break_value;
        
        name = ['results_case_21_tolvariation_func2_coarse',num2str(model.base_model.coarse_factor),...
            '_nsample_',num2str(n_sample),'.mat'];
        save(fullfile(fdir,name),'results','opt_data_red_ar','opt_data_det_ar');
        
        
        %email senden
        setpref('Internet','SMTP_Server','wap04.mathematik.uni-stuttgart.de')
        sendmail('dihlmann','case21 fertig','Sich verändernde toleranzen fertig. Functional1');
        
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % changing the functionals: we change the position of the time interval
        % where the average concentration is calculated. Thereby we obtain
        % (hopefully) steep an flat functionals.
        % See if the breaking condition can be fullfilled for all type of
        % functionals.
        % Plan: change functional --> conduct optimization --> calculate constants
        % --> examine braking value
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
    case 22
        
        %set time intervals
        t1_vec = [0.75,0.5];%,0.25,0,0];
        t2_vec = [1,0.75];%,0.5,0.25,0.1];
        
        n_samples = length(t1_vec);
        opt_data_ar = cell(n_samples,1);
        %set optimization parameters
        load(fullfile(fdir,'coarse2','basisNmax200','reduced_data_complete_2ndfunctional.mat'));
        %load(fullfile(fdir,'coarse16','reduced_data_complete_2ndfunctional_coarse16.mat'));
        
        %set functional fields
        model.verbose = 4;
        model.base_model.verbose = 4;
        model.name_output_functional = 'box_mean_time';
        model.get_output = @lin_evol_opt_output_time_integrated;
        model.output_function_ptr = @output_function_box_mean;
        model.sbox_xmin = 0.5;
        model.sbox_xmax = 1;
        model.sbox_ymin = 0.5;
        model.sbox_ymax = 1;
        
        model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.sbox_xmin = 0.5;
        model.base_model.sbox_xmax = 1;
        model.base_model.sbox_ymin = 0.5;
        model.base_model.sbox_ymax = 1;
        
        model.base_model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.base_model.sbox_xmin = 0.5;
        model.base_model.base_model.sbox_xmax = 1;
        model.base_model.base_model.sbox_ymin = 0.5;
        model.base_model.base_model.sbox_ymax = 1;
        
        
        
        %set optimzation parameters
        model.optimization.tol = 0.0005;
        model.base_model.optimization.tol = 0.0005;
        model.base_model.base_model.optimization.tol =0.0005;
        
        model.optimization.init_params =[0.5,0.5];
        model.base_model.optimization.init_params =[0.5,0.5];
        model.base_model.base_model.optimization.init_params =[0.5,0.5];
        
        model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.base_model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        
        model.optimization.stepsize_rule = 'armijo';
        model.optimization.sigma = 0.8;
        model.optimization.stepsize_factor = 0.75;
        model.optimization.initial_stepsize = 10;
        model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.optimization.sigma = 0.8;
        model.base_model.optimization.stepsize_factor = 0.75;
        model.base_model.optimization.initial_stepsize = 10;
        model.base_model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.base_model.optimization.sigma = 0.8;
        model.base_model.base_model.optimization.stepsize_factor = 0.75;
        model.base_model.base_model.optimization.initial_stepsize = 10;
        
        %set constants very high so that a detailed constnat calculation is
        %triggered at the end of the optimization
        model.gamma_H = @(model)dummy_constant(model,100);
        model.base_model.gamma_H = @(model)dummy_constant(model,100);
        model.base_model.base_model.gamma_H =@(model)dummy_constant(model,100);
        model.L_H = @(model)dummy_constant(model,100);
        model.base_model.L_H =@(model)dummy_constant(model,100);
        model.base_model.base_model.L_H = @(model)dummy_constant(model,100);
        
        model.optimization.allow_constant_calculation = 1;
        model.base_model.optimization.allow_constant_calculation = 1;
        model.base_model.base_model.optimization.allow_constant_calculation = 1;
        
        
        for i=1:n_samples
            disp('-------------------')
            disp(['Sample nr',num2str(i),' of ',num2str(n_samples)])
            disp('==============================')
            %time frame settings
            model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt*t1_vec(i)));
            model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*t2_vec(i)));
            model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
            model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
            model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt*t1_vec(i)));
            model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*t2_vec(i)));
            model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
            model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
            model.base_model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt*t1_vec(i)));
            model.base_model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*t2_vec(i)));
            model.base_model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
            model.base_model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
            
            opt_data_ar{i} = optimize(model, reduced_data);
            
            
            
        end
        
        results.n_samples = n_samples;
        results.t1_vec = t1_vec;
        results.t2_vec = t2_vec;
        
        %save results
        name = ['results_case_22_functionalchange_coarse',num2str(model.base_model.coarse_factor),...
            '_nsample_',num2str(n_samples),'.mat'];
        save(fullfile(fdir,name),'results','opt_data_ar');
        
        
        %email senden
        setpref('Internet','SMTP_Server','wap04.mathematik.uni-stuttgart.de')
        sendmail('dihlmann','case22 fertig','aenderung von zeitintervallen');
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
        % vergleich von iterativem und nichtiterativem Fehlerschätzer
        %
        % von einem festen Stratpunkt aus werden ohne schrittweitensteuerung
        % gradientenoptimierungsverfahren durchgeführt.
        % Die werte beider Fehlerschätzer bei jedem Gradientenschritt werden
        % verglichen.
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
    case 23
        
        %lade model und reduzierte daten
        load(fullfile(fdir,'coarse2','reduced_data_complete_N200.mat'));
        
        %set optimization parameters
        model.optimization.init_params = [0.5,0.5];
        model.optimization.opt_mode = 'reduced';
        model.base_model.optimization.init_params = [0.5,0.5];
        model.base_model.optimization.opt_mode = 'reduced';
        model.base_model.base_model.optimization.init_params = [0.5,0.5];
        model.base_model.base_model.optimization.opt_mode = 'reduced';
        
        model.optimization.tol = 0.001;
        model.base_model.optimization.tol = 0.001;
        model.base_model.base_model.optimization.tol =0.001;
        
        model.base_model.optimization.optimizer = @gradient_opt;
        model.base_model.base_model.optimization.optimizer = @gradient_opt;
        model.optimization.stepsize_rule = 'quotient';
        model.base_model.optimization.stepsize_rule = 'quotient';
        model.base_model.base_model.optimization.stepsize_rule = 'quotient';
        model.optimization.initial_stepsize = 1;
        model.base_model.optimization.initial_stepsize = 1;
        model.base_model.base_model.optimization.initial_stepsize = 1;
        model.stepsize_coefficient = 3;
        model.base_model.stepsize_coefficient = 3;
        model.base_model.base_model.stepsize_coefficient = 3;
        
        
        %output fields
         model.output_function_ptr = @output_function_box_mean;
        model.sbox_xmin = 1;
        model.sbox_xmax = 2;
        model.sbox_ymin = 0;
        model.sbox_ymax = 0.5;
        model.s_l2norm = sqrt(2);%1/sqrt((model.sbox_ymax-model.sbox_ymin)*(model.sbox_xmax-model.sbox_xmin));%0.0078*coarse_factor;%analytisch: <=1/sqrt(0.5);
        model.name_output_functional = 'box_mean';
        model.get_output = @lin_evol_opt_output_time_independent;
        
        model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.sbox_xmin = 1;
        model.base_model.sbox_xmax = 2;
        model.base_model.sbox_ymin = 0;
        model.base_model.sbox_ymax = 0.5;
        model.base_model.s_l2norm = sqrt(2);%1/sqrt((model.sbox_ymax-model.sbox_ymin)*(model.sbox_xmax-model.sbox_xmin));%0.0078*coarse_factor;%analytisch: <=1/sqrt(0.5);
        model.base_model.name_output_functional = 'box_mean';
        model.base_model.get_output = @lin_evol_opt_output_time_independent;
        
        model.base_model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.base_model.sbox_xmin = 1;
        model.base_model.base_model.sbox_xmax = 2;
        model.base_model.base_model.sbox_ymin = 0;
        model.base_model.base_model.sbox_ymax = 0.5;
        model.base_model.base_model.s_l2norm = sqrt(2);%1/sqrt((model.sbox_ymax-model.sbox_ymin)*(model.sbox_xmax-model.sbox_xmin));%0.0078*coarse_factor;%analytisch: <=1/sqrt(0.5);
        model.base_model.base_model.name_output_functional = 'box_mean';
        model.base_model.base_model.get_output = @lin_evol_opt_output_time_independent;
        
        
        
        
        %conduct optimization with iterative error estimator
        opt_data_iter = optimize(model,reduced_data);
        
        %conduct optimzation wit noniterative error estimator
        model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.base_model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.gamma_H = @calculate_gamma_const;
        model.base_model.gamma_H = @calculate_gamma_const;
        model.base_model.base_model.gamma_H = @calculate_gamma_const;
        model.L_H = @calculate_Lipschitz_const;
        model.base_model.L_H = @calculate_Lipschitz_const;
        model.base_model.base_model.L_H = @calculate_Lipschitz_const;
%              model.gamma_H = @my_get_gamma_H_constant;
%                  model.base_model.gamma_H = @my_get_gamma_H_constant;
%                  model.base_model.base_model.gamma_H = @my_get_gamma_H_constant;
%                  model.L_H = @my_get_L_H_constant;
%                  model.base_model.L_H = @my_get_L_H_constant;
%                  model.base_model.base_model.L_H = @my_get_L_H_constant;
%         
        
        
        opt_data_noniter = optimize(model, reduced_data);
        
        save(fullfile(fdir,'results_case_23.mat'),'opt_data_iter','opt_data_noniter');
        
        
        %compare the results with case 51

        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%like case23: comparison between iterative and noniterative error estimator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
    case 24
        
        %set optimization parameters
        load(fullfile(fdir,'coarse2','basisNmax200','reduced_data_complete_2ndfunctional.mat'));
       %load(fullfile(fdir,'coarse2','reduced_data_complete_2ndfunctional.mat'));
       
        %set 2nd functional fields
        model.base_model.save_time_indices = 0:model.base_model.nt;
        model.base_model.base_model.save_time_indices = 0:model.base_model.nt;
        model.verbose = 4;
        model.base_model.verbose = 4;
        model.name_output_functional = 'box_mean_time';
        model.get_output = @lin_evol_opt_output_time_integrated;
        model.output_function_ptr = @output_function_box_mean;
        model.sbox_xmin = 0.5;
        model.sbox_xmax = 1;
        model.sbox_ymin = 0.5;
        model.sbox_ymax = 1;
        model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.sbox_xmin = 0.5;
        model.base_model.sbox_xmax = 1;
        model.base_model.sbox_ymin = 0.5;
        model.base_model.sbox_ymax = 1;
        model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.base_model.sbox_xmin = 0.5;
        model.base_model.base_model.sbox_xmax = 1;
        model.base_model.base_model.sbox_ymin = 0.5;
        model.base_model.base_model.sbox_ymax = 1;
        model.base_model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.base_model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.base_model.base_model.s_l2norm = sqrt(4)/sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.base_model.s_l1norm = sqrt(4)/((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        
        %set optimzation parameters
        model.optimization.tol = 0.0005;
        model.base_model.optimization.tol = 0.0005;
        model.base_model.base_model.optimization.tol =0.0005;
        
        model.optimization.init_params = [0.2,0.2];
        model.optimization.opt_mode = 'reduced';
        model.base_model.optimization.init_params = [0.2,0.2];
        model.base_model.optimization.opt_mode = 'reduced';
        model.base_model.base_model.optimization.init_params = [0.2,0.2];
        model.base_model.base_model.optimization.opt_mode = 'reduced';
        
         model.base_model.optimization.optimizer = @gradient_opt;
        model.base_model.base_model.optimization.optimizer = @gradient_opt;
        model.optimization.stepsize_rule = 'quotient';        
        model.base_model.optimization.stepsize_rule = 'quotient';
        model.base_model.base_model.optimization.stepsize_rule = 'quotient';
        model.optimization.initial_stepsize = 2;
        model.base_model.optimization.initial_stepsize = 2;
        model.base_model.base_model.optimization.initial_stepsize = 2;
        model.stepsize_coefficient = 100;
        model.base_model.stepsize_coefficient = 100;
        model.base_model.base_model.stepsize_coefficient = 100;
        
         %conduct optimization with iterative error estimator
        opt_data_iter = optimize(model,reduced_data);
        
        %conduct optimization with noniterative error estimator
        model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.base_model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
       model.gamma_H = @calculate_gamma_const;
        model.base_model.gamma_H = @calculate_gamma_const;
        model.base_model.base_model.gamma_H = @calculate_gamma_const;
        model.L_H = @calculate_Lipschitz_const;
        model.base_model.L_H = @calculate_Lipschitz_const;
        model.base_model.base_model.L_H = @calculate_Lipschitz_const;
        

        
        opt_data_noniter = optimize(model, reduced_data);
        
        save(fullfile(fdir,'results_case_24.mat'),'opt_data_iter','opt_data_noniter');
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% plot the funcitonal surface / Landscape generated with case 123 and 124
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 50
        load('Landschaft_func1_coarse2_numpts_25.mat')
        h=figure;
        surf(x_vec,y_vec,func1_det,'FaceColor','interp');
        xlabel('\mu_1');
        ylabel('\mu_2');
        zlabel('J(u)');
        print(h,'-dpng','Functional1.png');
        print(h,'-deps','Functional1.eps');
        
        load('Landschaft_func2_coarse2_numpts_25.mat')
        h=figure;
        surf(x_vec,y_vec,func2_det,'FaceColor','interp');
        xlabel('\mu_1');
        ylabel('\mu_2');
        zlabel('J(u)');
        print(h,'-dpng','Functional2.png');
        print(h,'-deps','Functional2.eps');
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5        
% evaluate and plot results fom case51:
% comparison of iterative and noniterative error estimator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 51
        
        %load result file from case 23
        %load(fullfile(fdir,'results_case_24.mat'));
        load(fullfile(fdir,'results_case_23.mat'));
        
        %plot the gradient way
        figure
        plot_PM(opt_data_iter);
        D=[];
        for i=1:size(opt_data_iter.Delta_mu,2)
            D(i) = norm(opt_data_iter.Delta_mu(:,i));
        end
        
        %trennen in break value >1 und <=1
        non_val_indx = find(opt_data_noniter.break_value_vec>1);
        val_indx = find(opt_data_noniter.break_value_vec<=1);
        
        if ~isempty(non_val_indx)
            non_val_indx = [non_val_indx,val_indx(1)];
            non_val_Delta_mu = opt_data_noniter.Delta_mu_vec(non_val_indx);
            non_val_grad_norm = opt_data_iter.norm_grad_f_x(non_val_indx);
        end
               
        val_Delta_mu = opt_data_noniter.Delta_mu_vec(val_indx);
        val_grad_norm = opt_data_iter.norm_grad_f_x(val_indx);
            
        
        figure
        loglog(opt_data_iter.norm_grad_f_x,D,'bo-','LineWidth',2)
        hold on
        loglog(val_grad_norm,val_Delta_mu,'ro-','LineWidth',2)
        if ~isempty(non_val_indx)
            loglog(non_val_grad_norm,non_val_Delta_mu,'o--r','LineWidth',2);
        end
        %axis([3.5e-4, 0.13,0.35e-3,4])%func2
        axis([8e-4, 0.4,4e-4,1.6])%func1
        title('Comparison of both error estimators')
        legend('incremental error est.','general error est. (valid)','general error est. (not valid)')
        xlabel('functional gradient norm')
        ylabel('\Delta_{\mu}')
        
        
        
        %with stepnumbers
        nr_steps = length(D);
        figure
        semilogy(0:nr_steps-1,D,'LineWidth',2)
        hold on
        semilogy(val_indx-1,val_Delta_mu,'r','LineWidth',2)
        if ~isempty(non_val_indx)
            semilogy(non_val_indx-1,non_val_Delta_mu,'r--','LineWidth',2)
        end
        %axis([3.5e-4, 0.13,0.35e-3,4])%func2
        title('Comparison of both error estimators')
        legend('incremental error est.','non-incremental error est. valid','non-incremental error est. not valid')
        xlabel('step number')
        ylabel('\Delta_{\mu}')
        

        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
%        Evaluation of case 20: Simplex optmization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5        
    case 52
        
        load('results_case_20_simplex_coarse2_nsample_30.mat')
        
        %The invalid error estimations (Delta_mu =1) are still included in
        %the evaluation... so evaluate again the results neglecting the
        %invalid error estimators
        
        %evaluate the resulte
        t_opt_red_samples = [];
        t_opt_det_samples = [];
        Delta_samples = [];
        err_samples = [];
        effectivity_samples = [];
        opti_steps_samples = [];
        nr_fct_calls_samples = [];
        for i = 1:length(opt_data_simpl_red_ar)
            if (opt_data_simpl_red_ar{i}.Delta_mu ~=-1)
                t_opt_red_samples = [t_opt_red_samples,opt_data_simpl_red_ar{i}.t_opt];
                t_opt_det_samples = [t_opt_det_samples,opt_data_simpl_det_ar{i}.t_opt];
                Delta_samples = [Delta_samples ,opt_data_simpl_red_ar{i}.Delta_mu(end)];
                err_samples = [err_samples,norm(opt_data_simpl_det_ar{i}.optimal_params - opt_data_simpl_red_ar{i}.optimal_params)];
                effectivity_samples = [effectivity_samples,Delta_samples(end) / err_samples(end)];
                opti_steps_samples = [opti_steps_samples,opt_data_simpl_red_ar{i}.nr_iterations];
                nr_fct_calls_samples = [nr_fct_calls_samples,opt_data_simpl_red_ar{i}.nr_fct_calls];
            end
        end
        
        
        results_simpl.t_opt_det_max = max(t_opt_det_samples);
        results_simpl.t_opt_det_min = min(t_opt_det_samples);
        results_simpl.t_opt_det_av = mean(t_opt_det_samples);
        results_simpl.t_opt_red_max = max(t_opt_red_samples);
        results_simpl.t_opt_red_min = min(t_opt_red_samples);
        results_simpl.t_opt_red_av = mean(t_opt_red_samples);
        results_simpl.Delta_max = max(Delta_samples);
        results_simpl.Delta_min = min(Delta_samples);
        results_simpl.Delta_av = mean(Delta_samples);
        results_simpl.err_max = max(err_samples);
        results_simpl.err_min = min(err_samples);
        results_simpl.err_av = mean(err_samples);
        results_simpl.effectivity_max = max(effectivity_samples);
        results_simpl.effectivity_min = min(effectivity_samples);
        results_simpl.effectivity_av = mean(effectivity_samples);
        results_simpl.opti_steps_max = max(opti_steps_samples);
        results_simpl.opti_steps_min = min(opti_steps_samples);
        results_simpl.opti_steps_av = mean(opti_steps_samples);
        results_simpl.nr_fct_calls_max = max(nr_fct_calls_samples);
        results_simpl.nr_fct_calls_min = min(nr_fct_calls_samples);
        results_simpl.nr_fct_calls_av = mean(nr_fct_calls_samples);
        
        
        %save the results
        name = ['results_evaluated_case_20_simplex30.mat'];
        save(fullfile(fdir,name),'results_simpl','opt_data_simpl_red_ar','opt_data_simpl_det_ar');
        
        keyboard
%%%%%%%%%%%%%%%%%%%%%%%5    
%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Tests and verifications
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Test the error estimators
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 100
        
        params.coarse_factor = 8;
        pp.plot_function = @plot_element_data;

        model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
        model.compute_derivative_info = 0;
        model.k=0.25;
        model.t=0;

        model_data = gen_model_data(model);

        mu= get_mu(model)

        sim_data = detailed_simulation(model, model_data);

        detailed_data = model_data;
        detailed_data.RB = [];

        detailed_data.RB = model.rb_init_data_basis(model, detailed_data);


        %PCA
        PCAtr = PCA_fixspace(sim_data.U, detailed_data.RB,model_data.W);
        detailed_data.RB = [detailed_data.RB,PCAtr];


        %check orthogonality
        RB = detailed_data.RB;
        oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
        loop = 0;
        while (oc>1e-14)&&(loop<5)
            RB = orthonormalize_gram_schmidt(RB, model_data.W);
            oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
            loop = loop +1;
        end

        detailed_data.RB = RB;

        reduced_data = gen_reduced_data(model, detailed_data);

        rb_sim_data = rb_simulation(model, reduced_data);

        %model.rb_reconstruction = @rb_reconstruction_default;
        %rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data);
        rb_sim_data = rb_reconstruction_derivative(model, detailed_data, rb_sim_data);
        
        U_rb = rb_sim_data.U;
        
        disp(' ');
        disp(['Geschaetzter Fehler: ',num2str(rb_sim_data.Delta(end))]);
        err = fv_l2_error(U_rb(:,model.save_time_indices+1), sim_data.U,model_data.grid,[]);
        disp(['Tatsaechlicher Fehler: ',num2str(err(end))]);
        plot_data_sequence(model, model_data.grid, U_rb(:,model.save_time_indices+1)-sim_data.U,pp);
        %plot_data_sequence(model, model_data.grid,U_rb-sim_data.U,plot_params);
        title('Fehlerplot');
        figure
        plot(err)
        title('Fehler RB über der Zeit');
               
        keyboard;
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Basis generation using different diffusivity constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 101
        
        params.coarse_factor = 4;
        
        model = convection_diffusion_fv_output_opt_model_rect(params);
        model.compute_derivative_info = 0;
        model.t = 0;
        
        if isfield(expar,'k')
            model.k = expar.k;
        else
            model.k = 0;
        end
        
        %RB related entries
        model.RB_stop_Nmax = 30;
        
        model_data = gen_model_data(model);
        
        detailed_data = gen_detailed_data(model, model_data);
        
        save(['DiffGreedybase_coarse',num2str(model.coarse_factor),'_Nmax_',num2str(model.RB_stop_Nmax),'_k_',num2str(model.k),'.mat'],...
            'model', 'detailed_data');
        
        keyboard;
        
    case 102
       
        fname = {'DiffGreedybase_coarse4_Nmax_30_k_0.05.mat','DiffGreedybase_coarse4_Nmax_30_k_0.01.mat',...
            'DiffGreedybase_coarse4_Nmax_30_k_0.005.mat','DiffGreedybase_coarse4_Nmax_30_k_0.001.mat'};
        
        nfiles = length(fname);
        
        k=cell(1,nfiles);
        tr_errs = zeros(nfiles,30);
        colr={'r','g','b','k'};
        
        figure
        for i=1:nfiles
            load(fname{i});
            k{i} = num2str(model.k);
            tr_errs(i,:) = detailed_data.RB_info.max_err_sequence;
            
            plot(detailed_data.RB_info.max_err_sequence,colr{i});
            hold on;
        end
        legend(k);
        keyboard;
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Test derivative calculation
% compare detailed derivative simulation with an fd-model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 103
        
        params.coarse_factor = 4;
        model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
        model.compute_derivative_info = 1;
        model.t = 0;
        
        model.k=0.25;
        
        mu=[0.1,0.3,0.9];
        model = set_mu(model, mu);
        
       model_data = gen_model_data(model);
       
       sim_data = detailed_simulation(model, model_data);
       U_der = sim_data.U_der;
       
       U_der_fd = lin_evol_opt_fd_derivative(model, model_data,mu);
       
       
       %Compare detailed derivative solutions
       U_diff = cell(1,length(U_der));
       for i=1:length(U_der)
           U_diff{i} = U_der{i} - U_der_fd{i};
           disp(['Maximum difference in ',num2str(i),'. derivative: ',num2str(max(max(U_diff{i})))]);
       end
       
       keyboard;
       
       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% verifiy derivative error estimators (PCA of sol+der trajectory)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 104
              
        
        params.coarse_factor = 16;
        %params.separate_bases = 1;
        %model = advection_fv_output_opt_model(params);
        model = convection_diffusion_fv_output_opt_model_rect_normed(params);
        model.compute_derivative_info = 1;
        model.t = 0;
        model.starting_time_step = 0;
        model.stopping_time_step = 120;
        
        model.k=0;
        
        mu=[0.1,0.3,0.9];
        model = set_mu(model, mu);
        
        %model.starting_time_step = 0;
        %model.stopping_time_step = 2;
                
        
       model_data = gen_model_data(model);
       
       sim_data = detailed_simulation(model, model_data);
       U_der = sim_data.U_der;
       U=sim_data.U;
       
       %construct mixed basis by PCA over solution snd derivative solutions
        detailed_data = model_data;
        detailed_data.RB = [];

        detailed_data.RB = model.rb_init_data_basis(model, detailed_data);


        %PCA
        traj = U;
        %for i = 1:length(U_der)
        %    traj = [traj,U_der{i}];
        %end
        PCAtr = PCA_fixspace(traj, detailed_data.RB,model_data.W);
        detailed_data.RB = [detailed_data.RB,PCAtr];
        

        %check orthogonality
        RB = detailed_data.RB;
        oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
        loop = 0;
        while (oc>1e-14)&&(loop<5)
            RB = orthonormalize_gram_schmidt(RB, model_data.W);
            oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
            loop = loop +1;
        end

        detailed_data.RB = RB;
        %detailed_data.RB_der{1} = RB;
        %detailed_data.RB_der{2} = RB;
        %detailed_data.RB_der{3} = RB;
        
        reduced_data = gen_reduced_data(model, detailed_data);

        rb_sim_data = rb_simulation(model, reduced_data);
        rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data);
        U_der_rb = rb_sim_data.U_der;
        
        %Auswertung
        nr_der = length(U_der);
        
        U_diff_der = cell(1,nr_der);
        err = zeros(nr_der,size(U,2));
        
        clr={'r','g','b','k','m'};
        
        disp(['Anzahl an Basisvektoren: ',num2str(size(detailed_data.RB,2))])
        
        figure(1);clf;
        figure(2);clf;
        figure(3);clf;
        for i=1:nr_der
            if ~isfield(model,'starting_time_step')
                U_diff_der{i} = U_der_rb(:,model.save_time_indices+1,i) - U_der{i};
                disp(['Maximum difference in ',num2str(i),'. rb derivative: ',num2str(max(max(U_diff_der{i})))]);
                err(i,:) = fv_l2_error(U_der_rb(:,model.save_time_indices+1,i), U_der{i}, model_data.grid,[]);
            else
                 U_diff_der{i} = U_der_rb(:,:,i) - U_der{i};
                disp(['Maximum difference in ',num2str(i),'. rb derivative: ',num2str(max(max(U_diff_der{i})))]);
                err(i,:) = fv_l2_error(U_der_rb(:,:,i), U_der{i}, model_data.grid,[]);
            end
            figure(1)
            plot(err(i,:),clr{i})
            hold on
            figure(2)
            plot(rb_sim_data.Delta_der(i,:),clr{i})
            hold on
            figure(3)
            plot(sqrt(rb_sim_data.res_u_der(i,:)),clr{i})
            hold on
        end
        figure(1)
        if ~isfield(model,'starting_time_step')
            plot(fv_l2_error(U,rb_sim_data.U(:,model.save_time_indices+1), model_data.grid,[]),'k')
        else
            plot(fv_l2_error(U,rb_sim_data.U, model_data.grid,[]),'k')
        end
        title('true error')
        legend('1','2','3','sol')
        figure(2)
        plot(rb_sim_data.Delta,'k')
        title('estimated error')
        legend('1','2','3','sol')
        figure(3)
        plot(sqrt(rb_sim_data.res_u),'k');
        title('residual')
        legend('1','2','3','sol')
        
       keyboard;
       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% verifiy derivative error estimators by separate basis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       
       
    case 105
        

        params.coarse_factor = 16;
        params.separate_bases = 1;
        %params.cone_const = 1;
        
        %model = advection_fv_output_opt_model(params);
        model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
        model.compute_derivative_info = 1;
        model.t = 0;
        
        
        
        
        model.k=0;
        
        mu=[0.1,0.3,0.9];
        model = set_mu(model, mu);
        
        %model.starting_time_step = 0;
        %model.stopping_time_step = 2;
                
        
       model_data = gen_model_data(model);
       
       sim_data = detailed_simulation(model, model_data);
       U_der = sim_data.U_der;
       U=sim_data.U;
       
       %construct mixed basis by PCA over solution snd derivative solutions
        detailed_data = model_data;
        detailed_data.RB = [];

        detailed_data.RB = model.rb_init_data_basis(model, detailed_data);


        %PCA_sol
        traj = U;
        PCAtr = model_orthonormalize_gram_schmidt(model, model_data,traj,1e-12);
        %PCAtr = PCA_fixspace(traj, detailed_data.RB,model_data.W,size(traj,2)-1);
        %PCAtr = model.PCA_fixspace(model, model_data, traj,detailed_data.RB);
        detailed_data.RB = [detailed_data.RB,PCAtr];
        
        
        %check orthogonality
        RB = detailed_data.RB;
        oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
        loop = 0;
        while (oc>1e-14)&&(loop<5)
            RB = orthonormalize_gram_schmidt(RB, model_data.W);
            oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
            loop = loop +1;
        end
        detailed_data.RB = RB;
        
        %PCA_der
        detailed_data.RB_der = cell(1,length(U_der));
        for i=1:length(U_der)
            detailed_data.RB_der{i} = model.rb_init_data_basis(model, detailed_data);
            traj = U_der{i};
            PCAtr = model_orthonormalize_gram_schmidt(model, model_data,traj,1e-12);
            %PCAtr = PCA_fixspace(traj, detailed_data.RB_der{i}, model_data.W,size(traj,2)-1);
            %PCAtr = model.PCA_fixspace(model, model_data, traj,detailed_data.RB,size(traj,2));
            detailed_data.RB_der{i} = [detailed_data.RB_der{i}, PCAtr];
            
            RB = detailed_data.RB_der{i};
            oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
            loop = 0;
            while (oc>1e-14)&&(loop<5)
                RB = orthonormalize_gram_schmidt(RB, model_data.W);
                oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
                loop = loop +1;
            end
            detailed_data.RB_der{i} = RB;
        end

                
        reduced_data = gen_reduced_data(model, detailed_data);

        rb_sim_data = rb_simulation(model, reduced_data);
        rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data);
        U_der_rb = rb_sim_data.U_der;
        
        %Auswertung
        nr_der = length(U_der);
        
        U_diff_der = cell(1,nr_der);
        err = zeros(nr_der,size(U,2));
        
        clr={'r','g','b','k'};
        
        disp(['Anzahl an Basisvektoren: ',num2str(size(detailed_data.RB,2))])
        
        figure(1);clf;
        figure(2);clf;
        figure(3);clf;
        for i=1:nr_der
            if ~isfield(model,'starting_time_step')
                U_diff_der{i} = U_der_rb(:,model.save_time_indices+1,i) - U_der{i};
                disp(['Maximum difference in ',num2str(i),'. rb derivative: ',num2str(max(max(U_diff_der{i})))]);
                err(i,:) = fv_l2_error(U_der_rb(:,model.save_time_indices+1,i), U_der{i}, model_data.grid,[]);
            else
                 U_diff_der{i} = U_der_rb(:,:,i) - U_der{i};
                disp(['Maximum difference in ',num2str(i),'. rb derivative: ',num2str(max(max(U_diff_der{i})))]);
                err(i,:) = fv_l2_error(U_der_rb(:,:,i), U_der{i}, model_data.grid,[]);
            end
            figure(1)
            plot(err(i,:),clr{i})
            hold on
            figure(2)
            plot(rb_sim_data.Delta_der(i,:),clr{i})
            hold on
            figure(3)
            plot(rb_sim_data.res_u_der(i,:),clr{i})
            hold on
        end
        figure(1)
        plot(fv_l2_error(U,rb_sim_data.U(:,model.save_time_indices+1), model_data.grid,[]),'k')
        title('true error')
        legend('1','2','3','sol')
        figure(2)
        plot(rb_sim_data.Delta,'k')
        title('estimated error')
        legend('1','2','3','sol')
        figure(3)
        plot(rb_sim_data.res_u,'k');
        title('residual')
        legend('1','2','3','sol')
        
       keyboard;
    
       
       
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Test derivative calculation with k as parameter
% compare detailed derivative simulation with an fd-model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 106
        
        params.coarse_factor = 8;
        params.k_const = 0;
        model = convection_diffusion_fv_output_opt_model_rect(params);
        model.compute_derivative_info = 1;
        model.t = 0;
        
        mu=[0.1,0.3,0.9,0.02];
        model = set_mu(model, mu);
        
       model_data = gen_model_data(model);
       
       sim_data = detailed_simulation(model, model_data);
       U_der = sim_data.U_der;
       
       U_der_fd = lin_evol_opt_fd_derivative(model, model_data,mu);
       
       
       %Compare detailed derivative solutions
       U_diff = cell(1,length(U_der));
       for i=1:length(U_der)
           U_diff{i} = U_der{i} - U_der_fd{i};
           disp(['Maximum difference in ',num2str(i),'. derivative: ',num2str(max(max(U_diff{i})))]);
       end
       
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% verifiy derivative error estimators (PCA of sol+der trajectory) with 4th
% parameter k
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 107
              
        
        params.coarse_factor = 8;
        params.k_const = 0;
        %params.separate_bases = 1;
        %model = advection_fv_output_opt_model(params);
        model = convection_diffusion_fv_output_opt_model_rect(params);
        model.compute_derivative_info = 1;
        model.t = 0;
        
       
        
        mu=[0.3,0.5,0.5,0.01];
        model = set_mu(model, mu);
        
        %model.starting_time_step = 0;
        %model.stopping_time_step = 2;
                
        
       model_data = gen_model_data(model);
       
       sim_data = detailed_simulation(model, model_data);
       U_der = sim_data.U_der;
       U=sim_data.U;
       
       %construct mixed basis by PCA over solution snd derivative solutions
        detailed_data = model_data;
        detailed_data.RB = [];

        detailed_data.RB = model.rb_init_data_basis(model, detailed_data);


        %PCA
        traj = U;
        for i = 1:length(U_der)
            traj = [traj,U_der{i}];
        end
        PCAtr = PCA_fixspace(traj, detailed_data.RB,model_data.W);
        detailed_data.RB = [detailed_data.RB,PCAtr];
        

        %check orthogonality
        RB = detailed_data.RB;
        oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
        loop = 0;
        while (oc>1e-14)&&(loop<5)
            RB = orthonormalize_gram_schmidt(RB, model_data.W);
            oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
            loop = loop +1;
        end

        detailed_data.RB = RB;
        %detailed_data.RB_der{1} = RB;
        %detailed_data.RB_der{2} = RB;
        %detailed_data.RB_der{3} = RB;
        
        reduced_data = gen_reduced_data(model, detailed_data);

        rb_sim_data = rb_simulation(model, reduced_data);
        rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data);
        U_der_rb = rb_sim_data.U_der;
        
        %Auswertung
        nr_der = length(U_der);
        
        U_diff_der = cell(1,nr_der);
        err = zeros(nr_der,size(U,2));
        
        clr={'r','g','b','k','m'};
        
        disp(['Anzahl an Basisvektoren: ',num2str(size(detailed_data.RB,2))])
        
        figure(1);
        figure(2);
        figure(3);
        for i=1:nr_der
            if ~isfield(model,'starting_time_step')
                U_diff_der{i} = U_der_rb(:,model.save_time_indices+1,i) - U_der{i};
                disp(['Maximum difference in ',num2str(i),'. rb derivative: ',num2str(max(max(U_diff_der{i})))]);
                err(i,:) = fv_l2_error(U_der_rb(:,model.save_time_indices+1,i), U_der{i}, model_data.grid,[]);
            else
                 U_diff_der{i} = U_der_rb(:,:,i) - U_der{i};
                disp(['Maximum difference in ',num2str(i),'. rb derivative: ',num2str(max(max(U_diff_der{i})))]);
                err(i,:) = fv_l2_error(U_der_rb(:,:,i), U_der{i}, model_data.grid,[]);
            end
             figure(1)
            plot(err(i,:),clr{i})
            hold on
            figure(2)
            plot(rb_sim_data.Delta_der(i,:),clr{i})
            hold on
            figure(3)
            plot(sqrt(rb_sim_data.res_u_der(i,:)),clr{i})
            hold on
        end
        figure(1)
        plot(fv_l2_error(U,rb_sim_data.U(:,model.save_time_indices+1), model_data.grid,[]),'k')
        title('true error')
        legend('1','2','3','4','sol')
        figure(2)
        plot(rb_sim_data.Delta,'k')
        title('estimated error')
        legend('1','2','3','4','sol')
        figure(3)
        plot(sqrt(rb_sim_data.res_u),'k');
        title('residual')
        legend('1','2','3','4','sol')
        
       keyboard;
       
       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% verifiy derivative error estimators by separate basis (with k as
% parameter)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       
       
    case 108
        

        params.coarse_factor = 8;
        params.k_const = 0;
        params.separate_bases = 1;
        %model = advection_fv_output_opt_model(params);
        model = convection_diffusion_fv_output_opt_model_rect(params);
        model.compute_derivative_info = 1;
        model.t = 0;
        
        model.k=0;
        
        mu=[0.3,0.5,0.5,0.5];
        model = set_mu(model, mu);
        
        %model.starting_time_step = 0;
        %model.stopping_time_step = 2;
                
        
       model_data = gen_model_data(model);
       
       sim_data = detailed_simulation(model, model_data);
       U_der = sim_data.U_der;
       U=sim_data.U;
       
       %construct mixed basis by PCA over solution snd derivative solutions
        detailed_data = model_data;
        detailed_data.RB = [];

        detailed_data.RB = model.rb_init_data_basis(model, detailed_data);


        %PCA_sol
        traj = U;
        PCAtr = PCA_fixspace(traj, detailed_data.RB,model_data.W);
        detailed_data.RB = [detailed_data.RB,PCAtr];
        
        
        %check orthogonality
        RB = detailed_data.RB;
        oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
        loop = 0;
        while (oc>1e-14)&&(loop<5)
            RB = orthonormalize_gram_schmidt(RB, model_data.W);
            oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
            loop = loop +1;
        end
        detailed_data.RB = RB;
        
        %PCA_der
        detailed_data.RB_der = cell(1,length(U_der));
        for i=1:length(U_der)
            detailed_data.RB_der{i} = model.rb_init_data_basis(model, detailed_data);
            traj = U_der{i};
            PCAtr = PCA_fixspace(traj, detailed_data.RB_der{i}, model_data.W);
            detailed_data.RB_der{i} = [detailed_data.RB_der{i}, PCAtr];
            
            RB = detailed_data.RB_der{i};
            oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
            loop = 0;
            while (oc>1e-14)&&(loop<5)
                RB = orthonormalize_gram_schmidt(RB, model_data.W);
                oc = max(max(RB'*model_data.W*RB-eye(size(RB,2), size(RB,2))))
                loop = loop +1;
            end
            detailed_data.RB_der{i} = RB;
        end

                
        reduced_data = gen_reduced_data(model, detailed_data);

        rb_sim_data = rb_simulation(model, reduced_data);
        rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data);
        U_der_rb = rb_sim_data.U_der;
        
        %Auswertung
        nr_der = length(U_der);
        
        U_diff_der = cell(1,nr_der);
        err = zeros(nr_der,size(U,2));
        
        clr={'r','g','b','m','k'};
        
        disp(['Anzahl an Basisvektoren: ',num2str(size(detailed_data.RB,2))])
        
        figure(1);
        figure(2);
        figure(3);
        for i=1:nr_der
            if ~isfield(model,'starting_time_step')
                U_diff_der{i} = U_der_rb(:,model.save_time_indices+1,i) - U_der{i};
                disp(['Maximum difference in ',num2str(i),'. rb derivative: ',num2str(max(max(U_diff_der{i})))]);
                err(i,:) = fv_l2_error(U_der_rb(:,model.save_time_indices+1,i), U_der{i}, model_data.grid,[]);
            else
                 U_diff_der{i} = U_der_rb(:,:,i) - U_der{i};
                disp(['Maximum difference in ',num2str(i),'. rb derivative: ',num2str(max(max(U_diff_der{i})))]);
                err(i,:) = fv_l2_error(U_der_rb(:,:,i), U_der{i}, model_data.grid,[]);
            end
            figure(1)
            plot(err(i,:),clr{i})
            hold on
            figure(2)
            plot(rb_sim_data.Delta_der(i,:),clr{i})
            hold on
            figure(3)
            plot(rb_sim_data.res_u_der(i,:),clr{i})
            hold on
        end
        figure(1)
        plot(fv_l2_error(U,rb_sim_data.U(:,model.save_time_indices+1), model_data.grid,[]),'k')
        title('true error')
        legend('1','2','3','4','sol')
        figure(2)
        plot(rb_sim_data.Delta,'k')
        title('estimated error')
        legend('1','2','3','4','sol')
        figure(3)
        plot(rb_sim_data.res_u,'k');
        title('residual')
        legend('1','2','3','4','sol')
        
       keyboard;
    
       
       
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Konstantenbestimmung mit diffusion
   %
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
       
   
    case 109
        
        params.coarse_factor = 16;
        params.k_const = 0;
        
        M_test = [0,0,0,0;...
                    1,0,0,0;...
                    0,1,0,0;...
                    0,0,1,0;...
                    0,0,0,1;...
                    1,1,1,1;...
                    1,1,1,1];
        
        
        model = convection_diffusion_fv_output_opt_model_rect_normed(params);
        
        model.t=0;
        model.k=0.05;
        model.compute_derivative_info = 1;
        model.decomp_mode = 0;
        
        model_data = gen_model_data(model);
        
        [L_I,L_E,b] = model.operators_ptr(model, model_data);
        
        %check if linear cimobination and decomp mode is the same
        model.decomp_mode = 1;
        [L_I_comp,L_E_comp,b_comp] = model.operators_ptr(model, model_data);
        
        model.decomp_mode = 2;
        [sL_I,sL_E,sb] = model.operators_ptr(model, model_data);
        
        L_I2 = lincomb_sequence(L_I_comp,sL_I);
        L_E2 = lincomb_sequence(L_E_comp,sL_E);
        b2 = lincomb_sequence(b_comp,sb);
        
        diffI=L_I - L_I2;
        diffE = L_E - L_E2;
        diffb = b-b2;
        disp('0 0 0?')
        max(max(diffI))
        max(max(diffE))
        max(max(diffb))
        
        C_LEmax = 0;
        C_LImax = 0;
        
        disp(['Norm der L_I_comp bei coarse_factor',num2str(model.coarse_factor),...
             'ist',num2str(norm(full(L_I_comp{1})))])
        
        
        for i=1:size(M_test,1)
            mu=M_test(i,:);
            model = set_mu(model,mu);
            model.decomp_mode = 2;
            [sL_I,sL_E,sb] = model.operators_ptr(model, model_data);
            L_I2 = lincomb_sequence(L_I_comp,sL_I);
            L_E2 = lincomb_sequence(L_E_comp,sL_E);
            b2 = lincomb_sequence(b_comp,sb);
            
        %%%%
        % Konstante ||Id+dt*L_E||
                
            C_LE = norm(eye(size(L_E,1))+model.dt*full(L_E));
            if C_LE >C_LEmax
                C_LEmax = C_LE;
            end
            
        
        %%%%
        % Konstante ||(Id-dt*L_I)^(-1)||
        
        
             C_LI = norm(inv(eye(size(L_I2,1))-model.dt*full(L_I2)));
             if C_LI > C_LImax
                 C_LImax = C_LI;
             end
        
            
        end
        
        disp(['C_E Konstante ist: ',num2str(C_LEmax)]);
        disp(['C_I Konstante ist: ',num2str(C_LImax)]);
        %Beide Konstanten sind <=1
        
        
        %%%%
        % Derivative Konstanten C_dE und C_dI
        
        
        
        for p=1:4
            C_LEdmax=0;
           C_LId = zeros(1,size(M_test,1));
        for i=1:size(M_test,1)
           [sd_LI, sd_LE, sd_b] = model.rb_derivative_operators_coefficients_ptr(model);
           
            L_Id = lincomb_sequence(L_I_comp,sd_LI(p));
            L_Ed = lincomb_sequence(L_E_comp,sd_LE(p,:));
            bd = lincomb_sequence(b_comp,sd_b(p,:));
                
            C_LId(i) = norm(full(L_Id));
            C_LEd = norm(full(L_Ed));
            
            %if C_LId > C_LIdmax
            %    C_LIdmax = C_LId;
            %end
            if C_LEd > C_LEdmax
                C_LEdmax = C_LEd;
            end
        end
        C_LIdmax = max(C_LId);
        disp(['Ableitungsoperatorkonstanten für ',num2str(p),'.Abl. und coarse_factor ',num2str(model.coarse_factor)]);
        disp(['C_Lid : ',num2str(C_LIdmax)]);
        disp(['C_Led : ', num2str(C_LEdmax)]);
        disp(['Ableitungskontanten für ',num2str(p),'. Abl. inkl. dt:']);
        disp(['C_Lid : ',num2str(model.dt*C_LIdmax)]);
        disp(['C_Led : ', num2str(model.dt*C_LEdmax)]);
        end
        disp('Änderung der C_LId Konstante mit Parameter k:')
        C_LId
        
        
        %Operatornorm des Rieszrepresentanten des outputs ||J||
        % |J(u)|<=||J||||u||
        %
        % Es gilt ||v|| = sup_{||w||=1;w\in\R^H} v*K*K^{-0.5}*w;%
        % z = v*K*K^{-0.5}
        %--> ||v|| = norm(z)
        disp('')
        disp('Bestimmung von ||J||');
        disp('')
        coarse_factor = 1;
        model = convection_diffusion_fv_output_opt_model_rect_normed(params);
        model.decomp_mode = 1;
        model.t=0;
        model_data = gen_model_data(model);
        
        v=model.operators_output(model, model_data);
        K=model_data.W;
        z=v{1}'*sqrt(inv(K));
        disp(['||v||= ',num2str(norm(z))]);
        
        keyboard;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% lipschitz constants for functional gradient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
    case 1092 %obsolet!!!!!!!!!!!!!!!!!!!!!!!!
        
        %Lipschitz Konstanten für verschiedene Fälle (zu optimierende
        %Variablen, coarse Faktoren) in arrays gespeichert und als File
        %geschrieben.
        % Das Model sollte dieses File laden und die für sich relevante
        % Lipschitzkonstante auslesen.
        % Zur verbesserung der Genauigkeit können auch loka gültige
        % Lipschitzkonstanten gespeichert werden.
        
        %fix model_type here
        
        %model_type = 'normal';
        %model_type = 'vx_const';
        model_type = 'cone_const';
        %time_type = 'time_const';
        time_type = 'time_var';
        
        if (nargin<3)||(~isfield(expar,'coarse_factor'))
            params.coarse_factor =2;
        end    
        
        %fix here grid density for localized lipschitz constant
        grid_dens=10;
        %grid density of subgrid on each lipschitzgrid element
        sub_grid_dens=1;
        
        if strcmp(time_type,'time_var')
            params.time_const=0;
        else
            params.time_const=1;
        end
        
        switch model_type
            case 'vx_const'
                params.vx_const=1;
                params_to_optimize=[1,0,1];
                grid_params.range={[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1];
            case 'cone_const'
                params.cone_const=1;
                params_to_optimize=[0,1,1];
                grid_params.range={[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1];
            otherwise
                params.normal_model=1;
                params_to_optimize=[1,1,1];
                grid_params.range={[0,1],[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1,1];
        end
        
        model = convection_diffusion_fv_output_opt_model_rect_normed(params);
        model.t=0;
        model.k=0.25;
        model.compute_derivative_info = 1;
        model.decomp_mode = 0;
        model.optimization.params_to_optmize= params_to_optimize;
        model_data = gen_model_data(model);
        
        %generate lipschitz parameter grid
        lip_grid = cubegrid(grid_params);
        
        lipschitz_const = zeros(1,lip_grid.nelements);
        
        model.decomp_mode=0;
        v=model.operators_output(model, model_data);
        
        for i=1:lip_grid.nelements
            %for every element i in the grid we have to calculate the
            %lipschitz constant
            
            %generating subgrid
            sub_par.range=get_ranges_of_element(lip_grid,i);
            sub_par.numintervals = sub_grid_dens*ones(1,lip_grid.dimension);
            subgrid=cubegrid(sub_par);
            Mtest = subgrid.vertex;
            
            
            Cl = zeros(1,size(Mtest,1));
            for j=1:size(Mtest,1)
                 model = set_mu(model, Mtest(j,:));
                 %2nd derivative
                 U_der2 = lin_evol_opt_fd_second_derivative(model, model_data, Mtest(j,:));
                 functional = zeros(length(U_der2),1);
                 %functional of 2nd der.
                 for m=1:length(U_der2)
                     functional(m)=v'*U_der2{m}(:,end);
                 end
                 Cl(j) = norm(functional);
                 
            end
            
            lipschitz_const(i) = max(Cl);
                               
        end
               
        lipschitz_data=[];
        lipschitz_data.lip_grid = lip_grid;
        lipschitz_data.lipschitz_const=lipschitz_const;
        
        save(fullfile(fdir,['lipschitz_data_diff_new',num2str(model.k),'_',model_type,'_',time_type,'_coarse_',num2str(params.coarse_factor)]),'lipschitz_data')
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculating Lipschitz constants of the functional gradient for every mu-direction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 1093
        
        %Lipschitz Konstanten für verschiedene Fälle (zu optimierende
        %Variablen, coarse Faktoren) in arrays gespeichert und als File
        %geschrieben.
        % Das Model sollte dieses File laden und die für sich relevante
        % Lipschitzkonstante auslesen.
        % Zur verbesserung der Genauigkeit können auch loka gültige
        % Lipschitzkonstanten gespeichert werden.
        
        %fix model_type here
        
        %model_type = 'normal';
        %model_type = 'vx_const';
        model_type = 'cone_const';
        %time_type = 'time_const';
        time_type = 'time_var';
        
        %params.functional = 'box_mean';
        params.functional = 'box_mean_time';
        
        if (nargin>1)&&(isfield(expar,'coarse_factor'))
            params.coarse_factor = expar.coarse_factor;
        else
            params.coarse_factor =2;
        end
        
        %fix here grid density for localized lipschitz constant
        grid_dens=10;
        %grid density of subgrid on each lipschitzgrid element
        sub_grid_dens=1;
        
        if strcmp(time_type,'time_var')
            params.time_const=0;
        else
            params.time_const=1;
        end
        
        switch model_type
            case 'vx_const'
                params.vx_const=1;
                params_to_optimize=[1,0,1];
                grid_params.range={[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1];
            case 'cone_const'
                params.cone_const=1;
                params_to_optimize=[0,1,1];
                grid_params.range={[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1];
            otherwise
                params.normal_model=1;
                params_to_optimize=[1,1,1];
                grid_params.range={[0,1],[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1,1];
        end
        
        model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
        model.t=0;
        model.k=0.25;
        model.compute_derivative_info = 1;
        model.decomp_mode = 0;
        %model.optimization.params_to_optmize= params_to_optimize;
        model_data = gen_model_data(model);
        
        %generate lipschitz parameter grid
        lip_grid = cubegrid(grid_params);
        
        lipschitz_const = zeros(1,lip_grid.nelements);
        
        model.decomp_mode=0;
        v=model.operators_output(model, model_data);
        
        for i=1:lip_grid.nelements
            %for every element i in the grid we have to calculate the
            %lipschitz constant
            
            %generating subgrid
            sub_par.range=get_ranges_of_element(lip_grid,i);
            sub_par.numintervals = sub_grid_dens*ones(1,lip_grid.dimension);
            subgrid=cubegrid(sub_par);
            Mtest = subgrid.vertex;
            
            
            %Cl = zeros(sum(params_to_optimize),sum(params_to_optimize));
            Hes_norm =0;
            for j=1:size(Mtest,1)
                model = set_mu(model, Mtest(j,:));
                %Hessian
                Hes = filecache_function(@lin_evol_opt_fd_Hessian,model, model_data, Mtest(j,:));
                
                
                %functional of every Hessian entry.... find biggest Hessian norm
                help_ar = zeros(lip_grid.dimension,lip_grid.dimension);
                if strcmp(params.functional, 'box_mean_time')
                    for k=1:size(Hes,1)
                        for l=1:size(Hes,2)
                            model.compute_derivative_info = 0;
                            sim_data_dummy.U = Hes{k,l};
                            sim_data_dummy = lin_evol_opt_output_time_integrated(model, sim_data_dummy,model_data);
                            help_ar(k,l) = sim_data_dummy.y(end);
                        end
                    end
                    
                else
                    for m=1:size(Hes,1)
                        for n=1:size(Hes,2)
                            help_ar(m,n)=v'*Hes{1}(:,end);
                        end
                    end
                end
                buf = norm(help_ar);%norm of hessian matrix for this test parameter
                
                if buf>Hes_norm
                    Hes_norm = buf;
                end
                
            end
            
            lipschitz_const(i) = Hes_norm;
            
        end
        
        lipschitz_data=[];
        lipschitz_data.lip_grid = lip_grid;
        lipschitz_data.lipschitz_const=lipschitz_const;
        
        if strcmp(model.name_output_functional,'box_mean')
            save(fullfile(fdir,['lipschitz_data_diff',num2str(model.k),'_Hessian_',model_type,'_',time_type,'_coarse_',num2str(params.coarse_factor),'.mat']),'lipschitz_data')
        elseif strcmp(model.name_output_functional,'box_time_int')
            save(fullfile(fdir,['lipschitz_data_diff',num2str(model.k),'_timeInt_',model_type,'_',time_type,'_coarse_',num2str(params.coarse_factor),'.mat']),'lipschitz_data')
        else
            error('no functional defined')
        end
        
        
        %keyboard;
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% separate bases basis generation using only small parameter space
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

    case 110
        
        params.coarse_factor = 2;
        params.cone_const = 1;
        params.separate_bases = 1;
        
        model = convection_diffusion_fv_output_opt_model_rect(params);
        model.RB_solution_Nmax = 200;
        model.RB_solution_epsilon = 1e-6;
        model.RB_extension_algorithm = @RB_extension_PCA_fixspace_without_filecache;
        
        
        model.k=0.25;
        model.compute_derivative_info=1;
        model.RB_numintervals = [3,3];
        model.RB_stop_Nmax = 100;
        model.RB_stop_epsilon=5e-4;
        model.mu_ranges={[0.8,0.9],[0.6,0.7]};
        %model.starting_time_step = 0;
        %model.stopping_time_step = 50;
        
        model_data = gen_model_data(model);
        detailed_data = gen_detailed_data(model, model_data);
        
        keyboard;
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% separate bases generation
% 2 parameters + diffusion and t-partition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        case 111
        
        params.coarse_factor = 16;
        params.cone_const = 1;
        params.separate_bases = 1;
        
        model = convection_diffusion_fv_output_opt_model_rect_normed(params);
        model.RB_extension_algorithm = @RB_extension_PCA_fixspace_without_filecache;
        model.verbose = 8;
        model.RB_solution_Nmax = 12;
        model.RB_solution_epsilon = 5e-8;
        model.RB_stop_Nmax = 8;
        model.RB_stop_epsilon = 5e-7;
        %t-part params
        nr_tparts = 2;
        t_part_range = cell(1,nr_tparts);
        for i=1:nr_tparts
            t_part_range{i} = [(i-1)*floor((model.nt)/nr_tparts),(i)*floor((model.nt)/nr_tparts)];
        end
        t_part_range{nr_tparts}(2) = model.nt;
        t_params.t_part_range=t_part_range;
        t_params.t_part_rb_generation_mode = 'fixed';
                    
        model.k=0.25;
        model.compute_derivative_info=1;
        model.RB_numintervals = [3,3];
        
        model.mu_ranges={[0.4,0.5],[0.4,0.5]};
                
        model = t_part_model(model,t_params);
        
        
        model_data = gen_model_data(model);
        detailed_data = gen_detailed_data(model, model_data);
        
        filename=['newDiffGreedyBase_tpart',num2str(nr_tparts),'_2p_coarse',num2str(model.coarse_factor),...
            '_Nmax',num2str(model.RB_stop_Nmax),'_k_',num2str(model.k),'_small0405.mat'];
        save(fullfile(fdir,filename),'model','detailed_data','model_data');
        
        %keyboard;
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% separate bases generation (p+t) for selected partitions
% 2 parameters + diffusion + ppartition and t-partition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        case 112
        
        filecache_clear;
            
        params.coarse_factor = 2; 
        params.cone_const = 1;
        params.separate_bases = 1;
        
        p_params.p_part_numintervals = [10,10];
        p_params.p_part_generation_mode = 'uniform';
        
        model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
        model.verbose = 3;
        %t-part params
        %!!!!!!
        nr_tparts = 5;
        %!!!!!!!!!!!!!
        t_part_range = cell(1,nr_tparts);
        for i=1:nr_tparts
            t_part_range{i} = [(i-1)*floor((model.nt)/nr_tparts),(i)*floor((model.nt)/nr_tparts)];
        end
        t_part_range{nr_tparts}(2) = model.nt;
        t_params.t_part_range=t_part_range;
        t_params.t_part_rb_generation_mode = 'fixed';
                    
        model.k=0.25;
        model.compute_derivative_info=1;
        model.RB_numintervals = [3,3];
        model.RB_solution_Nmax = 150; 
        model.RB_solution_epsilon = 1e-7;
        model.RB_stop_Nmax = 150;
        model.RB_stop_epsilon = 5e-7;
        model.mu_ranges={[0,1],[0,1]};
        model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
        
        model = t_part_model(model,t_params);
        
        model = p_part_model(model,p_params);
        %!!!!!!!!!
        model.gen_detailed_data = @p_part_gen_detailed_data_parallel_selected2;
        %!!!!!!!!!!!!
        
        if isfield(expar,'selected_partitions')
            model.selected_partitions = expar.selected_partitions;
        else
            model.selected_partitions = 1:100;%[23,34,45,58,59,67,68,69];
            disp('no partitions selected')
        end
        model.verbose = 3;
        
        model_data = gen_model_data(model);
        detailed_data = gen_detailed_data(model, model_data);
        
        fdir = fullfile(fdir,'coarse2');
        filename=['DiffGreedyBase_ptpart',num2str(nr_tparts),'_2p_coarse',num2str(model.base_model.coarse_factor),...
            '_Nmax',num2str(model.base_model.RB_stop_Nmax),'_nrtparts_',num2str(nr_tparts),'_k_',num2str(model.base_model.k),...
            '_selected',num2str(model.selected_partitions),'.mat'];
        fn = fullfile(fdir,filename);
        save(fn,'-v7.3','model','detailed_data', 'model_data');
        
        %keyboard;        
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% calculate average dimension of tp-basis (all equally weighted)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 113
        
        load(fullfile(fdir,'coarse2','reduced_data_complete.mat'));
        
        %load('Copy_of_DiffGreedyBase_ptpart10_2p_coarse4_Nmax30_k_0.25_selected.mat')
        
        part_used = 1:100;%[34,45,67,68,68,69,69,69];
        
        sumN=[0,0,0];
        
        nr_tparts = length(reduced_data.parts_reduced_data{part_used(1)}.t_part_reduced_data);
        for i=1:length(part_used)
            rd=reduced_data.parts_reduced_data{part_used(i)};
            for j=1:length(rd.t_part_reduced_data)
                sumN(1) = sumN(1)+rd.t_part_reduced_data{j}.N;
                sumN(2) = sumN(2)+length(rd.t_part_reduced_data{j}.c0{1}{1});
                sumN(3) = sumN(3)+length(rd.t_part_reduced_data{j}.c0{2}{1});
            end
        end
        
        N=size(part_used,2)*nr_tparts;
        sumN = sumN./N;
        
        disp(['Durchschnittliche Basenzahl: ',num2str(sumN)]);
        keyboard;
        
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generiere eine p_part Trickbasis 
% PCA von Trajektorien zu vorgegebenen Parametern werden berechnet
% anschließend werden diese PCAs als reduzierte Basen in ein p_part_Gitter
% eingeordnet.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
    case 114
        
        params.coarse_factor = 2;
        params.cone_const = 1;
        params.separate_bases = 1;
        
               
        model = convection_diffusion_fv_output_opt_model_rect_normed(params);
        model.compute_derivative_info = 1;
        model_data = gen_model_data(model);
        %what are the given parameters?
        if isfield(expar,'M_train')
            M_train = expar.M_train;
        else
            M_train = get_mu(model)';
        end
        
        %how does the p_part grid look like?
        if isfield(expar,'p_part_numintervals')
            p_params.p_part_numintervals = expar.p_part_numintervals;
        else
            p_params.p_part_numintervals = [10,10];
        end
        p_params.p_part_generation_mode = 'uniform';
        
        
        %preparations t-part_model
        nr_tparts = 10;
        
        t_part_range = cell(1,nr_tparts);
        nt_ind = length(model.save_time_indices);
        for i=1:nr_tparts
            t_part_range{i} = [(i-1)*floor((nt_ind)/nr_tparts),(i)*floor((nt_ind)/nr_tparts)];
        end
        t_part_range{1}(1) = 0;
        t_part_range{nr_tparts}(2) = nt_ind-1;
        t_params.t_part_range=t_part_range;
        t_params.t_part_rb_generation_mode = 'fixed';
        model = t_part_model(model,t_params);
        
        %concluding model
        model = p_part_model(model, p_params);
        
        %generate p-grid
        g_params.range = model.p_part_ranges;
        g_params.numintervals = model.p_part_numintervals;
        pgrid = cubegrid(g_params);
        
        %generate empty detailed_data structure
        detailed_data = [];
        detailed_data.pgrid = pgrid;
        detailed_data.parts_detailed_data = cell(get(pgrid,'nelements'),1);
        
        %initialize t-part-detailed_data
        t_part_detailed_data = cell(nr_tparts,1);
        for i=1:nr_tparts
            t_part_detailed_data{i}.W = model_data.W;
            t_part_detailed_data{i}.grid = model_data.grid;
            t_part_detailed_data{i}.t_part_range = t_part_range;
        end
        
        %init RB basis
        RB_init = model.base_model.rb_init_data_basis(model.base_model, model_data);
        
        for i=1:size(M_train,1)
           %loop over all elements in M_train
           model = set_mu(model, M_train(i,:));
           %detailed simulation for this parameter and PCA generation
           sim_data = lin_evol_opt_detailed_der_simulation_t_part(model.base_model, model_data);
           %RB_ext = PCA_fixspace(sim_data.U,RB_init,model_data.W);
                     
           %if Basis in this parameter partition exists already-->extend
           %else create new 
           %only save solutions; orthonormalization and PCA later
           p=coord2leaf_element(pgrid,M_train(i,:));
           if isempty(detailed_data.parts_detailed_data{p})
               for j=1:nr_tparts
                    t_part_detailed_data{j}.RB = sim_data.U(:,(t_part_range{j}(1):t_part_range{j}(2))+1);
                    for d=1:length(sim_data.U_der)
                        t_part_detailed_data{j}.RB_der{d} = sim_data.U_der{d}(:,(t_part_range{j}(1):t_part_range{j}(2))+1);
                    end
               end
               %attach initial conditions in first t_part
               t_part_detailed_data{1}.RB = [RB_init,t_part_detailed_data{1}.RB];
               for d=1:length(t_part_detailed_data{1}.RB_der)
                    t_part_detailed_data{1}.RB_der{d} = [RB_init,t_part_detailed_data{1}.RB_der{d}];
               end
               %detailed_data.parts_detailed_data{p}.RB = model_orthonormalize_qr(model.base_model,...
               %    model_data,[RB_init,RB_ext]);
               detailed_data.parts_detailed_data{p}.t_part_detailed_data = t_part_detailed_data;
               detailed_Data.parts_detailed_data{p}.t_part_range = t_part_range;
               detailed_data.parts_detailed_data{p}.RB_info.M_train =  M_train(i,:);
               detailed_data.parts_detailed_data{p}.grid = model_data.grid;
               detailed_data.parts_detailed_data{p}.W = model_data.W;
           else
               %RB = model_orthonormalize_qr(model.base_model, model_data,...
               %    [detailed_data.parts_detailed_data{p}.RB, RB_ext]);
               %detailed_data.parts_detailed_data{p}.RB = RB;
               for j=1:nr_tparts
                    t_part_detailed_data{j}.RB = [detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB,...
                        sim_data.U(:,(t_part_range{j}(1):t_part_range{j}(2))+1)];
                    for d=1:length(sim_data.U_der)
                        t_part_detailed_data{j}.RB_der{d} = [detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB_der{d},...
                            sim_data.U_der{d}(:,(t_part_range{j}(1):t_part_range{j}(2))+1)];
                    end
               end
               detailed_data.parts_detailed_data{p}.t_part_detailed_data = t_part_detailed_data;
               detailed_data.parts_detailed_data{p}.RB_info.M_train = ...
                            [detailed_data.parts_detailed_data{p}.RB_info.M_train;M_train(i,:)];
           end
           
        end
        
        %make a PCA in all p-parts and t-parts over the trajectories
        for p=1:length(detailed_data.parts_detailed_data)
            if ~isempty(detailed_data.parts_detailed_data{p})
                for j=1:nr_tparts
                   PCAtr = PCA_fixspace(detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB(:,2:end),...
                       detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB(:,1),model_data.W);
                   detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB = ...
                       [detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB(:,1),PCAtr];
                   detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB = model_orthonormalize_qr(model.base_model, model_data,...
                       detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB);
                   for d=1:length(detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB_der)
                        PCAtr = PCA_fixspace(detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB_der{d}(:,2:end),...
                            detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB_der{d}(:,1),model_data.W);
                        detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB_der{d} = ...
                            [detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB_der{d}(:,1),PCAtr];
                        detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB_der{d} = ...
                            model_orthonormalize_qr(model.base_model, model_data,detailed_data.parts_detailed_data{p}.t_part_detailed_data{j}.RB_der{d});
                   end
                end
            end
        end
        
        
        
        reduced_data = gen_reduced_data(model, detailed_data);
        
        keyboard;
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ZUsammenfügen der einzelnen Partitionsn, die in case 112 erzeugt wurden
% !achtung: in dem Verzeichnis dürfen nur die benötigten dateien liegen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
        
    case 115
        
        %fdir2 = fullfile(fdir,'coarse4');
        %cd /home/dihlmann/matlab/rbmatlab/rbasis/lin_evol_opt/convdiff/coarse2
        cd /home/dihlmann/matlab/rbmatlab-core/rbasis/lin_evol_opt/convdiff/coarse2
        
        s = what;
        f = s.mat;
        
        load(f{1});
        dd = detailed_data;
        
        disp('start appending bases')
        
        for i=2:length(f);
            disp(i)
            load(f{i});
            for p=1:length(detailed_data.parts_detailed_data)
                if ~isempty(detailed_data.parts_detailed_data{p})
                    dd.parts_detailed_data{p} = detailed_data.parts_detailed_data{p};
                end
            end
        end
                
        
        detailed_data = dd;
        disp('saving the complete basis')
        save -v7.3 complete_TP_basis_coarse2.mat model detailed_data    
        
        disp('Generating the reduced data');
        reduced_data = gen_reduced_data(model, detailed_data);
        disp('saving reduced_data')
        save -v7.3 reduced_data_complete.mat model reduced_data
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% generiere Landschaftsplots für Funktional, Gradient und Hessian
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 116
        
        if (nargin>1)&&(isfield(expar,'coarse_factor'))
            params.coarse_factor = expar.coarse_factor;
        else
            params.coarse_factor = 16;
        end
        
        params.cone_const =1;
        %params.functional = 'box_mean';
        %params.functional = 'box_mean_time';
        params.functional = 'box_mean_time2';
        
        %model preparation
        model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
        model.t=0;
        model.k=0.25;
        model.compute_derivative_info = 1;
        model.decomp_mode = 0;
        model_data = gen_model_data(model);
        
        %set landscape grid options (density...)
        if (nargin>1)&&(isfield(expar,'num_pts'))
            num_pts = expar.num_pts;
        else
            num_pts = 6;
        end
        p_range = [0,1];
        delt = (p_range(2)-p_range(1))/(num_pts-1);
        mu_test = p_range(1):delt:p_range(2);
        
        
        %Ouput functional
        v=model.operators_output(model, model_data);
        
        %landscape values - initialization
        J = zeros(num_pts,num_pts); %functional values
        grad1 = zeros(num_pts, num_pts); %gradient first entry
        grad2 = zeros(num_pts, num_pts); %gradient second entry
        H11 = zeros(num_pts, num_pts); %Hessian 11entry
        H12 = zeros(num_pts, num_pts); %Hessian 12 entry
        H22 = zeros(num_pts, num_pts); % Hessian 22 entry
        
        %loop over parameter grid
        for i=1:num_pts
            for j=1:num_pts
                H= zeros(2,2);
               mu = [mu_test(i),mu_test(j)];
               model = set_mu(model, mu);
               sim_data = detailed_simulation(model, model_data);
               H_sol = lin_evol_opt_fd_Hessian(model, model_data,mu);
                for k=1:2
                    for l=1:2
                        H(k,l)=v'*H_sol{k,l}(:,end);
                    end
                end

                %save results
                J(i,j) = sim_data.y(end);
                grad1(i,j) = sim_data.y_der(1,end);
                grad2(i,j) = sim_data.y_der(2,end);
                H11(i,j) = H(1,1);
                H12(i,j) = H(1,2);
                H22(i,j) = H(2,2);
            end
        end
        
        figure;surf(mu_test,mu_test,J);title(['J_',num2str(params.coarse_factor)]);figure;
        surf(mu_test,mu_test,grad1);title(['grad1',num2str(params.coarse_factor)]);figure;
        surf(mu_test,mu_test,grad2);title(['grad2',num2str(params.coarse_factor)]);figure;
        surf(mu_test,mu_test,H11);title(['H11',num2str(params.coarse_factor)]);figure;
        surf(mu_test,mu_test,H12);title(['H12',num2str(params.coarse_factor)]);figure;
        surf(mu_test,mu_test,H22);title(['H22',num2str(params.coarse_factor)])
        
        keyboard;
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% test lipschitz constant grid-zuweisung... für case 1093
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 117
        
        %fix model_type here
        
        %model_type = 'normal';
        %model_type = 'vx_const';
        model_type = 'cone_const';
        %time_type = 'time_const';
        time_type = 'time_var';
        
        if (nargin>1)&&(isfield(expar,'coarse_factor'))
            params.coarse_factor = expar.coarse_factor;
        else
            params.coarse_factor =16;
        end    
        
        %fix here grid density for localized lipschitz constant
        grid_dens=10;
        
        
        
        if strcmp(time_type,'time_var')
            params.time_const=0;
        else
            params.time_const=1;
        end
        
        switch model_type
            case 'vx_const'
                params.vx_const=1;
                params_to_optimize=[1,0,1];
                grid_params.range={[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1];
            case 'cone_const'
                params.cone_const=1;
                params_to_optimize=[0,1,1];
                grid_params.range={[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1];
            otherwise
                params.normal_model=1;
                params_to_optimize=[1,1,1];
                grid_params.range={[0,1],[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1,1];
        end
        
        model = convection_diffusion_fv_output_opt_model_rect_normed(params);
        model.t=0;
        model.k=0.25;
        model.compute_derivative_info = 1;
        model.decomp_mode = 0;
        model.optimization.params_to_optmize= params_to_optimize;
        model_data = gen_model_data(model);
        
        %generate lipschitz parameter grid
        %lip_grid = cubegrid(grid_params);
        
        
        delt = 1/grid_dens;
        lip_vec = 0:delt:1;
        numpts = length(lip_vec);
        lipschitz_const = zeros(numpts,numpts);
        
        model.decomp_mode=0;
        v=model.operators_output(model, model_data);
        
        for i=1:numpts
            for j=1:numpts
            %for every element i in the grid we have to calculate the
            %lipschitz constant
            
            mu=[lip_vec(i),lip_vec(j)];
            
                model = set_mu(model, mu);
                %Hessian
                Hes = lin_evol_opt_fd_Hessian(model, model_data, mu);
                
                
                %functional of every Hessian entry.... find biggest Hessian norm
                help_ar = zeros(2,2);
                 for m=1:size(Hes,1)
                     for n=1:size(Hes,2)
                         help_ar(m,n)=v'*Hes{1}(:,end);
                     end
                 end
                 Hes_norm = norm(help_ar);%norm of hessian matrix for this test parameter
                 lipschitz_const(i,j) = Hes_norm;
                 
            end
            
            
                               
        end
               
        keyboard;
        
        
    case 118
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % for the new non-iterative error estimator in optimization
        % calculate the third derivatives of functional (Hessian of
        % gradient) to determine the Lipschitz constant of the Hessian
        % ...L_H...
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %mu=[0.5,0.5];
        
        %parametertestpunkte festlegen
        nr_pts = 6;
        delt = 1/(nr_pts-1);
        mu_v = 0:delt:1;
        
        params.coarse_factor = 4;
        params.cone_const= 1;
        
        model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
        model.decomp_mode = 0;
        model.t = 0;
        
        model_data = gen_model_data(model);
        
        %output functional
        v = model.operators_output(model, model_data);
        L_H = zeros(nr_pts,nr_pts);
        for k=1:nr_pts
            for l=1:nr_pts
                mu = [mu_v(k),mu_v(l)];

                H_der = lin_evol_opt_fd_Hessian_der(model, model_data, mu);

                for i=1:size(H_der,1)
                    for j=1:size(H_der,2)
                        HgradJ(i,j) = v'*H_der{i,j}{1}(:,end);
                    end
                end
                 L_H(k,l) = norm(HgradJ);
                
            end
        end

        
        keyboard;
        
        
    case 119
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % for the non-iterative error estimator calculate
        % gamma-constant
        % this is the norm of the inverse of the hessian at a given
        % parameter value mu. We scan for the largest gamma = gamma_H
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %parametertestpunkte festlegen
        nr_pts = 6;
        delt = 1/(nr_pts-1);
        mu_v = 0:delt:1;
        
        params.coarse_factor = 4;
        params.cone_const= 1;
        
        model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
        model.decomp_mode = 0;
        model.t = 0;
        
        model_data = gen_model_data(model);
        
        %output functional
        v = model.operators_output(model, model_data);
        gamma_H = zeros(nr_pts,nr_pts);
        
        for k=1:nr_pts
            for l=1:nr_pts
                mu = [mu_v(k),mu_v(l)];

                H = lin_evol_opt_fd_Hessian(model, model_data, mu);
                for i=1:size(H,1)
                    for j=1:size(H,2)
                        HgradJ(i,j) = v'*H{i,j}(:,end);
                    end
                end
                 gamma_H(k,l) = norm(inv(HgradJ));
                
            end
        end

        
        keyboard;
        
        
        
        
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for the new non-iterative error estimator calculate the Lipschitz
% constant of the hessian (on a grid in parameter space)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 120
        
         %Lipschitz Konstanten für verschiedene Fälle (zu optimierende
        %Variablen, coarse Faktoren) in arrays gespeichert und als File
        %geschrieben.
        % Das Model sollte dieses File laden und die für sich relevante
        % Lipschitzkonstante auslesen.
        % Zur verbesserung der Genauigkeit können auch loka gültige
        % Lipschitzkonstanten gespeichert werden.
        
        %fix model_type here
        
        %model_type = 'normal';
        %model_type = 'vx_const';
        model_type = 'cone_const';
        %time_type = 'time_const';
        time_type = 'time_var';
        %functional = 'box_mean';
        if (nargin>1)&&(isfield(expar,'functional'))
            switch expar.functional
                case 1
                    functional = 'box_mean';
                case 2
                    functional = 'box_mean_time';
            end
        else
            functional = 'box_mean_time';
        end
        
        params.functional = functional;
        
        if (nargin>1)&&(isfield(expar,'coarse_factor'))
            params.coarse_factor = expar.coarse_factor;
        else
            params.coarse_factor =2;
        end    
        
        %fix here grid density for localized lipschitz constant
        if (nargin>1)&&(isfield(expar,'grid_dens'))
            grid_dens = expar.grid_dens;
        else
            grid_dens=3;
        end
        %grid density of subgrid on each lipschitzgrid element
        sub_grid_dens=1;
        
        if strcmp(time_type,'time_var')
            params.time_const=0;
        else
            params.time_const=1;
        end
        
        switch model_type
            case 'vx_const'
                params.vx_const=1;
                params_to_optimize=[1,0,1];
                grid_params.range={[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1];
            case 'cone_const'
                params.cone_const=1;
                params_to_optimize=[0,1,1];
                grid_params.range={[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1];
            otherwise
                params.normal_model=1;
                params_to_optimize=[1,1,1];
                grid_params.range={[0,1],[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1,1];
        end
        
        model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
        model.t=0;
        model.k=0.25;
        model.compute_derivative_info = 1;
        model.decomp_mode = 0;
        model.optimization.params_to_optmize= params_to_optimize;
        model_data = gen_model_data(model);
        
        %generate lipschitz parameter grid
        lip_H_grid = cubegrid(grid_params);
        
        lipschitz_H_const = zeros(1,lip_H_grid.nelements);
        
        model.decomp_mode=0;
        v=model.operators_output(model, model_data);
        
        matlabpool open
        parfor i=1:lip_H_grid.nelements
            %for every element i in the grid we have to calculate the
            %lipschitz constant
            disp('--------------------------------')
            disp(['sample nr. ',num2str(i)])
            disp('--------------------------------')
            disp('--------------------------------')
            
            %generating subgrid
            sub_par=[];
            sub_par.range=get_ranges_of_element(lip_H_grid,i);
            sub_par.numintervals = sub_grid_dens*ones(1,lip_H_grid.dimension);
            subgrid=cubegrid(sub_par);
            Mtest = subgrid.vertex;
            
           
            %Cl = zeros(sum(params_to_optimize),sum(params_to_optimize));
            L_H = 0;
            for j=1:size(Mtest,1)
                buf_model = set_mu(model, Mtest(j,:));
                %Hessian
                %H_der = filecache_function(@lin_evol_opt_fd_Hessian_der,buf_model, model_data, Mtest(j,:));
                H_der = lin_evol_opt_fd_Hessian_der(buf_model, model_data, Mtest(j,:));
                
                HgradJ = zeros(size(H_der));
                if strcmp(functional, 'box_mean_time')
                    for k=1:size(H_der,1)
                        for l=1:size(H_der,2)
                            buf_model.compute_derivative_info = 0;
                            sim_data_dummy.U = H_der{k,l};
                            sim_data_dummy = lin_evol_opt_output_time_integrated(buf_model, sim_data_dummy,model_data);
                            HgradJ(k,l) = sim_data_dummy.y(end);
                        end
                    end
                
                else    
                   for k=1:size(H_der,1)
                        for l=1:size(H_der,2)
                         HgradJ(k,l) = v'*H_der{k,l}(:,end);
                        end
                    end
                
                end
                
                
                
                
                buf = norm(HgradJ);
                
                 if buf>L_H
                     L_H = buf;
                 end
                 
            end
            
            lipschitz_H_const(i) = L_H;
                               
        end
        
        matlabpool close
               
        lipschitz_H_data=[];
        lipschitz_H_data.lip_H_grid = lip_H_grid;
        lipschitz_H_data.lipschitz_H_const=lipschitz_H_const;
        
        
        if strcmp(model.name_output_functional, 'box_time_int')
            func_type = 'timeInt';
        else
            func_type = '';
        end
        
        
        
        
        save(fullfile(fdir,['lipschitz_data_H',func_type,num2str(model.k),'_',model_type,'_',time_type,'_coarse_',num2str(params.coarse_factor),'.mat']),'lipschitz_H_data')
        
        setpref('Internet','SMTP_Server','wap04.mathematik.uni-stuttgart.de')
        sendmail('dihlmann','am08_experimente fertig','Lipschitz konstante H fertig. Good work Spock!');
        
                       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for the new non-iterative error estimator calculate the gamma
% constant of the hessian (on a grid in parameter space)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 121
        
        
        %fix model_type here
        
        %model_type = 'normal';
        %model_type = 'vx_const';
        model_type = 'cone_const';
        %time_type = 'time_const';
        time_type = 'time_var';
        %functional = 'box_mean';
        functional = 'box_mean_time';
        
        params.functional = functional;
                
        if (nargin>1)&&(isfield(expar,'coarse_factor'))
            params.coarse_factor = expar.coarse_factor;
        else
            params.coarse_factor =2;
        end    
        
        %fix here grid density for localized lipschitz constant
        if (nargin>1)&&(isfield(expar,'grid_dens'))
            grid_dens = expar.grid_dens;
        else
            grid_dens=3;
        end
        %grid density of subgrid on each lipschitzgrid element
        sub_grid_dens=1;
        
        if strcmp(time_type,'time_var')
            params.time_const=0;
        else
            params.time_const=1;
        end
        
        switch model_type
            case 'vx_const'
                params.vx_const=1;
                %params_to_optimize=[1,0,1];
                grid_params.range={[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1];
            case 'cone_const'
                params.cone_const=1;
                %params_to_optimize=[0,1,1];
                grid_params.range={[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1];
            otherwise
                params.normal_model=1;
                %params_to_optimize=[1,1,1];
                grid_params.range={[0,1],[0,1],[0,1]};
                grid_params.numintervals=grid_dens*[1,1,1];
        end
        
        model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
        model.t=0;
        model.k=0.25;
        model.compute_derivative_info = 1;
        model.decomp_mode = 0;
        %model.optimization.params_to_optimize= params_to_optimize;
        model_data = gen_model_data(model);
        
        %generate lipschitz parameter grid
        gamma_H_grid = cubegrid(grid_params);
        
        gamma_H_const = zeros(1,gamma_H_grid.nelements);
        
        model.decomp_mode=0;
        v=model.operators_output(model, model_data);
        
        for i=1:gamma_H_grid.nelements
            %for every element i in the grid we have to calculate the
            %lipschitz constant
            disp('--------------------------------')
            disp(['sample nr. ',num2str(i),' of ',num2str(gamma_H_grid.nelements)])
            disp('--------------------------------')
            disp('--------------------------------')
            %generating subgrid
            sub_par.range=get_ranges_of_element(gamma_H_grid,i);
            sub_par.numintervals = sub_grid_dens*ones(1,gamma_H_grid.dimension);
            subgrid=cubegrid(sub_par);
            Mtest = subgrid.vertex;
            
           
            %Cl = zeros(sum(params_to_optimize),sum(params_to_optimize));
            gamma_H = 0;
            for j=1:size(Mtest,1)
                model = set_mu(model, Mtest(j,:));
                %Hessian
                H = filecache_function(@lin_evol_opt_fd_Hessian,model, model_data, Mtest(j,:));
                
                
                if strcmp(functional, 'box_mean_time')
                    for k=1:size(H,1)
                        for l=1:size(H,2)
                            model.compute_derivative_info=0;
                            sim_data_dummy.U = H{k,l};
                            sim_data_dummy = lin_evol_opt_output_time_integrated(model, sim_data_dummy,model_data);
                            HJ(k,l) = sim_data_dummy.y(end);
                        end
                    end
                else    
                    for k=1:size(H,1)
                        for l=1:size(H,2)
                            HJ(k,l) = v'*H{k,l}(:,end);
                        end
                    end
                end
                
                buf = norm(inv(HJ));
                
                 if buf>gamma_H
                     gamma_H = buf;
                 end
                 
            end
            
            gamma_H_const(i) = gamma_H;
                               
        end
               
        gamma_H_data=[];
        gamma_H_data.gamma_H_grid = gamma_H_grid;
        gamma_H_data.gamma_H_const=gamma_H_const;
        
        
        if strcmp(model.name_output_functional, 'box_time_int')
            func_type = 'timeInt';
        else
            func_type = '';
        end
        
      
        save(fullfile(fdir,['gamma_data_H',func_type,num2str(model.k),'_',model_type,'_',time_type,'_coarse_',num2str(params.coarse_factor),'.mat']),'gamma_H_data')

        %%%%%%%%%%%%%%%%%%%%%%%5
        %generate reduced data for second functional
        %%%%%%%%%%%%%%%%%
    case 122
        
        %in verzeichnis wechseln wo reduzierte basis liegt
        cd /home/dihlmann/matlab/rbmatlab-core/rbasis/lin_evol_opt/convdiff/coarse2/basisNmax200
        
        load(fullfile(fdir,'coarse2','basisNmax200','complete_TP_basis_coarse2.mat'));
        %load(fullfile(fdir,'coarse16','complete_TP_basis_coarse16.mat'));
        %model fields für zweites funktional ändern
        %set 2nd functional fields
         model.verbose = 4;
        model.base_model.verbose = 4;
        model.name_output_functional = 'box_mean_time';
        model.get_output = @lin_evol_opt_output_time_integrated;
        model.output_function_ptr = @output_function_box_mean;
        model.sbox_xmin = 0.5;
        model.sbox_xmax = 1;
        model.sbox_ymin = 0.5;
        model.sbox_ymax = 1;
        model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.s_l2norm = sqrt(4)*sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.s_l1norm = sqrt(4);
        model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.sbox_xmin = 0.5;
        model.base_model.sbox_xmax = 1;
        model.base_model.sbox_ymin = 0.5;
        model.base_model.sbox_ymax = 1;
        model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.base_model.s_l2norm = sqrt(4)*sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.s_l1norm = sqrt(4);
        model.base_model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.base_model.sbox_xmin = 0.5;
        model.base_model.base_model.sbox_xmax = 1;
        model.base_model.base_model.sbox_ymin = 0.5;
        model.base_model.base_model.sbox_ymax = 1;
        model.base_model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt/2));
        model.base_model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt*3/4));
        model.base_model.base_model.s_l2norm = sqrt(4)*sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.base_model.s_l1norm = sqrt(4);
        %set optimzation parameters
        model.optimization.tol = 0.00001;
        model.base_model.optimization.tol = 0.00001;
        model.base_model.base_model.optimization.tol =0.00001;
        
       %new optimizer
        model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.base_model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        %new stepsize rule
        model.optimization.stepsize_rule = 'armijo';
        model.optimization.sigma = 0.5;
        model.optimization.stepsize_factor = 0.75;
        model.optimitation.initial_step = 50;
        model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.optimization.sigma = 0.5;
        model.base_model.optimization.stepsize_factor = 0.75;
        model.base_model.optimitation.initial_step = 50;
        
        %neue reduzierte daten generieren
        disp('generating reduced data for second functional')
        reduced_data = gen_reduced_data(model, detailed_data);
        disp('Saving reduced data')
        save -v7.3 reduced_data_complete_2ndfunctional.mat model reduced_data
        
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% generiere Landschaftsplot für funktional 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 123
        
        %set landscape grid options (density...)
        if (nargin>1)&&(isfield(expar,'num_pts'))
            num_pts = expar.num_pts;
        else
            num_pts = 25;
        end
        p_range = [0,1];
        delt = (p_range(2)-p_range(1))/(num_pts-1);
        mu_test = p_range(1):delt:p_range(2);
        
        %introduce variables
        n_grid = length(mu_test);
        x_vec = mu_test;
        y_vec = mu_test;
        func1_det = zeros(n_grid,n_grid);
        func1_red = zeros(n_grid,n_grid);
        
        %lade model und reduced data
        disp('loading data...')
        load(fullfile(fdir,'coarse2','reduced_data_complete_N200.mat'));
        %load(fullfile(fdir,'coarse16','reduced_data_complete_coarse16.mat'));
        
        
      disp('reduced simulations on the grid...') 
       
       %reduced_simulation on grid
       %matlabpool open
        for i=1:n_grid
            disp(['Zwischenstand: ',num2str(i),' von ',num2str(n_grid),'.']);
            for j=1:n_grid
                model_buf = set_mu(model,[x_vec(i),y_vec(j)]);
                func1_red(i,j) = model.optimization.objective_function(model_buf, reduced_data);
            end
       end
       %matlabpool close
       
       save(fullfile(fdir,['Landschaft_func1_reduced_coarse',num2str(model.base_model.coarse_factor),...
           '_numpts_',num2str(num_pts),'.mat']),'func1_red','x_vec','y_vec')
       
       
        %generiere model_data
        model_data = gen_model_data(model);
       
        disp('detailed simulations on grid...')
        %detailed simulations on grid
        matlabpool open
       parfor i=1:n_grid
           disp(['Zwischenstand: ',num2str(i),' von ',num2str(n_grid),'.']);
            for j=1:n_grid
                model_buf = set_mu(model,[x_vec(i),y_vec(j)]);
                func1_det(i,j) = model.optimization.objective_function(model_buf, model_data);
            end
       end
       matlabpool close
       save(fullfile(fdir,['Landschaft_func1_coarse',num2str(model.base_model.coarse_factor),...
           '_numpts_',num2str(num_pts),'.mat']),'func1_red','func1_det','x_vec','y_vec')
       
        
        

       
       
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% generiere Landschaftsplots für funktional 3 (nur detailliert)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 125
        
        
        %set landscape grid options (density...)
        if (nargin>1)&&(isfield(expar,'num_pts'))
            num_pts = expar.num_pts;
        else
            num_pts = 8;
        end
        p_range = [0,1];
        delt = (p_range(2)-p_range(1))/(num_pts-1);
        mu_test = p_range(1):delt:p_range(2);
        
        %introduce variables
        n_grid = length(mu_test);
        x_vec = mu_test;
        y_vec = mu_test;
        func3_det = zeros(n_grid,n_grid);
        %func3_red = zeros(n_grid,n_grid);
        
        %lade model und reduced data
        disp('loading data...')
        load(fullfile(fdir,'coarse2','basisNmax200','reduced_data_complete_2ndfunctional.mat'));
        %load(fullfile(fdir,'coarse16','reduced_data_complete_coarse16_2ndfunctional.mat'));
        
        %change fields in model for functional 2
        model.verbose = 4;
        model.base_model.verbose = 4;
        model.name_output_functional = 'box_mean_time';
        model.get_output = @lin_evol_opt_output_time_integrated;
        model.output_function_ptr = @output_function_box_mean;
         model.sbox_xmin = 0;
        model.sbox_xmax = 1;
        model.sbox_ymin = 0.5;
        model.sbox_ymax = 1;
        model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt*0.75));
        model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt));
        model.s_l2norm = sqrt(4)*sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.s_l1norm = sqrt(4);
        model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.sbox_xmin = 0;
        model.base_model.sbox_xmax = 1;
        model.base_model.sbox_ymin = 0.5;
        model.base_model.sbox_ymax = 1;
        model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt*0.75));
        model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt));
        model.base_model.s_l2norm = sqrt(4)*sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.s_l1norm = sqrt(4);
        model.base_model.base_model.name_output_functional = 'box_mean_time';
        model.base_model.base_model.get_output = @lin_evol_opt_output_time_integrated;
        model.base_model.base_model.output_function_ptr = @output_function_box_mean;
        model.base_model.base_model.sbox_xmin = 0;
        model.base_model.base_model.sbox_xmax = 1;
        model.base_model.base_model.sbox_ymin = 0.5;
        model.base_model.base_model.sbox_ymax = 1;
        model.base_model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt*0.75));
        model.base_model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt));
        model.base_model.base_model.s_l2norm = sqrt(4)*sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
        model.base_model.base_model.s_l1norm = sqrt(4);
        %set optimzation parameters
        model.optimization.tol = 0.0001;
        model.base_model.optimization.tol = 0.0001;
        model.base_model.base_model.optimization.tol =0.0001;
        
       
%        disp('reduced simulations on the grid...') 
%        
%        %reduced_simulation on grid
%        matlabpool open 
%        parfor i=1:n_grid
%             disp(['Zwischenstand: ',num2str(i),' von ',num2str(n_grid),'.']);
%             for j=1:n_grid
%                 model_buf = set_mu(model,[x_vec(i),y_vec(j)]);
%                 func2_red(i,j) = model.optimization.objective_function(model_buf, reduced_data);
%             end
%        end
%        matlabpool close
%        
%        save(fullfile(fdir,['Landschaft_func2_reduced_coarse',num2str(model.base_model.coarse_factor),...
%            '_numpts_',num2str(num_pts),'.mat']),'func2_red','x_vec','y_vec')
%        
       
        %generiere model_data
        model_data = gen_model_data(model);
       
        %detailed simulations on grid
        matlabpool open
       parfor i=1:n_grid
            for j=1:n_grid
                disp(['Zwischenstand: ',num2str(i),' von ',num2str(n_grid),'.']);
                model_buf = set_mu(model,[x_vec(i),y_vec(j)]);
                func3_det(i,j) = model.optimization.objective_function(model_buf, model_data);
            end
       end
       matlabpool close
       save(fullfile(fdir,['Landschaft_func3_coarse',num2str(model.base_model.coarse_factor),...
           '_numpts_',num2str(num_pts),'.mat']),'func3_det','x_vec','y_vec')
       

       
       
       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% mit verschiedenen Funktionalen herumspielen und sie plotten
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 126
       
         params.coarse_factor = 2;
         params.cone_const = 1;
         params.functional = 'box_mean_time2';
%         
         model = convection_diffusion_fv_output_opt_model_rect_normed_outflow(params);
         model.t=0;
         model.k=0.25;
         model.compute_derivative_info = 1;
         model.decomp_mode = 0;
         model_data = gen_model_data(model);
%         
%         %set landscape grid options (density...)
%         if (nargin>1)&&(isfield(expar,'num_pts'))
%             num_pts = expar.num_pts;
%         else
%             num_pts = 7;
%         end
%         p_range = [0,1];
%         delt = (p_range(2)-p_range(1))/(num_pts-1);
%         mu_test = p_range(1):delt:p_range(2);
%         
%         
%         %Ouput functional
%         v=model.operators_output(model, model_data);
%         
%         %landscape values - initialization
%         J = zeros(num_pts,num_pts); %functional values
%         
%         %loop over parameter grid
%         for i=1:num_pts
%             for j=1:num_pts
%                
%                mu = [mu_test(i),mu_test(j)];
%                model = set_mu(model, mu);
%                sim_data = detailed_simulation(model, model_data);
%                
%                 %save results
%                 J(i,j) = sim_data.y(end);
%                 
%             end
%         end
%         
%         figure;surf(mu_test,mu_test,J);title(['J_',num2str(params.coarse_factor)]);
%         
%         
%         keyboard
%         
        %part 2
        
        
         %adapt reduced data to new functional
         disp('loading detailed data')
          cd /home/dihlmann/matlab/rbmatlab-core/rbasis/lin_evol_opt/convdiff/coarse2
         load(fullfile(fdir,'coarse2','basisNmax200','complete_TP_basis_coarse2.mat'));
%         
        %optimization
        model.verbose = 4;
%         model.base_model.verbose = 4;
%         model.name_output_functional = 'box_mean_time';
%         model.get_output = @lin_evol_opt_output_time_integrated;
%         model.output_function_ptr = @output_function_box_mean;
%          model.sbox_xmin = 0;
%         model.sbox_xmax = 1;
%         model.sbox_ymin = 0.5;
%         model.sbox_ymax = 1;
%         model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt*0.75));
%         model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt));
%         model.s_l2norm = sqrt(4)*sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
%         model.s_l1norm = sqrt(4);
%         model.base_model.name_output_functional = 'box_mean_time';
%         model.base_model.get_output = @lin_evol_opt_output_time_integrated;
%         model.base_model.output_function_ptr = @output_function_box_mean;
%         model.base_model.sbox_xmin = 0;
%         model.base_model.sbox_xmax = 1;
%         model.base_model.sbox_ymin = 0.5;
%         model.base_model.sbox_ymax = 1;
%         model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt*0.75));
%         model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt));
%         model.base_model.s_l2norm = sqrt(4)*sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
%         model.base_model.s_l1norm = sqrt(4);
%         model.base_model.base_model.name_output_functional = 'box_mean_time';
%         model.base_model.base_model.get_output = @lin_evol_opt_output_time_integrated;
%         model.base_model.base_model.output_function_ptr = @output_function_box_mean;
%         model.base_model.base_model.sbox_xmin = 0;
%         model.base_model.base_model.sbox_xmax = 1;
%         model.base_model.base_model.sbox_ymin = 0.5;
%         model.base_model.base_model.sbox_ymax = 1;
%         model.base_model.base_model.stime_start = find(model.base_model.save_time_indices ==floor(model.base_model.nt*0.75));
%         model.base_model.base_model.stime_stop = find(model.base_model.save_time_indices == floor(model.base_model.nt));
%         model.base_model.base_model.s_l2norm = sqrt(4)*sqrt((model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices));
%         model.base_model.base_model.s_l1norm = sqrt(4);
%         %set optimzation parameters
        model.optimization.tol = 0.0001;
        model.base_model.optimization.tol = 0.0001;
        model.base_model.base_model.optimization.tol =0.0001;
        
        model.optimization.init_params =[0.3,0.4];
        model.base_model.optimization.init_params =[0.3,0.4];
        model.base_model.base_model.optimization.init_params =[0.3,0.4];
        
        model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        model.base_model.base_model.optimization.optimizer = @gradient_opt_non_iter_err;
        
        model.optimization.stepsize_rule = 'armijo';
        model.optimization.sigma = 0.8;
        model.optimization.stepsize_factor = 0.75;
        model.optimization.initial_stepsize = 10;
        model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.optimization.sigma = 0.8;
        model.base_model.optimization.stepsize_factor = 0.75;
        model.base_model.optimization.initial_stepsize = 10;
        model.base_model.base_model.optimization.stepsize_rule = 'armijo';
        model.base_model.base_model.optimization.sigma = 0.8;
        model.base_model.base_model.optimization.stepsize_factor = 0.75;
        model.base_model.base_model.optimization.initial_stepsize = 10;
        
        model.gamma_H = @(model)dummy_constant(model,100);
        model.base_model.gamma_H = @(model)dummy_constant(model,100);
        model.base_model.base_model.gamma_H =@(model)dummy_constant(model,100);
        model.L_H = @(model)dummy_constant(model,100);
        model.base_model.L_H =@(model)dummy_constant(model,100);
        model.base_model.base_model.L_H = @(model)dummy_constant(model,100);
       
        model.optimization.allow_constant_calculation = 1;
        model.base_model.optimization.allow_constant_calculation = 1;
        model.base_model.base_model.optimization.allow_constant_calculation = 1;
        
       %save reduced data
    disp('generating reduced data for second functional')
        reduced_data = gen_reduced_data(model, detailed_data);
        disp('Saving reduced data')
        save -v7.3 reduced_data_complete_3rdfunc.mat model reduced_data
        
       
        %optimzation and constant calulation
        opt_data = optimize(model,reduced_data);
        
        save(fullfile(fdir,'results_case126reduced.mat'),'opt_data');
        %keyboard;
         setpref('Internet','SMTP_Server','wap04.mathematik.uni-stuttgart.de')
        sendmail('dihlmann','optimierung case 126 abgeschlossen','look at it');
        
        %calculate gamma and lip at the optimum
%         mu = opt_data.optimal_params;
%         
%         gamma = calculate_gamma_const(model,model_data,mu)
%         lip = calculate_Lipschitz_const(model, model_data,mu)
%         
%         save(fullfile(fdir,'results_case126.mat'),'opt_data','gamma','lip');
%         
end