riccati_gen_detailed_data.m

function detailed_data = riccati_gen_detailed_data(model, model_data)
%function detailed_data = riccati_gen_detailed_data(model, model_data)
% 
%  function generates the reduced basis
%  the result is stored in the detailed_data
%
%  This function generates the following data structure:
%    detailed_data.
%      RB                  The reduced basis
%      RB_info             Information about the basis
%      
%  Andreas Schmidt, 2015
%
snapshots = [];
training_mus = {};
chosen_mu = [];

if ~isfield(model, 'RB_generation_mode')
    error(['Please specify model.RB_generation_mode '])
end

mode = model.RB_generation_mode;
dsim_mode = 'riccati';

switch mode
    case 'greedy_uniform_fixed'
        if ~isfield(model, 'RB_numintervals')
            error('Please specify model.RB_numintervals')
        end
        mu_ranges = model.mu_ranges;
        par.numintervals = model.RB_numintervals*ones(length(mu_ranges),1);
        par.range = model.mu_ranges;
        MMesh0 = cubegrid(par);
        training_mus = get(MMesh0, 'vertex');
        algorithm = @greedy;
    case 'greedy_fixed'
        if ~isfield(model, 'training')
            error('Please specify model.training')
        end
        training_mus = model.training;
        algorithm = @greedy;
    case 'greedy_rand'
        training_mus = rand_uniform(model.RB_samplepoints, model.mu_ranges)';
        algorithm = @greedy;
    case 'greedy_rand_log'
        training_mus = rand_log(model.RB_samplepoints, model.mu_ranges)';
        algorithm = @greedy;
    case 'pod'
        algorithm = @greedy;
        dsim_mode = 'pod';
        training_mus = rand_uniform(model.RB_samplepoints, model.mu_ranges)';
        
    otherwise
        error('Basis generation mode unknown')
end
detailed_data = algorithm(model, model_data, training_mus, dsim_mode);

end



%% First version: Greedy for detailed system
function detailed_data = greedy(model, model_data, training_mus, varargin)
% Implementation of the greedy algorithm
    RB_info.M_train = training_mus;
    
    bgenStart = tic;
    % Outer tolerance
    o_tol = 10e-6;
    if isfield(model, 'RB_o_tol')
        o_tol = model.RB_o_tol;
    end
    RB_info.o_tol = o_tol;
    
    % Errors (dummy values set for initial calculation)
    errors = ones(size(training_mus, 1), 1);
    
    % The reduced basis will be stored here:
    RB = [];
    
    % History of the error decay over the whole parameter domain:
    error_history = {};
    chosen_mu = [];
    
    % Maximum error decay:
    error_decay = [Inf];
    
    % The sizes of the low rank factors
    Z_sizes = [];
    
    % Number of inner iterations
    inner_elements_added = [];

    i = 1;
    z0 = 0;
    Z = cell(size(training_mus,1));

    % For the detailed simulation time
    dsim_times = [];
    
    mode = 'riccati';
    if nargin > 3
        mode = varargin{1};
    end
    
    % Adaptive algorithm or fixed-1 ?
    basis_extension_mode = model.RB_extension_mode;
    addmultiple = false;
    if strcmp(basis_extension_mode, 'addmultiple')
        addmultiple = true;
        i_tol = 0.8;
        if isfield(model, 'RB_i_tol')
            i_tol = model.RB_i_tol;
        end
        maxperstep = model.RB_extension_max_per_step;
        RB_info.i_tol = i_tol;
    end
    RB_info.addmultiple = addmultiple;
    [E,~,~,~] = riccati_assemble_system_matrices(model,model_data);
    
    while max(errors) > o_tol
        display(['LRFG, Iteration ', num2str(i)])
       
        if isempty(RB)
            % Use the first element in the training set for the
            % initalization of the basis
            theta = 1;
            maxErr = Inf;
        else
            % Select the worst approximated element
            [maxErr, theta] = max(errors);
        end
        
        display([' - Element ', num2str(theta), ' with error ', num2str(maxErr)])
        
        % Save the error decay
        chosen_mu = [chosen_mu theta];

        % Copy the model temporarily and perform a detailed simulation to
        % get the new Low Rank factor Z and cache it in the cell Z
        display(' - Getting detailed simulation')
        if isempty(Z{theta})
            mu = training_mus(theta,:);
            tmpModel = set_mu(model, mu);
            
            tStart = tic;
            if strcmp(mode, 'riccati')
                dsim = detailed_simulation(tmpModel, model_data);
            elseif strcmp(mode, 'pod')
                snapshots = [];
                [E,A,~,C] = riccati_assemble_system_matrices(tmpModel, model_data);
                parfor l = 1:size(C,1)
                    c = C(l,:)';
                    options = odeset('Jacobian', A, 'Mass', E');
                    [T,Y] = ode23s( @(t,x) (A'*x), [0,5], E\c, options);
                    snapshots = [snapshots Y'];
                end
                dsim.Z = snapshots;
            else
                error('Detailed data generation mode not known')
            end
            dsim_times = [dsim_times toc(tStart)];
            
            Z{theta} = dsim.Z;
        end
        thisZ = Z{theta};
        
        Z_sizes = [Z_sizes size(thisZ,2)];
        display([' - Low rank factor has ', num2str(size(thisZ,2)), ' elements'])
        if z0 == 0
            z0 = norm(thisZ(:,1));
        end
        
        % Extend the basis, i.e. run the inner greedy loop
        %[RB,c,ierror,newZ] = riccati_greedy_extend_basis(RB,thisZ,itol,z0);
        if length(RB) > 0
            thisZ = thisZ - RB*(RB'*thisZ);
        end
        
        if ~addmultiple
            [Vextend,S,~] = svds(thisZ,1);
        else
            [Vextend,S,~] = svd(thisZ,'econ');
            S = diag(S).^2;
            sumS = sum(S);
            
            for i = 1:length(diag(S))
                E = 1 - sum(S(1:i))/sumS;
                if E < i_tol || i == maxperstep
                    break
                end
            end
            Vextend = Vextend(:, 1:i);
        end

        c = size(Vextend, 2);
        ierror = max(diag(S));
       
        % Just make sure the basis is orthogonal in cases where numerical
        % errors occur:
        RB = orthonormalize_gram_schmidt([RB Vextend]);
        
        inner_elements_added = [inner_elements_added c];
        
        display([' - Added ', num2str(c), ' new basis vectors']);
        display([' - New basis size: ', num2str(size(RB,2))]);

        % Only recalculate if elements have been added
        % If no elements were added the same error bound can be used again
        if c>0
            display(' - Recalculating errors')
            errors = greedy_calc_errors(model, RB, training_mus);
        end
        
        last_error = error_decay(end);
        error_decay = [error_decay max(errors)];
        error_history{i} = errors;
        
        i = i + 1;
    end
    RB_info.basisgen_runtime = toc(bgenStart);
    
    detailed_data.RB = RB;
    
    % Start gamma procedure:
    gamma_mode = model.RB_gamma_mode;
    gamma.mode = gamma_mode;
    gamma_timer = tic;
    % Check the mode for gamma:
    if strcmp(gamma_mode, 'max')
        error('To implement')
    elseif strcmp(gamma_mode, 'kernel_interpolation')
        greedy_interpol_count = 100;
        if isfield(model, 'RB_gamma_samplesize')
            greedy_interpol_count = model.RB_gamma_samplesize;
        end
        mus = rand_uniform(greedy_interpol_count, model.mu_ranges)';
        greedy_gammas = zeros(greedy_interpol_count, 1);
        runtimes = zeros(greedy_interpol_count,1);
        parfor_progress(size(mus,1));
        parfor i = 1:size(mus, 1)
            mu = mus(i,:);
            tmpModel = set_mu(model, mu);
            dsim = detailed_simulation(tmpModel, model_data);
            gamma_timing = tic;
            greedy_gammas(i) = riccati_gamma(tmpModel, model_data, dsim);
            runtimes(i) = toc(gamma_timing);
            parfor_progress
        end
        
        gamma.interpol_mus = mus;
        gamma.interpol_gammas = greedy_gammas;
        gamma.kernel = 'gauss';
        gamma.runtimes = runtimes;
    elseif strcmp(gamma_mode, 'interpol')
        detailed_data.gamma_mode = 'interpol';
        display('Starting calculation of gamma values')
        
        % Build the mu-grid:
        Nmu = 6;
        ranges = model.mu_ranges;
        
        mug = cell(1,length(ranges));
        for j = 1:length(ranges)
            mug{j} = linspace(ranges{j}(1), ranges{j}(2), Nmu);
        end
        mu_grid = cell(1,numel(mug));
        
        [mu_grid{:}] = ndgrid(mug{:});
        gamma_values = zeros(size(mu_grid{1}));
        parfor_progress(numel(mu_grid{1}));
        reduced_data = gen_reduced_data(model, detailed_data);
        
        parfor k = 1:numel(mu_grid{1})
            mu = zeros(length(ranges),1);
            for j = 1:length(ranges)
                mu(j) = mu_grid{j}(k);
            end
            
            tmpModel = set_mu(model, mu);
            sim = rb_reconstruction(tmpModel, detailed_data, rb_simulation(tmpModel, reduced_data));
            gamma_values(k) = riccati_gamma(tmpModel, model_data, sim);
            
            parfor_progress;
        end
        parfor_progress(0);
        
        gamma.values = gamma_values;
        gamma.mu_grid = mu_grid;
        gamma.interpolant = griddedInterpolant( mu_grid{:}, gamma_values, 'linear' );
    end
    detailed_data.gamma = gamma;
    RB_info.gamma_runtime = toc(gamma_timer);
    
    % Fill the structure with all information from the basis building
    % process.
    inner.Z_sizes = Z_sizes;
    inner.iterations = inner_elements_added;
    RB_info.inner = inner;
    
    outer.error_decay = error_decay;
    outer.chosen_mu = chosen_mu;
    outer.tol = o_tol;
    outer.dsim_timings = dsim_times;
    RB_info.outer = outer;
    
    detailed_data.RB_info = RB_info;
end

% Calculate the errors for the basis stored in the matrix RB
function [errs, avg_time] = greedy_calc_errors(model, RB, training_mus)
    errs = zeros(size(training_mus,1),1);
    dd.RB = RB;
    dd.tol = 10e-12;
    rd = gen_reduced_data(model, dd);
    t = tic;
    parfor i = 1:size(training_mus,1)
        tmpModel = model;
        tmpModel = set_mu(tmpModel, training_mus(i, :));
        
        rsim = rb_simulation(tmpModel, rd);
        errs(i) = rsim.P_residual_normalized;
    end
    avg_time = toc(t);
end