general/+Postprocess/+StochasticAssessment/Output.m

classdef Output < handle
  % result class for computations executed by an
  % Postprocess.StochasticAssessment.Assessment object.
  %
  % This class provides methods to transform and visualize the gathered data
  % for publication in LaTeX/Tikz documents.

  properties

    % values computed by Postprocess.StochasticAssessment.Assessment
    %
    % Usually this is a cell array of structure arrays. For cells correspond to
    % the @ref Assessment.rsamples "rsamples" and the structure array
    % corresponds to the test samples given by @ref Assessment.M_test "M_test".
    values;

    % cell array of all condition numbers of reduced simulaton system matrices.
    rd_conds;

    rtime_est = {};     % time for the computation of the error estimator. By default this is not set.

    % the object of type Assessment with which the values are generated
    %
    % @todo Maybe we should cache this, too...
    assessment;

    % maxmimum error values to be plotted. Above this values error display is
    % cropped.
    %
    % @default = 1;
    error_bound = 1;

    % maxmimum estimator values to be plotted. Above this values error display
    % is cropped.
    %
    % @default = 100;
    estimator_bound = 100;

    % define stability-region as error being smaller than sqrt(diffmax * area),
    % e.g. diffmax = 4, area = 2e-7

    % stability region
    %
    % If the averaged error exceeds this value, the visualization is cropped
    % here. This happens especially for small reduced basis sizes. (Default =
    % '1e-2')
    stab_limit = 1e-2;

  end

  methods

    function output = Output(values, assessment)
      % function output = Output(values, assessment)
      % Constructor
      %
      % Parameters:
      %  values:     the structure with the gathered data by the
      %              Postprocess.StochasticAssessment.Assessment object
      %  assessment: the generating object of type Postprocess.StochasticAssessment.Assessment

      output.values     = values;
      output.assessment = assessment;

    end

    function plot_landscape(this)
      %function plot_landscape(this)
      % plots a two dimensional error landscape in case of two dimensional
      % @ref Assessment.plot_fields "plot_fields"
      %
      % Used attributes:
      %  - #error_bound and
      %  - #stab_limit

      errs = this.get_slice('max_err');
      i = logical(errs>this.error_bound);
      errs(i) = this.error_bound;
      i = logical(isnan(errs));
      errs(i) = this.error_bound;
      reshape(errs, size(this.values));
      C = ones(size(errs));
      i = logical(errs<this.stab_limit);
      C(i) = 2;

      figure;
      if length(this.assessment.samples) == 2
        if isempty(find(C==1,1)); % i.e. all stable
          surf(this.assessment.samples{1},this.assessment.tsamples{2},errs');
        else % some not stable
          surf(this.assessment.samples{1},this.assessment.tsamples{2},errs',C');
        end;
        shading interp;
        % figure, pcolor(Ms, Ns,C);
        set(gca,'Zscale','log');
        xlabel(this.assessment.pf_descr{1});
        ylabel(this.assessment.pf_descr{2});
        zlabel(this.assessment.error_label);
      else
        plot(this.assessment.rsamples, errs');
        xlabel(this.assessment.pf_descr{1});
        ylabel(this.assessment.error_label);
        set(gca,'Yscale','log');
      end
      run_name = strrep(this.assessment.run_name,'_',' ');
      title([run_name,' L-infty([0,T],L2) error']);

      figure;
      if length(this.assessment.samples) == 2
        surf(this.assessment.samples{1}, this.assessment.samples{2}, real(this.get_slice('max_err_inds')));
        xlabel(this.assessment.pf_descr{1});
        ylabel(this.assessment.pf_descr{2});
        zlabel('time index');
      else
        plot(this.assessment.samples, this.get_slice('max_err_inds'));
        xlabel(this.assessment.pf_descr{1});
        zlabel('time index');
      end

      title('time index of maximum l2-error');
      %set(gca,'Zscale','log');
    end

    function print_table(this, datafile)
      % function print_table(this, datafile)
      % prints out a CSV table with gathered data.
      %
      % An example for a pgfplotstable graphic that can interpret the table output
      % is shown in the following code snippet: 
      % @code
      % \pgfplotstabletypeset[
      %   font={\footnotesize},
      %   columns={dim1,dim2,{time_avg},{err_max},{offtime_total}},
      %   columns/{dim1}/.style={
      %     column name={N},
      %     postproc cell content/.append code={%
      %       \ifnum\pgfplotstablerow=0
      %       \pgfkeyssetvalue{/pgfplots/table/@cell content}{$H={}$#2}%
      %       \fi
      %     }%
      %   },
      %   columns/{dim2}/.style={
      %     column name={M},
      %     postproc cell content/.append code={%
      %       \ifnum\pgfplotstablerow=0
      %       \pgfkeyssetvalue{/pgfplots/table/@cell content}{\ensuremath{-}}%
      %       \fi
      %     },
      %   },
      %   columns/{time_avg}/.style={
      %     column name={\o-run--time[s]}
      %   },
      %   columns/{err_max}/.style={
      %     column name={max. error}
      %   },
      %   columns/{offtime_total}/.style={
      %     column name={offline time[h]},
      %     divide by=3600
      %   },
      %   every even row/.style={
      %     before row={\rowcolor[gray]{0.9}}},
      %   every head row/.style={
      %     before row=\toprule,after row=\midrule},
      %   every last row/.style={
      %     after row=\bottomrule},
      %   display columns/2/.style={dec sep align},
      %   display columns/3/.style={zerofill,precision=2},
      %   empty cells with={\ensuremath{-}}
      % ]{table.dat};
      % @endcode
      %
      % Parameters:
      %  datafile: a string specifying the path and filename of the data table to be stored
      %
      % This function prints out a table with data for the average, min and max
      % values of
      %  - reduced simulation time,
      %  - reduced simulation reconstruction time,
      %  - reduced simulation error estimator computation time,
      %  - maximum error estimator and the
      %  - maximum "true" error
      % Furthermore, the offline time is printed out split up into
      %  - the total time spent in the offline phase,
      %  - the time spent for (pre-)computation of detailed simulations
      %  - the time spent in the EI-greedy algorithm
      %  - the time spent in the reduced basis generation algorithm
      ratime    = this.reduce_slice(@mean, 'rtime');
      rmtime    = this.reduce_slice(@min, 'rtime');
      rMtime    = this.reduce_slice(@max, 'rtime');
      rratime   = this.reduce_slice(@mean, 'rrtime');
      rrmtime   = this.reduce_slice(@min, 'rrtime');
      rrMtime   = this.reduce_slice(@max, 'rrtime');
      restatime = this.reduce_slice(@mean, 'rtime_est');
      restmtime = this.reduce_slice(@min, 'rtime_est');
      restMtime = this.reduce_slice(@max, 'rtime_est');
      etaa      = this.reduce_slice(@mean, 'max_Delta');
      etam      = this.reduce_slice(@min, 'max_Delta');
      etaM      = this.reduce_slice(@max, 'max_Delta');
      erra      = this.reduce_slice(@mean, 'max_err');
      errm      = this.reduce_slice(@min, 'max_err');
      errM      = this.reduce_slice(@max, 'max_err');

      dimension = ones(length(ratime),2);
      ratio = ones(length(ratime),1);
      simulation = cell(length(ratime),1);
      offt  = ones(length(ratime), 4);
      asm = this.assessment;
      rbdd = get_by_description(asm.detailed_data, 'rb');

      for comb = 1:length(asm.rsamples)
        pf_val = asm.rsamples(:,comb);
        for i = 1:length(pf_val)
          asm.update_rbmodel(asm.plot_fields{i}, asm.rsamples(i, comb));
        end
        N = asm.rbmodel.N;
        M = asm.rbmodel.M;
        simulation{comb} = 'Reduced';
        dimension(comb,1) = N;
        dimension(comb, 2) = M;
        if length(pf_val) == 1
          ratio(comb) = asm.rsamples(comb);
        else
          ratio(comb) = comb;
        end
        [offt(comb,1), offt(comb,2), offt(comb,3)] = offtime(asm.detailed_data, asm.rbmodel);
      end
      if offt(1,2) < 50
        offt(:,1) = 90;
      elseif offt(1,2) < 65
        offt(:,1) = 117;
      elseif offt(1,2) < 1300
        offt(:,1) = 887;
      else
        offt(:,1) = 1394;
      end
      offt(:,4) = sum(offt, 2);
      datime    = this.reduce_slice(@mean, 'dtime', [], 1);
      dmtime    = this.reduce_slice(@min, 'dtime', [], 1);
      dMtime    = this.reduce_slice(@max, 'dtime', [], 1);
      Hsize = asm.detailed_data.model_data.grid.nelements;
      ddimension = Hsize;
      dratio = -1;
      print_datatable(datafile,...
        {...
         'ratio', 'dim1', 'dim2', ...
         'time_avg', 'time_max', 'time_min', 'rec_time_avg', 'rec_time_min', 'rec_time_max', ...
         'est_time_avg', 'est_time_min', 'est_time_max', 'eta_avg', 'eta_min', 'eta_max', ...
         'err_avg', 'err_min', 'err_max', ...,
         'offtime_dsims', 'offtime_ei', 'offtime_rb', 'offtime_total' ...
        }, ...
        [...
         dratio, 0, ddimension, ...
         datime, dmtime, dMtime, -1, -1 , -1, ...
         -1, -1, -1, 0, 0, 0, ...
         0, 0, 0, ...
         0, 0, 0, 0; ...
         ratio(:), dimension, ...
         ratime(:), rmtime(:), rMtime(:), rratime(:), rrmtime(:), rrMtime(:), ...
         restatime(:), restmtime(:), restMtime(:), etaa(:), etam(:), etaM(:), ...
         erra(:), errm(:), errM(:), ...
         offt...
        ]');

    end

    function print_landscape(this, datafile, fieldname)
      % function print_landscape(this, datafile, fieldname)
      % equivalent to the plot_landscape() function put prints out a table
      % which can be interpreted e.g. by pgfplots/Tikz.
      %
      % An example for a pgfplots graphics that can interpret the table output
      % is shown in the following code snippet: 
      % @code
      % \begin{axis}[zmode=log,log basis z=10,axis on top,
      %              view={-50}{-30}, axis lines=box,
      %              xlabel={$N$}, ylabel={$M$},
      %              zlabel={$\eta_{N,M,M'}$}]
      %  \addplot3[surf]
      %  table[x=N, y=M, z=error] {error_landscape.dat};
      % @endcode
      %
      % Parameters:
      %  datafile:  a string specifying the path and filename of the data table to be stored
      %  fieldname: the fieldname in the #values attribute that shall be
      %             plotted @default 'max_err'
      if nargin <=2
        fieldname = 'max_err';
      end
      errs = this.reduce_slice(@max, fieldname);
      i = logical(errs>this.error_bound);
      errs(i) = this.error_bound;
      i = logical(isnan(errs));
      errs(i) = this.error_bound;

      rsamples    = this.assessment.rsamples;
      blocklength = length(this.assessment.samples{1});

      print_datatable(datafile, [this.assessment.plot_fields, {'error'}], ...
                                [rsamples; errs(:)'], blocklength);
%      C = ones(size(errs));
%      i = logical(errs<this.stab_limit);
%      C(i) = 2;
    end

    function print_3d_curve(this, datafile, fieldname)
      % function print_landscape(this, datafile, fieldname)
      % equivalent to the print_landscape() function put prints out a table
      % which can be interpreted for a 3d-curve plot e.g. by pgfplots/Tikz.
      %
      % This function can be used to display certain curves on the error
      % landscape generated by print_landscape().
      %
      % An example for a pgfplots graphics that can interpret the table output
      % is shown in the following code snippet: 
      % @code
      % \begin{axis}[zmode=log,log basis z=10,axis on top,
      %              view={-50}{-30}, axis lines=box,
      %              xlabel={$N$}, ylabel={$M$},
      %              zlabel={$\eta_{N,M,M'}$}]
      %  \addplot3[mark=0]
      %  table[x=N, y=M, z=error] {error_curve.dat};
      % @endcode
      %
      % Parameters:
      %  datafile:  a string specifying the path and filename of the data table to be stored
      %  fieldname: the fieldname in the #values attribute that shall be
      %             plotted @default 'max_err'
      if nargin <=2
        fieldname = 'max_err';
      end
      errs = this.reduce_slice(@max, fieldname);
      i = logical(errs>this.error_bound);
      errs(i) = this.error_bound;
      i = logical(isnan(errs));
      errs(i) = this.error_bound;

      Ns = ones(length(errs), 1);
      Ms = ones(length(errs), 1);
      rsamples = this.assessment.rsamples;
      for i=1:length(Ns)
        c = rsamples(i);
        this.assessment.update_rbmodel([], c);
        Ns(i) = this.assessment.rbmodel.N;
        Ms(i) = this.assessment.rbmodel.M;
      end

      print_datatable(datafile, {'c', 'N', 'M', 'error'}, ...
                                [rsamples(:), Ns, Ms, errs(:)]');
    end

    function slic = get_slice(this, fieldname, cell_range, struct_range)
      % function slic = get_slice(this, fieldname, cell_range, out_range)
      % helper function slicing a field in the #values attribute
      %
      % Parameters:
      %  fieldname:    name of the field in the #values attributed to be sliced
      %  cell_range:   cell indices to be extracted from the field.
      %                If empty all cell entries are returned
      %                @default []
      %  struct_range: Usually the #values attribute is a cell of structure arrays.
      %                arrays. This range defines a slice in the structure
      %                array. If empty all array entries are returned.
      %                @default []
      %
      % Return values:
      %  slic:  the sliced #values attribute
      if nargin <= 2 || isempty(cell_range)
        cell_range = 1:length(this.values(:));
      end
      if nargin <= 3 || isempty(struct_range)
        struct_range = 1:length(this.values{1}(:));
      end
      if isequal(fieldname, 'rtime_est')
        slic = cellfun(@(x) x(struct_range), this.rtime_est(cell_range), 'UniformOutput', false);
      else
        slic = cellfun(@(x) [x(struct_range).(fieldname)], this.values(cell_range), 'UniformOutput', false);
      end
      slic = cell2mat(slic');
    end

    function rslice = reduce_slice(this, red_func, fieldname, cell_range, struct_range, direction)
      % function rslice = reduce_slice(this, red_func, fieldname, cell_range, struct_range, direction)
      % extends the get_slice() function by also applying a reduction function
      % on each matrix entry.
      %
      % Parameters:
      %  fieldname:    name of the field in the #values attributed to be sliced
      %  red_func:     function pointer who is applied to each matrix entry in
      %                the slice, e.g. '@@max' or '@@mean'
      %  cell_range:   cell indices to be extracted from the field. If empty
      %                all cell entries are returned.
      %                @default []
      %  struct_range: Usually the #values attribute is a cell of structure arrays.
      %                arrays. This range defines a slice in the structure
      %                array. If empty all array entries are returned.
      %                @default []
      %  direction:    integer specifying in which direction a matrix shall be
      %                reduced. In case it is equal to '2', the matrix is
      %                transposed before 'red_func' is applied.
      %                @default 1
      %
      % Return values:
      %  rslice:  the sliced and reduced #values attribute
      if nargin <=3
        cell_range = [];
      end
      if nargin <= 4
        struct_range = [];
      end
      if nargin <= 5
        direction = 1;
      end
      slic = get_slice(this, fieldname, cell_range, struct_range);
      if direction == 2
        slic   = slic';
        rslice = zeros(1,size(slic,1));
      else
        rslice = zeros(size(slic,1),1);
      end
      for i = 1:size(slic,1)
        rslice(i) = red_func(slic(i,:));
      end
    end


    function merge(this, other)
      % function merge(this, other)
      % helper function merging this Output object with another one.
      %
      % The following rules apply during merging:
      %  -# If 'values.max_Delta' is all 'Inf', copy 'max_Delta' and
      %     'max_Delta_ind' from 'other'.
      %  -# If 'values.max_err' is all 'Inf', copy 'max_err' and
      %     'max_err_ind' from 'other'.
      %  -# If one of 'rconds', 'dconds' or 'rtime_est' is empty, copy it
      %     from 'other'
      %  -# If one of 'rtime', 'rrtime' or 'dtime' is set to all-'NaN', copy it
      %     from 'other'
      %  -# If 'this' and 'other' have valid 'rtime' (i.e. not-'Nan') and only
      %     one of both has computed estimators, compute 'rtime_est' by
      %     subtracting one from the other.
      %
      % Parameters:
      %  other: another object of type Postprocess.StochasticAssessment.Output to be
      %         merged into this one.
      cpfields = {};
      rtime_exists = false;
      % check what we ar missing:
      if all([this.values{end}(:).max_Delta] == Inf)
        cpfields = [ cpfields, {'max_Delta', 'max_Delta_ind'}];
      end
      if all([this.values{end}(:).max_err] == Inf)
        cpfields = [ cpfields, {'max_err', 'max_err_ind'}];
      end
      if isempty([this.values{end}(:).rconds])
        cpfields = [ cpfields, {'rconds'}];
      end
      if isempty([this.values{end}(:).dconds])
        cpfields = [ cpfields, {'dconds'}];
      end
      if all(isnan([this.values{end}(:).rtime]))
        cpfields = [ cpfields, {'rtime'}];
      else
        rtime_exists = true;
      end
      if all(isnan([this.values{end}(:).rrtime]))
        cpfields = [ cpfields, {'rrtime'}];
      end
      if all(isnan([this.values{end}(:).dtime]))
        cpfields = [ cpfields, {'dtime'}];
      end

      for i = 1:length(cpfields)
        for j = 1:length(this.values(:))
          [this.values{j}(:).(cpfields{i})] = other.values{j}(:).(cpfields{i});
        end
      end

      if rtime_exists && any(~isnan([other.values{end}(:).rtime]))
        if this.assessment.compute_estimates ~= other.assessment.compute_estimates
          this.rtime_est = cell(size(this.values));
          for i = 1:length(this.values(:))
            this.rtime_est{i} = abs([this.values{i}(:).rtime] - [other.values{i}(:).rtime]);
          end
        end
      end

      if isempty(this.rtime_est)
        this.rtime_est = other.rtime_est;
      end

      if isempty(this.rd_conds)
        this.rd_conds = other.rd_conds;
      end
    end
  end
end