Algorithm.m

classdef Algorithm < handle 
%ALGORITHM The main implementation of the LRFG algorithm
%   This algorithm is an implementation of the low-rank factor greedy
%   algorithm as it was presented in the paper by A. Schmidt and B.
%   Haasdonk (2016). It is a generic class which can be applied to all OOP
%   models that implement a certain interface.
%

    properties(GetAccess=protected, SetAccess=protected)
        model;         % Pointer to the ReducedModel class
        model_data;    % The reference to the ModelData class
    end
    properties
        M_train;        % Save the M_train class for subsequent use
        
        error_indicator = 'normalized_residual'; % Use the normalized residual as 
                                        % error indicator
                                        
        orthonormalize_E = true;        % Orthonormalize the basis wrt. to the mass matrix
    end
    
    methods
        function this = Algorithm(model, model_data, M_train)
            % ALGORITHM Constructor for the LRFG algorithm
            %
            this.model = model;
            this.model_data = model_data;
            
            if ~exist('M_train', 'var')
                M_train = ParameterSampling.Uniform(9);
            end
            
            this.M_train = M_train;
        end
        
        function set.M_train(this, M_train)
            if ~isa(M_train, 'ParameterSampling.Interface')
                error('M_train must be of type ParameterSampling.Interface')
            end
            
            this.M_train = M_train;
        end
        
    end
    
    methods(Access=protected)

        function this = initialize_basis(this, model, model_data)
            % This function should initialize the basis with a first
            % element. The default implementation takes the solution
            if ~isempty(this.RB_V) || ~isempty(this.RB_W)
                error('The basis is not empty, so the initialize_basis method cannot be called');
            end
            
            mus = this.M_train.sample;
            model.set_mu(mus(1,:));
            
            sim = detailed_simulation(model, model_data);
            
            [V,~,~] = svds(sim.Z, 1);
            this.RB_V = V;
            this.RB_W = V;
        end
        
    end
    
    % =====================================================================
    % Public methods
    methods
                        
        % =================================================================
        function [V,W,info] = run(this)
            % RUN 
            %    Implementation of the LRFG algorithm.
            %    This function actually performs the algorithm and fills
            %    the corresponding fields in detailed_data
            %
            if init_required(this.M_train)
              init_sample(this.M_train, this.model);
            end
            
            % Cache the mus since we need them a few times below and
            % abbreviate this.rmodel 
            mus = this.M_train.sample;
            model = copy(this.model);
                        
            t_start = tic;
            
            % This class will later be populated into the DetailedData.info
            % property
            info = Greedy.LRFG.BasisInfo;
            
            info.greedy_tolerance = model.RB_greedy_tolerance;
            gtol = info.greedy_tolerance;
            this.orthonormalize_E = model.RB_orthonormalize_E;
            info.orthonormalize_E = this.orthonormalize_E;
            this.error_indicator = model.RB_error_indicator;
            info.error_indicator = this.error_indicator;
            info.pod_max_extension = model.RB_pod_max_extension;
            
            % To store the error decay:
            errors = ones(size(mus, 1), 1);

            % The reduced basis will be stored here:
            RB_V = [];
            RB_W = [];

            % Certain arrays that contain statistical information about the
            % LRFG run:
            chosen_mu = [];
            error_decay = [];
            Z_sizes = [];
            vectors_added = [];
            num_solves = 0;
            
            % Some additional helper variables:
            i = 1; z0 = 0; sims = cell(size(mus,1),1);

            % For the detailed simulation time
            dsim_times = [];
            dsim_times_real = [];
    
            % Inner tolerance for the POD:
            podtol = model.RB_pod_tolerance;
            info.pod_tolerance = podtol;

            % Grab the mass matrix, such that we can normalize
            E = speye(model.n);
            if this.orthonormalize_E
                E = model.mass_matrix(this.model_data);
            end
            n = model.n;

            display('================================================================')
            display('Starting Basis Generation')
            display('----------------------------------------------------------------')
            display('Basis generation settings:')
            display(['  - Number of parameters               ', num2str(length(mus))])
            display(['  - Inner tolerance                    ', num2str(podtol)])
            display(['  - Greedy tolerance                   ', num2str(gtol)])
            display(['  - System dimension                   ', num2str(n)])
            display('+------+---------------+-------+------------+--------+----------+')
            display('|    i |      max. err.|   idx.|    Z size  |  added |    time  |')
            display('+------+---------------+-------+------------+--------+----------+')

            % Here starts the main LRFG loop:
            while max(errors) > gtol
                if isempty(RB_V) && isempty(RB_W)
                    % Use the first element in the training set for the
                    % initalization of the basis
                    theta = 1;
                    max_err = Inf;
                else
                    % Select the worst approximated element
                    [max_err, theta] = max(errors);
                end

                % Save the chosen parameter index:
                chosen_mu = [chosen_mu theta]; 

                % Calculate the full dimensional solution
                model.set_mu(mus(theta,:));
                t = tic;
                [dsim,cached] = filecache_function(@detailed_simulation, model, this.model_data);
                time = (cached<1)*toc(t);
                num_solves = num_solves + 1;
                dsim_times_real = [dsim_times_real time];
                dsim_times = [dsim_times dsim.time];

                Z_sizes = [Z_sizes size(dsim.Z, 2)];
                if ~isempty(RB_V) && ~isempty(RB_W)
                    Z = dsim.Z - RB_W*((RB_V'*RB_W)\(RB_V'*dsim.Z));
                else
                    Z = dsim.Z;
                end

                % Perform the basis extension:
                [Phi, Lambda, ~] = svd(Z, 'econ');
                dLambda = sqrt(diag(Lambda));
                sumS = sum(dLambda);
                for k = 1:length(dLambda)
                    eE = sum(dLambda(1:k))/sumS;
                    if eE >= podtol || (info.pod_max_extension > 0 && k == info.pod_max_extension)
                        break
                    end
                end
                
                Vextend = Phi(:, 1:k);
                c = size(Vextend, 2);

                % Just make sure the basis is orthogonal in cases where numerical
                % errors occur:
                RB_W = orthonormalize_gram_schmidt([RB_W Vextend], E);
                RB_V = RB_W;

                vectors_added = [vectors_added c];

                % Recalculate the errors:
                [errs_normalized,errs_absolute] = this.calc_error_indicator(RB_W, RB_V, mus);

                errors = errs_normalized;
                switch info.error_indicator
                    case 'residual'
                        errors = errs_absolute;
                    case 'normalized_residual'
                        errors = errs_normalized;
                    case 'relative_change'
                        if i == 1
                            initial_residuals = errs_absolute;
                            errors = ones(1,length(errs_absolute));
                        else
                            errors = errs_absolute ./ initial_residuals;
                        end
                    otherwise
                        error(['Error indicator mode "', info.error_indicator, '" not known'])
                end
                
                error_decay = [error_decay max(errors)];
                error_history_normalized{i} = errs_normalized;
                error_history_absolute{i} = errs_absolute;
                
                display(sprintf('| %4d | %.7d | %5d | %10d | %6d | %7f |', i, max(errors), theta, size(dsim.Z,2), c, time))
                i = i + 1;
            end
            display('+------+---------------+-------+------------+--------+----------+')

            V = RB_V;
            W = RB_W;
            
            % Fill in the remaining information:
            info.Z_sizes = Z_sizes;
            info.vectors_added = vectors_added;
            info.detailed_simulation_times = dsim_times;
            info.detailed_simulation_times_real = dsim_times_real;
            info.chosen_mu = chosen_mu;
            info.error_history_normalized = error_history_normalized;
            info.error_history_absolute = error_history_absolute;
            info.error_decay = error_decay;
            info.num_solves = num_solves;
            info.M_train = this.M_train;
            
            % End of basisgen: save the complete offline time
            info.runtime = toc(t_start);
        end
        
        
        % Calculate the errors for the basis stored in the matrix RB
        function [errs_normalized, errs_absolute, avg_time] = calc_error_indicator(this, RB_W, RB_V, mus)
            errs_normalized = zeros(size(mus,1),1);
            errs_absolute = zeros(size(mus,1),1);

            % Create a DetailedData object and a ReducedData object for
            % temporary use:
            pt = this.model.problem_type();
            dd = eval([pt, '.DetailedData(this.model, this.model_data)']);
            dd.RB_V = RB_V;
            dd.RB_W = RB_W;

            rd = gen_reduced_data(this.model, dd);
            t = tic;
            model = copy(this.model);
            
            for i = 1:size(mus,1)
                tmpModel = model;
                tmpModel.set_mu(mus(i, :));
                
                % Disable the error estimator but enable the calculation of
                % the residual:
                tmpModel.enable_error_estimator = false;
                tmpModel.calc_residual = true;
  
                % Simulate and get the errors:
                rsim = rb_simulation(tmpModel, rd);
                errs_normalized(i) = rsim.nresidual;
                errs_absolute(i) = rsim.residual;
            end

            avg_time = toc(t)/size(mus,1);
        end
    end
end