ISeparableOperator.m

classdef ISeparableOperator < handle

  properties (Constant, Abstract)
    id;

    ret_vars;

    arg_vars;

  end

  methods (Abstract)

    res = apply(this, descr, model_data, sim_data);

    [LL, bb]               = full_matrix(this, descr, model_data);

    [LL_comps, bb_comps]   = components(this, descr, model_data);

    [LL_coeffs, bb_coeffs] = coefficients(this, descr);

    A = inner_product_matrix(this, descr, model_data);
  end

  methods
    function [LL, bb] = gram_projection(this, model, detailed_data)
      first_dd = get_by_description(detailed_data, this.arg_vars{1});
      model_data = first_dd.model_data;

      [L,bb] = full_matrix(this, descr, model_data);

      RB = cell(1, length(this.arg_vars));
      for i = 1:length(this.arg_vars)
        dd_i = get_by_description(detailed_data, this.arg_vars{i});
        RB{i} = dd_i.RB;
      end
      RB = cat(1, RB{:});

      LL = L * RB;

      if ~isempty(this.ret_vars)
        A = inner_product_matrix(this, model, model_data);

        RB = cell(1, length(this.ret_vars));
        for i = 1:length(this.ret_vars)
          dd_i = get_by_description(detailed_data, this.ret_vars{i});
          RB{i} = dd_i.RB;
        end
        RB = cat(1, RB{:});
        LL = RB' * A * LL;
        bb = bb * A * RB;
      end
    end

    function [LL_comps, bb_comps] = gram_project_components(this, rmodel, detailed_data)
      first_dd = get_by_description(detailed_data, this.arg_vars{1});
      model_data = first_dd.model_data;

      RB = cell(1, length(this.arg_vars));
      for i = 1:length(this.arg_vars)
        dd_i = get_by_description(detailed_data, this.arg_vars{i});
        N_i = get_by_description(rmodel.N, this.arg_vars{i});
        RB{i} = dd_i.RB(:,1:N_i);
      end
      RB = cat(1, RB{:});

      [L_comps, bb_comps] = components(this, rmodel.descr, model_data);
      LL_comps = cellfun(@(X) X * RB, L_comps, 'UniformOutput', false);

      if ~isempty(this.ret_vars)
        A = inner_product_matrix(this, rmodel, model_data);

        RB = cell(1, length(this.ret_vars));
        for i = 1:length(this.ret_vars)
          dd_i = get_by_description(detailed_data, this.ret_vars{i});
          N_i = get_by_description(rmodel.N, this.ret_vars{i});
          RB{i} = dd_i.RB(:,1:N_i);
        end
        RB = cat(1, RB{:});

        LL_comps = cellfun(@(X) RB' * A * X, ...
                      LL_comps, 'UniformOutput', false);
        bb_comps = cellfun(@(X) RB' * A * X, ...
                      bb_comps, 'UniformOutput', false);
      end

    end
  end
end