general/+Postprocess/+StochasticAssessment/Assessment.m

classdef Assessment < handle
  % a class used to compute reduced several reduced simulations over a huge
  % parameter sample extracting useful information
  %
  % The class generates one or two dimensional cell arrays of these data fields
  % storing information for a variation of some attributes of the reduced
  % model, which can be freely chosen by specification of #plot_fields .

  properties(Access=public)
    % is a vector of field names with at most 2 elements over which the error
    % landscape is computed or a cell with empty entries.
    %
    % A reasonable choice would be '{ \'N\', \'M\' }'. The error landscape is
    % plotted over a variation of these fields specified by samples.
    % In case of an empty cell entry, 'N' and 'M' are coupled by the fixed
    % ratio #M_by_N_ratio and samples should be a vector of coupling constant
    % between 0 and 1. See also: update_rmodel()
    %
    % @default = '{[]}'
    plot_fields = {[]};

    % is a cell array of scalar row vectors. The scalar values are 'plot_field'
    % values for which the error is computed and tested.
    %
    % @default = '{[0.1:0.1:1]}'
    samples = {(0.1:0.1:1)};

    % an parameter sampling object of type ParameterSampling.Interface
    %
    % @default is a random sampling with 10 elements and seed 654321;
    %
    % See also: init_random_sampling()
    M_test;

    % is a cell array of description texts for the plot fields.  If this field
    % is not set it is set to 'plot_fields'.
    plot_field_descr = {};


    % character string specifying this test run. The string should be unique
    % for every parameter combination.
    %
    % @default rmodel.name + '_stoachastic_assessment';
    run_name;

    % boolean value determining whether error estimates shall be computed
    compute_estimates  = false;

    % boolean value determining whether the "true" error `\|u_h(\mu) -
    % u_{\text{red}}(\mu)\|_{{\cal W}_h}` shall be computed
    compute_errors     = false;

    % boolean value determining whether the condition numbers of the system
    % matrices shall be computed
    compute_conditions = false;

    % boolean value indicating whether the refinement steps of the reduced
    % basis generation by Greedy.TrainingSetAdaptation shall be taken into account
    %
    % This option is only useful if
    %   -# #M_test is set to the training sample of the parameter space
    %     generation and
    %   -# the #plot_fields include the basis dimensions
    % In this case the parameter sampling used during basis generation for the
    % given basis dimension is selected.
    follow_refinement_steps = false;

    % if M and N are couple this value specifies the constant for the coupling.
    %
    % A value of zero means `M_{\max} / N_{\max}`
    M_by_N_ratio = 1;

    % the underlying reduced model of type IReducedModel
    rmodel;

    % an object of type CacheableObject storing the reduced model and the
    % detailed data
    cache_object;

  end

  properties(Dependent)
    % the "real" collection of vectors used for the #plot_fields variation
    %
    % See also: #samples
    rsamples;

    % the underlying detailed data of type Greedy.DataTree.Detailed.INode
    %
    % @todo make this cached...
    detailed_data;
  end

  methods

    function dd = get.detailed_data(this)
      dd = this.cache_object.obj.detailed_data;
    end

    function rsamples = get.rsamples(this)
      if length(this.samples) == 2
        rsamples = combvec(this.samples{1},this.samples{2});
      else
        rsamples = this.samples{1};
      end
    end

    function sta = Assessment(matfile, rmodel, mu_set_size, seed)
      % function sta = Assessment(matfile, rmodel, mu_set_size, seed)
      % constructor
      %
      % Parameters:
      %   matfile:     name of result file where a IReducedModel object and a
      %                Greedy.DataTree.Detailed.INode object must be stored.
      %   mu_set_size: Optional argument determining the number of random
      %                parameters in the validation sample. (Default = 10)
      %   seed:        a random seed used for initialization of
      %                ParameterSampling.Random object for the #M_test parameter
      %                sample set. (Default = 654321)

      if nargin == 1 && isa(matfile, 'Postprocess.StochasticAssessment.Assessment')
        cp = matfile;
        m = metaclass(cp);
        fns = m.Properties(cellfun(@(x) ~isequal(x.SetAccess, 'none') && ~x.Dependent, m.Properties));
        fns = cellfun(@(x) x.Name, fns, 'UniformOutput', false);
        fns = setdiff(fns, {'cache_object'});
        for i = 1:length(fns);
          sta.(fns{i}) = cp.(fns{i});
        end
        sta.cache_object = CacheableObject(cp.cache_object.matfile, 'detailed_data');
        sta.rmodel       = copy(cp.rmodel);
      else
        sta.rmodel       = copy(rmodel);
        sta.cache_object = CacheableObject(matfile, 'detailed_data');
        if nargin <= 2
          mu_set_size = 10;
        end
        if nargin <= 3
          seed = 654321;
        end
        init_random_sampling(sta, mu_set_size, seed);

        if isempty(sta.plot_field_descr)
          sta.plot_field_descr = sta.plot_fields;
        else
          assert(length(sta.plot_field_descr) == length(sta.plot_fields));
        end

        sta.run_name = [rmodel.detailed_model.name, '_stochastic_assessment'];
      end
    end

    function output=compute(this)
      % function compute(this)
      % stochastic estimation of error between reduced and detailed simulation over a
      % test sample of `\mu` vectors
      %
      % This function stochastically estimates the error between reduced and detailed
      % simulations for given `\mu`-vectors and various reduced and collateral basis
      % sizes. The results are visualized in a surface plot for problems with CRB.
      %
      % If required by the users, averaged time measurements for the reduced and the
      % detailed simulations are computed, too.
      %
      % Return values:
      %   output: an Output object

      %if isequal(params.mode, 'error_ei')
      %  model.detailed_simulation = model.detailed_ei_simulation;
      %end

      %if isequal(params.mode, 'error_ei_rb_proj')
      %  model.detailed_simulation = model.detailed_ei_rb_proj_simulation;
      %end

      reduced_data = gen_reduced_data(this.rmodel, this.detailed_data);

      if this.compute_estimates
        assert(this.rmodel.Mstrich > 0);
        assert(this.rmodel.enable_error_estimator);
      else
        assert(~this.rmodel.enable_error_estimator);
        assert(this.rmodel.Mstrich == 0);
      end

      basis_gen = this.detailed_data.bg_algorithm;
      while ~isa(basis_gen, 'Greedy.Algorithm')
        switch(class(basis_gen))
          case 'Greedy.Combined'
            basis_gen = get_rb_basis_generator(basis_gen);
          case 'Greedy.TrainingSetAdaptation'
            basis_gen = basis_gen.child;
          otherwise
            error('nonono');
        end
      end
      generator = basis_gen.detailed_gen;

      if ~generator.is_valid
        generator = SnapshotsGenerator.Trajectories(this.rmodel.detailed_model);
      end

      % default value (only need for mode == 'check' when no 'mu_set_size' is given)

      savefile = [this.run_name];
      tmpfile  = [savefile,'_tmp.mat'];
      savefile = [savefile,'.mat'];

      rps_size = length(this.rsamples);

      complete_out = cell(1, rps_size);

      if this.compute_conditions
        rd_conds = cell(1, rps_size);
      else
        rd_conds = {};
      end

      for comb = 1:rps_size

        pf_val=this.rsamples(:,comb);

        for i=1:length(pf_val)
          update_rmodel(this, this.plot_fields{i}, pf_val(i));
        end
        temp_rd = extract_reduced_data_subset(this.rmodel, reduced_data);

        if this.compute_conditions
          rd_conds{comb} = get_conds(temp_rd);
        end

        disp(['Processing combination ', num2str(comb), '/', num2str(rps_size)]);

        M = this.get_test_space_sample;
        generator.prepare(this.rmodel.detailed_model, this.detailed_data.model_data, M);
        mu_size = size(M,1);

        part_out = struct('max_Delta'    , cell(1, mu_size), ...
                          'max_Delta_ind', cell(1, mu_size), ...
                          'max_err'      , cell(1, mu_size), ...
                          'max_err_ind'  , cell(1, mu_size), ...
                          'rtime'        , cell(1, mu_size), ...
                          'rrtime'       , cell(1, mu_size), ...
                          'dtime'        , cell(1, mu_size), ...
                          'rconds'       , cell(1, mu_size), ...
                          'dconds'       , cell(1, mu_size));
        %for mu_ind = 1:mu_size
        parfor mu_ind = 1:mu_size
          tthis        = this;
          trm          = this.rmodel;
          tdm          = trm.detailed_model;
          trd          = temp_rd;
          tgen         = generator;

          disp(['processing parameter vector ',num2str(mu_ind),...
            '/',num2str(size(M,1))]);

          set_mu(trm, M(mu_ind,:));
          set_mu(tdm, M(mu_ind,:));

          compute_out      = compute_impl(tthis, tgen, trm, trd);

          part_out(mu_ind) = compute_out;
        end

        complete_out{comb} = part_out;
%        disp(['saving temporary workspace in ', tmpfile]);
%        rmodel       = this.rmodel;
%        detailed_data = this.detailed_data;
%        save(fullfile(rbmatlabtemp,tmpfile), 'complete_out', 'part_out', 'this', 'rmodel', 'detailed_data', 'reduced_data');

      end

      if length(this.samples) == 2
        ns1 = length(this.samples{1});
        ns2 = rps_size / ns1;
        complete_out = reshape(complete_out, ns1, ns2);
        if this.compute_conditions
          rd_conds = reshape(rd_conds, ns1, ns2);
        end
      end

      output = Postprocess.StochasticAssessment.Output(complete_out, this);
      output.rd_conds = rd_conds;

      save(fullfile(rbmatlabresult, savefile));

    end

    function init_random_sampling(this, mu_set_size, seed)
      % function init_random_sampling(this, mu_set_size, seed)
      % generates the initial random sample set #M_test
      %
      % Parameters:
      %   mu_set_size: @copydoc ParameterSampling.Random.nparameters
      %   seed       : @copydoc ParameterSampling.Random.seed

      detailed_model = this.rmodel.detailed_model;
      this.M_test = ParameterSampling.Random(mu_set_size, seed);
      if this.M_test.init_required()
        init_sample(this.M_test, detailed_model);
      end
    end

    function other = copy(this)
      % function other = copy(this)
      % deep copy of this object
      %
      % Return values:
      %   other:  copied Assessment object
      other = Postprocess.StochasticAssessment.Assessment(this);
    end

    function equal_distribution_samples(this, min, max, sample_size)
      % function equal_distribution_samples(min, max, sample_size)
      % adapts the #samples attribute such that they are equally distributed in
      % a given range.
      %
      % Parameters:
      %   min:         is a vector of minimum values of 'plot_fields' variables.
      %   max:         is a vector for maximum values of 'plot_fields'
      %                variables.
      %   sample_size: specifies the number of sample values between 'min' and
      %                'max' for which the error is computed and plotted. The
      %                default value is 'max-min'.

      if nargin == 2
        sample_size = max-min;
      end


      assert(all(max > min));

      this.samples = cell(1,length(max));
      for i = 1:length(max)
        this.samples{i} = min(i) + ...
          round((0:sample_size(i)-1) .* ...
          (max(i)-min(i))./(sample_size(i)-1));
      end
    end

    function update_rmodel(this, plot_field, new_value)
      % function update_rmodel(this, plot_field, new_value)
      % updates the reduced model, i.e. applies the changes defined by the
      % #plot_fields attribute
      %
      % Parameters:
      %   plot_field:  an entry in #plot_fields, if it is empty, a coupling of
      %                'N' and 'M' is assured.
      %   new_value:   the new value it should be set to.

      if ~isempty(plot_field)
        this.rmodel.(plot_field) = new_value;
      else
        couple_N_and_M_by_c(this, new_value);
      end

      rmodel = this.rmodel;
      disp(['N = ', num2str(rmodel.N), ' M = ', num2str(rmodel.M), ' Mstrich = ', num2str(rmodel.Mstrich)]);
    end


  end


  methods(Static)

    function [rb_sim_data, tictoc] = rb_simulation_tictoc(varargin)
      % function [rb_sim_data, tictoc] = rb_simulation_tictoc(varargin)
      % a wrapper around the IReducedModel.rb_simulation() measuring the
      % execution time
      %
      % Return values:
      %   rb_sim_data:  the reduced simulation data
      %   tictoc:       execution time
      [rb_sim_data, tictoc] = tictoc_wrapper(@rb_simulation, varargin{:});
    end


  end

  methods(Access=private)

    function output =  compute_impl(this, generator, rmodel, reduced_data)

      fprintf('[');

      model_data = this.detailed_data.model_data;

      max_Delta     = NaN;
      max_Delta_ind = 0;
      max_err       = NaN;
      max_err_ind   = 0;
      dtime         = NaN;
      rtime         = NaN;
      rrtime        = NaN;
      rconds        = [];
      dconds        = [];
      rmodel.compute_conditions = this.compute_conditions;
      try
        [rb_sim_data,tictoc] = this.rb_simulation_tictoc(rmodel, reduced_data);

        fprintf('N');
        if this.compute_estimates
          Delta = rb_sim_data.Delta;
          [max_Delta, max_Delta_ind] = max(Delta);
        else
          max_Delta     = Inf;
        end

        if this.compute_conditions
          rconds = rb_sim_data.conds;
        else
          % condition computations may make time measurements bogus.
          rtime                = tictoc;
        end

        if this.compute_errors
          dmodel = rmodel.detailed_model;
          if this.compute_conditions
            fields = {'conds'};
            dmodel.compute_conditions = true;
          else
            fields = [];
          end
          [U_H, dtic_toc, opts] = generate(generator, dmodel, model_data, fields);
          fprintf('H');
          if this.compute_conditions
            dconds = opts.conds;
          else
            % condition computations may make time measurements bogus.
            dtime = dtic_toc;
          end

          tic;
          rb_sim_data = rb_reconstruction(rmodel, this.detailed_data, rb_sim_data);
          rrtime = toc;
          fprintf('R');

          errors      = dmodel.l2_error_sequence_algorithm(U_H, rb_sim_data.U, ...
            model_data);
          [max_err,max_err_ind] = max(errors);
        else
          max_err     = Inf;
        end

      catch ME
        warning(['catched an error: ', ME.message]);
        disp('setting error to NaN');
      end
      fprintf(']');
      output.max_Delta     = max_Delta;
      output.max_Delta_ind = max_Delta_ind;
      output.max_err       = max_err;
      output.max_err_ind   = max_err_ind;
      output.rtime         = rtime;
      output.rrtime        = rrtime;
      output.dtime         = dtime;
      output.rconds        = rconds;
      output.dconds        = dconds;
    end


    function couple_N_and_M_by_c(this, c)
      % function model = couple_N_and_M_by_c(model, c)
      % modifies the reduced basis size fields of 'model' by a single variable.
      %
      % Parameters:
      %  c:  sets the reduced basis sizes to
      %      `(N,M) = c (N_{\mbox{max}}, c_{MbyN} N_{\mbox{max}})`. For time-adaptive
      %      schemes with several collateral reduced basis spaces, the following
      %      formula is applied: `(N,M^k) = c (N_{\mbox{max}}), c_{MbyN}
      %      \frac{M^k_{\mbox{max}}}{\mbox{max}_{k=0,...,K} M^k_{\mbox{max}}}
      %      N_{\mbox{max}}.`

      Nmax      = get_rb_size(this.detailed_data, this.rmodel);
      Mmax      = get_ei_size(this.detailed_data, this.rmodel);
      maxMmax   = max(Mmax)-this.rmodel.Mstrich;
      this.rmodel.N = max(round(c*Nmax),1);
      if this.M_by_N_ratio == 0
        M_by_N_r = maxMmax / Nmax;
      elseif this.M_by_N_ratio == -1
        leaf_dd = get_active_leaf(this.detailed_data, this.rmodel);
        assert(isa(leaf_dd, 'Greedy.DataTree.Detailed.PODEI'));
        ext_ind        = this.rmodel.N;
        ext_counts = get_field(leaf_dd, 'ei_ext_count');
        ei_ext_steps = get_field(leaf_dd, 'ei_ext_steps');
        rb_real_ext_steps = get_field(leaf_dd, 'rb_real_ext_steps');
        if isempty(ext_counts)
          n_ei_exts  = length(ei_ext_steps);
          ext_counts = ones(size(ei_ext_steps))*maxMmax/n_ei_exts;
        end
        Mcomp = round(sum((ei_ext_steps <= rb_real_ext_steps(ext_ind)).*ext_counts));
        this.rmodel.M = min(maxMmax, max(Mcomp, 1));
        return
      else
        M_by_N_r = this.M_by_N_ratio;
      end
      this.rmodel.M = arrayfun(@(x) min(x, max(round(c*x/maxMmax*M_by_N_r*Nmax),1)), maxMmax);
    end

    function M = get_test_space_sample(this)
      if this.follow_refinement_steps
        info = active_info(this.detailed_data);
        ref_level = get_ref_level(info, this.rmodel);
        M = get_sample_at_ref_level(this.M_test, ref_level);
      else
        M = this.M_test.sample;
      end
    end

  end
end