ReducedData.m

classdef ReducedData < AbstractModel.ReducedData
    
    properties
        A_comp;
        B_comp;
        C_comp;
        E_comp;
        R_comp;
        Q_comp;
        
        estim;
        
        gamma_function;
        gamma_mode;
        
        RB_V;
        RB_W;
        
        RB_size;
    end
    
    methods
        function obj = ReducedData(model, detailed_data)
            % REDUCEDDATA Constructor for the class
            model_data = detailed_data.model_data;
            
            % Extract the reduced basis matrices
            V = detailed_data.RB_V;
            W = detailed_data.RB_W;
            
            obj.RB_V = V;
            obj.RB_W = W;
            obj.RB_size = size(W,2);
            
            for k = 1:length(model_data.A_comp)
                obj.A_comp{k} = W'*model_data.A_comp{k}*V;
            end

            for k = 1:length(model_data.B_comp)
                obj.B_comp{k} = W'*model_data.B_comp{k};
            end

            for k = 1:length(model_data.C_comp)
                obj.C_comp{k} = model_data.C_comp{k}*V;
            end
            
            for k = 1:length(model_data.E_comp)
                obj.E_comp{k} = W'*model_data.E_comp{k}*V;
            end
            
            obj.R_comp = model_data.R_comp;
            obj.Q_comp = model_data.Q_comp;
            
            % Calculate all the reduced quantities:
            estim.M1a = {};
            for i = 1:length(model_data.C_comp)
                for j = 1:length(model_data.A_comp)
                    estim.M1a{i,j} = model_data.C_comp{i}*model_data.A_comp{j}'*W;
                end
            end
            for i = 1:length(model_data.E_comp)
                for j = 1:length(model_data.C_comp)
                    estim.M1b{i,j} = W'*model_data.E_comp{i}*model_data.C_comp{j}';
                end
            end
            
            estim.M3 = {};
            for i=1:length(model_data.C_comp)
                for j=1:length(model_data.C_comp)
                    for k=1:length(model_data.E_comp)
                        for l=1:length(model_data.E_comp)
                            estim.M4{i,j,k,l} = (W'*model_data.E_comp{k}*model_data.C_comp{i}')*(model_data.C_comp{j}*model_data.E_comp{l}'*W);
                        end
                    end
                end
            end

            estim.M5 = {};
            estim.M6 = {};
            for i = 1:length(model_data.A_comp)
                for j = 1:length(model_data.E_comp)
                    estim.M3{i,j} = W'*model_data.E_comp{j}*(model_data.A_comp{i})'*W;
                    estim.M6{i,j} = estim.M3{i,j}';
                end
                for j = 1:length(model_data.A_comp)
                    estim.M5{i,j} = W'*model_data.A_comp{i}*(model_data.A_comp{j}')*W;
                end
            end
            
             for i = 1:length(model_data.C_comp)
                 for j = 1:length(model_data.C_comp)
                     estim.M7{i,j} = model_data.C_comp{i}*model_data.C_comp{j}';
                 end
             end

            for i = 1:length(model_data.E_comp)
                for j = 1:length(model_data.E_comp)
                    estim.M8{i,j} = W'*model_data.E_comp{i}*(model_data.E_comp{j}')*W;
                end
            end

            
            
            if length(model_data.E_comp) > 1
                error('We cannot deal with more than one component in the matrices E,B,R,Q')
            end
            
            % Finally the last components for the error estimator:
            [~,nE] = evalc('svds(model_data.E_comp{1}, 1);');
            estim.normE = nE;
            nEinv = [];
            while isempty(nEinv)
              [~,nEinv] = evalc('1/svds(model_data.E_comp{1},1,0);');
            end
            estim.normEinv = nEinv;
            
            % Prepare the value for ||BR^{-1}B^T||
            estim.normBBT = {};
            estim.normBBT_use2norm = false;
            if length(model_data.B_comp) == 1
                estim.normBBT = norm(full(model_data.B_comp{1}));
                estim.normBBT_use2norm = true;
            else
                % We use the Frobenius norm as an upper bound for the
                % 2-norm
                for i = 1:length(model_data.B_comp)
                    for j = 1:length(model_data.B_comp)
                        for k = 1:length(model_data.B_comp)
                            for l = 1:length(model_data.B_comp)
                                estim.normBBT{i,j,k,l} = trace(model_data.B_comp{k}'*model_data.B_comp{l});
                            end
                        end
                    end
                end
            end

            estim.WtW = W'*W;
            
            % For the calculation of norm A:
            normA = {};
            for i = 1:length(model_data.A_comp)
                for j = 1:length(model_data.A_comp)
                    normA{i,j} = trace(model_data.A_comp{i}'*model_data.A_comp{j});
                end
            end
            estim.normA = normA;
            
            normB = {};
            for i = 1:length(model_data.B_comp)
                for j = 1:length(model_data.B_comp)
                    normB{i,j} = trace(model_data.B_comp{i}'*model_data.B_comp{j});
                end
            end
            estim.normB = normB;
            
            obj.estim = estim;
            if ~isempty(detailed_data.gamma)
                obj.gamma_function = detailed_data.gamma_calculation.gen_reduced_data(model, detailed_data, obj);
            end
            obj.gamma_mode = detailed_data.gamma_mode;
        end
        
    end
    
  
end