DetailedData.m

classdef DetailedData < Greedy.User.IDetailedData
  % class generating the reduced basis space for the LinEvol problem with a
  % Greedy algorithm.
  %
  % Parameters given by the constructor argument 'bg_descr'
  % can control the
  %   - trainings sample (#train_sample_mode, #train_num_intervals, #train_seed)
  %   - trainings sample refinement (#refinement_mode, #max_refinement_level,
  %                      #refinement_theta, #val_size, #val_seed, #stop_max_val_train_ratio)
  %   - general EI extension controls (#noof_ei_extensions) and
  %   - general greedy controls (#stop_Nmax, #stop_Mmax, #stop_timeout,
  %                              #stop_epsilon, #stop_max_val_train_ratio)

  properties (Dependent, SetAccess=private)
    % a basis generation algorithm of type Greedy.Interface used to generated
    % the reduced basis space.
    bg_algorithm;
  end

  properties (Access=private)
    rb_generator;

    ei_generators = {};
  end

  properties
    % flag indicating whether the parameter sampling for the training set
    % `M_{\text{train}}` type.
    %
    % Possible values:
    %   - 'uniform': (default) for a ParameterSampling.Uniform
    %   - 'random':  for a ParameterSampling.Random
    %   - 'single':  constructs a basis from a single trajectory of solution
    %                snapshots (cog of parameter sampling set).
    train_sample_mode        = 'uniform';

    % number of intervals in the ParameterSampling.Uniform respectively number
    % of vectors in the ParameterSampling.Random for the training set
    % `M_{\text{train}}` of the parameter space `{\cal M}`.
    train_num_intervals      = 3;

    % If 'train_sample_mode == random',  a ParameterSampling.Random is generated
    % for the training set `M_{\text{val}}` with size #train_num_intervals and
    % seed 'train_seed'.
    train_seed               = 4321;

    % sets the Greedy.TrainingSetAdaptation.refinement_mode option
    %
    % Possible values are:
    %  - 'uniform': the parameter sampling is refined uniformly
    %  - 'adaptive': the parameter sampling is refined adaptively, such that
    %                only regions of the parameter space with high estimated
    %                errors are refined.
    %  - 'none':     no Greedy.TrainingSetAdaptation object is created
    refinement_mode          = 'adaptive';

    % sets the Greedy.TrainingSetAdaptation.max_refinement_level option
    %
    % @copydetails Greedy.TrainingSetAdaptation.max_refinement_level
    max_refinement_level     = 5;

    % sets the Greedy.TrainingSetAdaptation.refinement_theta option
    %
    % @copydetails Greedy.TrainingSetAdaptation.refinement_theta
    refinement_theta         = 0.2;

    % If this value is non-zero, a ParameterSampling.Random is generated for the
    % validation set `M_{\text{val}}` with size #val_size and seed #val_seed.
    val_size                 = 10;

    % If the value of 'val_size' is non-zero,  a ParameterSampling.Random is
    % generated for the validation set `M_{\text{val}}` with size #val_size
    % and seed 'val_seed'.
    val_seed                 = 1234;

    % sets the Greedy.Plugin.EI.stop_Mmax option
    stop_Mmax                 = 200;

    % sets the Greedy.Plugin.PODEI.stop_Nmax option
    stop_Nmax                 = 50;

    % sets the Greedy.Algorithm.stop_epsilon option
    stop_epsilon              = 1e-5;

    % sets the Greedy.Algorithm.stop_max_val_train_ratio option
    stop_max_val_train_ratio  = 1.2;

    % sets the Greedy.Algorithm.stop_timeout option
    stop_timeout              = 12*60*60;

    % sets the number of ei extensions
    %
    % See Greedy.Plugin.EI.noof_ei_extensions for details.
    noof_ei_extensions       = 1;

    pca                      = false;
  end

  methods
    function dd = DetailedData(rmodel, model_data, checkpoint)
      % function dd = DetailedData(rmodel, model_data)
      % constructor constructing the reduced basis spaces and storing it in a
      % DataTree.
      %
      % Parameters:
      %   rmodel: of type LinEvol.ReducedModel
      
      if nargin == 2
        checkpoint = [];
      end

      dd              = dd@Greedy.User.IDetailedData(rmodel.bg_descr, model_data);

      dd.rb_generator = SnapshotsGenerator.Trajectories(rmodel.detailed_model);

      dmodel = rmodel.detailed_model;
      descr = dmodel.descr;

      Lsat_ei_gen = SnapshotsGenerator.SpaceOpEvals(dmodel, 'L_saturation', ...
                                                    dd.rb_generator, descr.L_saturation);
      Lvel_ei_gen = SnapshotsGenerator.SpaceOpEvals(dmodel, 'L_velocity', ...
                                                    dd.rb_generator, descr.L_velocity);

      Ldiv_rb_gen = SnapshotsGenerator.SpaceOpEvals(dmodel, 'L_divergence', ...
                                                    dd.rb_generator, descr.L_divergence);
      dd.ei_generators = { Lsat_ei_gen, Lvel_ei_gen, Ldiv_rb_gen };

      if isequal(dd.train_sample_mode, 'single_deim')
        dd.datatree = gen_single_deim_detailed_data(dd, rmodel, model_data, checkpoint);
      elseif isequal(dd.train_sample_mode, 'single')
        dd.datatree   = gen_single_detailed_data(dd, rmodel, model_data, checkpoint);
      else
        error('no basis generation algorithm for RB implemented yet, use train_sample_mode = {single_deim, single}!');
%        dd.datatree     = gen_detailed_data(dd.bg_algorithm, rmodel, model_data, checkpoint);
      end
    end

  end

  methods

    function datatree = gen_single_deim_detailed_data(this, rmodel, model_data, checkpoint)
      if nargin <=3 || isempty(checkpoint)
        [ei_dd_map, checkpoint] = Greedy.Checkpoint.restore_latest_if_available(rmodel, 'ei_data');
      else
        ei_dd_map = restore_detailed_data(checkpoint, rmodel);
      end
      
      M_train = ParameterSampling.Single;
      M_train.init_sample(rmodel.detailed_model);

      set_mu(rmodel, M_train.sample(end,:));


    end
    
    function datatree = gen_single_detailed_data(this, rmodel, model_data, checkpoint)

      % M_train = ParameterSampling.Uniform(10);

      if nargin <=3 || isempty(checkpoint)
        [ei_dd_map, checkpoint] = Greedy.Checkpoint.restore_latest_if_available(rmodel, 'ei_data');
      else
        ei_dd_map = restore_detailed_data(checkpoint, rmodel);
      end

      M_train = ParameterSampling.Single;
      M_train.init_sample(rmodel.detailed_model);

      set_mu(rmodel, M_train.sample(end,:));

      if ~get(checkpoint, 'ei_generated', false)
        ei_plugin_sat = Greedy.Plugin.EI(this.ei_generators{1});
        ei_plugin_vel = Greedy.Plugin.EI(this.ei_generators{2});
        ei_plugin_sat.ei_target_error = 'interpol';
        ei_plugin_vel.ei_target_error = 'interpol';
        ei_greedy_sat = Greedy.Algorithm(ei_plugin_sat, M_train, []);
        ei_greedy_vel = Greedy.Algorithm(ei_plugin_vel, M_train, []);
        ei_greedy_sat.stop_epsilon = 1e-7;
        ei_greedy_vel.stop_epsilon = 1e-7;

        ei_dd_vel = gen_detailed_data(ei_greedy_vel, rmodel, model_data);
        ei_dd_sat = gen_detailed_data(ei_greedy_sat, rmodel, model_data);

        ei_dd_map = Greedy.DataTree.Detailed.IdMapNode({'L_saturation', 'L_velocity'});
        set(ei_dd_map, 1, ei_dd_sat);
        set(ei_dd_map, 2, ei_dd_vel);
      end
      checkpoint.store(rmodel, ei_dd_map, 'ei_data', struct('ei_generated', true));

      disp(num2str(get_mu(rmodel)));
      [sim_data, tictoc] = this.rb_generator.generate(rmodel.detailed_model, model_data);
      [Pdiv, tictoc] = this.ei_generators{3}.generate(rmodel.detailed_model, model_data);

%      RBS = sim_data.S;
%      RBU = sim_data.U(:,2:end);
%      RBP = Pdiv;

      if this.pca
        RBS = PCA_fixspace([sim_data.S], [], model_data.W, inf, 'qr', 1e-7);
        RBU = PCA_fixspace([sim_data.U], [], model_data.diamondWinv, inf, 'qr', 1e-6);
        RBP = PCA_fixspace([sim_data.P, Pdiv], [], model_data.W, inf, 'qr', 1e-6);
      else
        RBS = orthonormalize_qr([sim_data.S], model_data.W, 1e-6);
        RBU = orthonormalize_qr([sim_data.U], model_data.diamondWinv, 1e-6);
        RBP = orthonormalize_qr([sim_data.P, Pdiv], model_data.W, 1e-6);
      end

      rb_dd_sat = Greedy.DataTree.Detailed.RBLeafNode(model_data, 'saturation');
      rb_dd_vel = Greedy.DataTree.Detailed.RBLeafNode(model_data, 'velocity');
      rb_dd_prs = Greedy.DataTree.Detailed.RBLeafNode(model_data, 'pressure');
%      rb_dd_sat.RB = [ei_dd_sat.QM, PCA_fixspace(sim_data.S(:,1), ei_dd_sat.QM, model_data.W, inf, 'qr')];
%      rb_dd_vel.RB = ei_dd_vel.QM;
      rb_dd_sat.RB = RBS;
      rb_dd_vel.RB = RBU;
      rb_dd_prs.RB = RBP;

      rb_dd_map = Greedy.DataTree.Detailed.IdMapNode({'saturation', 'velocity', 'pressure'});
      set(rb_dd_map, 1, rb_dd_sat);
      set(rb_dd_map, 2, rb_dd_vel);
      set(rb_dd_map, 3, rb_dd_prs);

      datatree = TwoPhaseFlow.DetailedLeafNode(model_data);
      datatree.ei =  ei_dd_map;
      datatree.rb =  rb_dd_map;

      rmodel.N = create_tree(rb_dd_map, DataTree.ScalarCreator(@(X) get_rb_size(X)));
      rmodel.M = create_tree(ei_dd_map, DataTree.ScalarCreator(@(X) get_ei_size(X)));
    end

    function basis_gen = get.bg_algorithm(this)
      % generates a default @ref Greedy.Interface "basis generation object" for
      % linear evolution problems as given by a DetailedModel.
      %
      % This function generates a Greedy.Interface implementation, with either
      % a subsequent generation of empirical interpolation bases spaces for all
      % operators and a reduced basis space afterwards, @b or a combined
      % generation of all reduced spaces with the Greedy.Plugin.PODEI variant
      % of the Greedy.Algorithm.
      %

      basis_gen = [];
    end

    function rb_size = get_rb_size(this, rmodel)
      leaf_dd = get_leaf(this, rmodel);
      rb_size = get_rb_size(leaf_dd);
    end

  end
end