ReducedData.m

classdef ReducedData < IReducedData
  % Reduced data implementation for linear stationary problems with finite
  % element discretizations.
  %
  % @todo Describe the matrix components, or add a link to a paper...


  properties
    Q_A;            % number of components in matrix `A`
    Q_f;            % number of components in vector `f`
    Q_l;            % number of components in vector `l`
    AN_comp;        % reduced components `A^q`
    fN_comp;        % reduced components `f^q`
    lN_comp;        % reduced components `l^q`
    Gff;            % estimator matrix components `G_{ff}`
    Gfa;            % estimator matrix components `G_{fA}`
    Gaa;            % estimator matrix components `G_{AA}`
    G;              % estimator matrix components `G`
  end

  properties (SetAccess = private, Dependent)
    % number of reduced basis vectors stored in this data node.
    N;
  end

  methods

    function res = ReducedData(rmodel, detailed_data)
      % computes the reduced basis vectors and matrices @b or implements a copy
      % constructor.
      %
      % Parameters:
      %   rmodel:        of type LinStat.ReducedModel
      %   detailed_data: of type LinStat.DetailedData

      if nargin == 0
        error('LinStat.ReducedData constructor needs an argument');
      elseif nargin == 2 && isa(rmodel, 'IReducedModel') && isa(detailed_data, 'IDetailedData')
        fill_fields(res, rmodel, detailed_data);
      elseif nargin == 1 && isa(rmodel, 'IReducedData')
        fns = properties(rmodel);
        fns = setdiff(fns,'N');
        for i = 1:length(fns)
          res.(fns{i}) = rmodel.(fns{i});
        end
        %        res.handle_to_last_subset = arg1.handle_to_last_subset;
      else
        error('Did not find constructor for your arguments');
      end
    end

    function N = get.N(this)
      if ~isempty(this.AN_comp{1})
        N = size(this.AN_comp{1}, 2);
      else
        N = 0;
      end
    end

    function reduced_data_subset = extract_reduced_data_subset(this, rmodel)
      % function reduced_data_subset = extract_reduced_data_subset(this, rmodel)
      % @copybrief IReducedModel.extract_reduced_data_subset()
      %
      % Parameters:
      %  rmodel: reduced model of type LinStat.ReducedModel indicating the
      %          number of reduced basis vectors to be used in the
      %          'reduced_data_subset'
      %
      % Return values:
      %  reduced_data_subset: object of type LinStat.ReducedData which is a
      %                       deep copy of this data node but with a number of
      %                       reduced basis vectors indicated by 'rmodel.N'. If
      %                       the number did not change, only a handle copy is
      %                       returned.

      if this.N ~= rmodel.N
        % extract correct N-sized submatrices and subvectors from reduced_data
        if rmodel.N > this.N
          error('N too large for current size of reduced basis!');
        end

        % create a copy
        reduced_data_subset = LinStat.ReducedData(this);

        N = rmodel.N;

        reduced_data_subset.AN_comp   = ...
          subblock_sequence(this.AN_comp,1:N,1:N);
        reduced_data_subset.fN_comp = subblock_sequence(this.fN_comp,1:N);
        reduced_data_subset.lN_comp = subblock_sequence(this.lN_comp,1:N);

        % inner product matrix of Riesz-representers of residual components:
        reduced_data_subset.Gff = this.Gff;
        rQ_f = length(reduced_data_subset.fN_comp);
        Q_a  = length(reduced_data_subset.AN_comp);
        reduced_data_subset.G = this.G(1:(rQ_f+N*Q_a),1:(rQ_f+N*Q_a));
        reduced_data_subset.Q_f = rQ_f;
      else
        reduced_data_subset = this;
      end
    end

    function conds = get_conds(this)

      disp('get_conds needs to be implemented for LinStat.ReducedData');
      conds = [];
    end
  end

  methods(Access=private)

    function fill_fields(reduced_data, dmodel, detailed_data)
      % function fill_fields(reduced_data, dmodel, detailed_data)
      % actually fills the reduced data fields

      model = dmodel.descr;

      model.decomp_mode = 1; % == components

      [A_comp,f_comp] = model.operators(model,detailed_data.model_data);

      Q_A = length(A_comp);
      reduced_data.Q_A = Q_A;
      Q_f = length(f_comp);
      reduced_data.Q_f = Q_f;

      reduced_data.AN_comp = cell(1,length(A_comp));
      for q = 1:Q_A
        reduced_data.AN_comp{q}  = ...
          detailed_data.RB'*A_comp{q}*detailed_data.RB;
      end;

      reduced_data.fN_comp = cell(1,length(f_comp));
      for q = 1:Q_f
        reduced_data.fN_comp{q} = ...
          detailed_data.RB' * f_comp{q};
      end;

      if model.compute_output_functional
        % assumption: nonparametic output functional, then simple RB
        % evaluation possible  
        l_comp = ...
          model.operators_output(model,detailed_data.model_data);
        Q_l = length(l_comp);
        reduced_data.lN_comp = cell(1,Q_l);
        for q = 1:Q_l
          reduced_data.lN_comp{q} = detailed_data.RB' * l_comp{q}; 
        end;
        reduced_data.Q_l = Q_l;
      end;

      % plus error estimation quantities

      % G = (v_r^q, v_r^q) = v_r^q' * K * v_r^q
      % with {v_r^q}_q = (v_f^q, v_a^qq')_{qq'}
      % G = [Gff, Gfa; Gfa', Gaa];
      %
      % (v_f^q,v_f^q') = v_f^q' * K * v_f^q

      % matrices of coefficient vectors of Riesz-representers:
      % K * v_f^q = f^q      (coefficient vector equations)
      K = model.get_inner_product_matrix(detailed_data.model_data);
      % search solution in H10, i.e. set dirichlet DOFs

      N = get_rb_size(detailed_data);

      v_f = zeros(size(detailed_data.RB,1),Q_f);
      K_v_f = zeros(size(detailed_data.RB,1),Q_f);
      for q = 1:Q_f
        K_v_f(:,q) = f_comp{q};
        v_f(:,q) = K \ f_comp{q};
      end;

      v_a = cell(N,1);
      K_v_a = cell(N,1);
      for n = 1:N
        K_v_a{n} = zeros(size(K,1),Q_A);
        v_a{n} = zeros(size(K,1),Q_A);
        for q = 1:Q_A
          K_v_a{n}(:,q) = (A_comp{q}*detailed_data.RB(:,n));
          v_a{n}(:,q) = K \ K_v_a{n}(:,q);
        end;
      end;

      % compute matrix G components = [G_ff, G_fa{n}; G_fa{n}', G_aa{n}];
      reduced_data.Gff = v_f' * K_v_f;
      reduced_data.Gfa = cell(1,N);
      for n = 1:N
        reduced_data.Gfa{n} = v_f' * K_v_a{n};
      end;
      reduced_data.Gaa = cell(N,N);
      for n1 = 1:N
        for n2 = 1:N
          reduced_data.Gaa{n1,n2} = v_a{n1}' * K_v_a{n2};
        end;
      end;

      % finally assemble G = [G_ff, G_fa{n}; G_fa{n}', G_aa{n}];
      Q_r = Q_f + N * Q_A;
      G = zeros(N,N);
      G(1:Q_f,1:Q_f) = reduced_data.Gff;
      for n = 1:N
        G(1:Q_f,Q_f+(n-1)*Q_A +(1:Q_A)) = ...
          reduced_data.Gfa{n};
        G(Q_f+(n-1)*Q_A +(1:Q_A),1:Q_f) = ...
          reduced_data.Gfa{n}';
      end;
      for n1 = 1:N
        for n2 = 1:N
          G(Q_f+(n1-1)*Q_A+(1:Q_A),Q_f+(n2-1)*Q_A +(1:Q_A)) = ...
            reduced_data.Gaa{n1,n2};
        end;
      end;
      reduced_data.G = G;

    end
  end
end