POD.m

classdef POD < Greedy.Plugin.PODCommon
  % extension class implementing the POD-greedy algorithm extension routines
  %

  properties(Constant)
    % cell array of field names which shall be copied to the detailed data node
    % structure on creation in init_basis() method.
    info_fields = {'stop_Nmax', 'noof_pca_modes'};
  end

  properties
    % maximum number of generated reduced basis vectors.
    stop_Nmax      = inf;

    % number of most significant POD modes to be added to the reduced basis in
    % each extension step.
    %
    % See also: basis_extension() for more information.
    noof_pca_modes = 1;
  end

  properties(Constant)
    % @copybrief ::GreedyInterfacegenerated_basis_type
    %
    % @copydetails ::GreedyInterfacegenerated_basis_type
    %
    generated_basis_type = 'rb';
  end

  methods

    function rbpcafe = POD(generator, id)
      % function rbpcafe = POD(generator)
      % constructor for the POD-greedy extension algorithms.
      %
      % Parameters:
      %   generator: object of type SnapshotsGenerator.Trajectories generating
      %              the (high dimensional) basis functions, i.e. solutions of
      %              the numerical scheme for the parametrized partial
      %              differential equation.

      rbpcafe = rbpcafe@Greedy.Plugin.PODCommon(generator);
    end

    function detailed_data = init_basis(this, rmodel, model_data, M_train)
      % function detailed_data = init_basis(this, rmodel, model_data, M_train)
      % @copybrief ::GreedyPlugin::Interfaceinit_basis()
      %
      % @copydetails GreedyPlugin::Interfaceinit_basis()
      %
      % Parameters:
      %   rmodel: of type ::GreedyUser::IReducedModel
      %
      % Return values:
      %   detailed_data: of type ::GreedyDataTreeDetailed.RBLeafNode
      %
      % Computes the initial data functions for all parameters from the
      % training set 'M_train' and initializes the reduced basis with an
      % orthonormalization of these DOF vectors. The result is stored in a
      % detailed data tree leaf node object.
      if isa(model_data, 'Greedy.DataTree.Detailed.INode')
        model_data = model_data.model_data;
      end
      dmodel = rmodel.detailed_model;
      nmus = size(M_train, 1);
      RBinit = zeros(model_data.grid.nelements, nmus);
      for i = 1:nmus
        dmodel = set_mu(dmodel, M_train(i,:));
        dmodel.decomp_mode = 0;

        U0 = init_values_algorithm(dmodel, model_data);
        RBinit(:,i) = U0;
      end

      RBinit = orthonormalize(dmodel, model_data, RBinit);

      if size(RBinit, 2) == 0
        RBinit = ones(size(RBinit,1),1);
      end
      if isstruct(model_data)
        detailed_data    = Greedy.DataTree.Detailed.RBLeafNode(model_data, this.id);
        set_fields(detailed_data, this, Greedy.Plugin.POD.info_fields);
      end
      detailed_data.RB = RBinit;
      disp(['found ',num2str(size(RBinit,2)),' basis functions for',...
            ' init data variation.']);
    end

    function prepare_reduced_data(this, rmodel, detailed_data)
      % function prepare_reduced_data(this, rmodel, detailed_data)
      % @copybrief GreedyPlugin::Interfaceprepare_reduced_data()
      %
      % @copydetails GreedyPlugin::Interfaceprepare_reduced_data()
      %
      % Parameters:
      %   rmodel: of type ::GreedyUser::IReducedModel
      %   detailed_data: of type ::GreedyDataTreeDetailedRBLeafNode

%      if isempty(this.last_detailed_data_handle_for_prepare_reduced) ...
%          || this.last_detailed_data_handle_for_prepare_reduced ~= detailed_data
        this.reduced_data = gen_reduced_data(rmodel, detailed_data);
%        this.last_detailed_data_handle_for_prepare_reduced = detailed_data;
%      end
    end

    function Uapprox = generate_reduced(this, rmodel, reduced_data, detailed_data, U)
      % function Uapprox = generate_reduced(this, rmodel, reduced_data, detailed_data, U)
      % @copybrief GreedyPlugin::Interfacegenerate_reduced()
      %
      % @copydetails GreedyPlugin::Interfacegenerate_reduced()
      %
      % Parameters:
      %   rmodel: of type ::GreedyUser::IReducedModel
      %   detailed_data: of type ::GreedyUser::IDetailedData
      %   reduced_data: of type ::GreedyUser::IReducedData
      if isempty(reduced_data)
        error('prepare_reduced_data() may need to be executed before generate_reduced()!');
      end
      rmodel.N = get_rb_size(detailed_data);
      rmodel.M = get_ei_size(detailed_data) - rmodel.Mstrich;

      rb_sim_data = rb_simulation(rmodel, reduced_data);
      rb_sim_data = rb_reconstruction(rmodel, detailed_data, rb_sim_data);
      Uapprox = rmodel.get_dofs_from_sim_data(rb_sim_data);
    end

    function [max_errs, max_err_sequence, max_mu_index] = ...
        error_estimators(this, rmodel, detailed_data, M_train)
      % function [max_errs, max_err_sequence, max_mu_index] = error_estimators(this, rmodel, detailed_data, M_train)
      % @copybrief GreedyPlugin::Defaulterror_estimators()
      %
      % @copydetails GreedyPlugin::Defaulterror_estimators()
      %
      % Parameters:
      %   rmodel: of type ::GreedyUser::IReducedModel
      %   detailed_data: of type ::GreedyUser::IDetailedData
      %
      % The estimator `\eta_{N,M,M'}` used here is described in @ref DHO11
      % "[DHO11], Chapter 5".
      %
      nmus = size(M_train, 1);

      tmp_max_errs = cell(nmus, 1);
      tmp_max_err  = zeros(nmus, 1);

      reduced_data = this.reduced_data;
      if rmodel.crb_enabled
        if isa(reduced_data.M, 'DataTree.INode')
          rmodel.M = create_scalar_tree(reduced_data.M, @(x) get(x,1) - rmodel.Mstrich);
        else
          rmodel.M = reduced_data.M - rmodel.Mstrich;
        end
      end

      if ~FakeMPI.MPI_IsActive
        ME.message = '';
        parfor i = 1:nmus
        %for i = 1:nmus
          trmodel = rmodel;
          if trmodel.verbose >= 5
            fprintf('.');
          end
          set_mu(trmodel, M_train(i,:));
          try
            sim_data        = rb_simulation(trmodel, reduced_data);
            tmp_max_errs{i} = trmodel.get_estimators_from_sim_data(sim_data);
            tmp_max_err(i)  = max(tmp_max_errs{i});
          catch ME
            tME = ME;
            warning('RBmatlab:simulation_error', ...
                    ['Method rb_simulation threw an error: ', tME.message]);
            tmp_max_errs{i} = NaN;
            tmp_max_err(i)  = NaN;
          end
        end
        fprintf('\n');
      else
        buf.reduced_data = reduced_data;
        buf.rmodel       = rmodel;
        buf.M_train      = M_train;
        FakeMPI.MPI_Deploy(buf);
        buf    = [];

        min_mu = 0;
        nnodes = min(FakeMPI.MPI_Size, nmus-min_mu);
        width  = max((nmus-min_mu)/(2*nnodes),1);

        ids  = 1:nnodes;
        bufs = cell(1,nnodes);
        for id=1:nnodes
          bufs{id}.st   = round((id-1)*width)+1+min_mu;
          bufs{id}.en   = min(round(id*width)+min_mu, nmus);
%          disp(['Send: id = ', num2str(id), ' st = ', num2str(buf.st), ' en = ', num2str(buf.en)]);
        end
        FakeMPI.MPI_Send(ids, bufs);
        min_mu = bufs{end}.en;

        running=true(1,nnodes);
        countdown = 0;

        pause(0.5);

        while (any(running))
          probed_ids = FakeMPI.MPI_Probe('*');
          id_mask = find(running);
          for id=intersect(probed_ids, id_mask)
            ret_struct = FakeMPI.MPI_Recv(id);
            try
              st   = ret_struct.st;
              en   = ret_struct.en;
              assert(bufs{id}.st == st)
              assert(bufs{id}.en == en)
              tmp_max_errs(st:en) = ret_struct.tmp_max_errs(:);
              tmp_max_err(st:en)  = ret_struct.tmp_max_err(:);
            catch exception
              disp('received corrupt data');
              disp(exception.message);
              disp('resend this job');
              FakeMPI.MPI_ReDeploy(id);
              FakeMPI.MPI_Send(id, bufs{id});
              pause(0.05);
              continue;
            end
            if min_mu < nmus
              if countdown == 0
                width = max((nmus-min_mu)/(2*nnodes), 1);
                countdown = nnodes;
              end
              buf.st = min_mu + 1;
              buf.en = min(min_mu + round(width), nmus);
              bufs{id} = buf;
%              disp(['Send: id = ', num2str(id), ' st = ', num2str(buf.st), ' en = ', num2str(buf.en)]);
              FakeMPI.MPI_ReDeploy(id);
              FakeMPI.MPI_Send(id, buf);
              pause(0.05);
              min_mu = buf.en;
              countdown = countdown - 1;
            else
              running(id) = 0;
              disp(['stopped parCompErrs #', num2str(id)]);
            end
          end
          pause(0.05);
        end
        if getenv('PBS_JOBID')
          pause(0.5);
        end

      end

      [dummy, max_mu_index] = max(tmp_max_err);
      max_errs              = tmp_max_err;
      max_err_sequence      = tmp_max_errs{max_mu_index};
    end

    function [breakloop, reason] = pre_check_for_end(this, rmodel, detailed_data)
      % function [breakloop, reason] = pre_check_for_end(this, rmodel, detailed_data)
      % @copybrief Greedy.Plugin.Interface.pre_check_for_end()
      %
      % @copydetails Greedy.Plugin.Interface.pre_check_for_end()
      %
      % Parameters:
      %   rmodel: of type .Greedy.User.IReducedModel
      %   detailed_data: of type .Greedy.DataTree.Detailed.RBLeafNode

      dd_leaf = get_active_leaf(detailed_data, rmodel);
      reason = '';
      breakloop = false;
      if stopped_on_active_child(detailed_data, 'stopped_on_empty_extension', rmodel)
        breakloop = true;
        reason = 'RB extension stopped on empty extension';
      elseif size(detailed_data.RB, 2) >= get_field(dd_leaf, 'stop_Nmax')
        set_stop_flag(dd_leaf, 'stopped_on_max_basis_size');
        breakloop = true;
        reason = 'RB extension reached maximum basis size';
      end
    end


%    breakloop = post_check_for_end(this, detailed_data);

    function ret_detailed_data = basis_extension(this, rmodel, detailed_data, max_err_seq, mu)
      % function ret_detailed_data = basis_extension(this, rmodel, detailed_data, max_err_seq, mu)
      % @copybrief Greedy.Plugin.Interface.basis_extension()
      %
      % @copydetails Greedy.Plugin.Interface.basis_extension()
      %
      % Parameters:
      %   rmodel: of type .Greedy.User.IReducedModel
      %   detailed_data: of type .Greedy.DataTree.Detailed.RBLeafNode
      %
      % Return values:
      %   detailed_data: of type .Greedy.DataTree.Detailed.RBLeafNode
      %
      % Implementation details:
      %
      % - This method computes first POD modes of this sequence of discrete
      % functions and adds them to the reduced basis space.
      % - The number of POD modes used can be specified by the #noof_pca_modes
      % property.
      % - The POD is computed by the PCA_fixspace() function.

      % extend rb basis if you receive rbeidetailedata
      if isa(detailed_data, 'Greedy.DataTree.Detailed.PODEILeafNode')
        ret_detailed_data = detailed_data;
        detailed_data = ret_detailed_data.rb;
      else
        ret_detailed_data = detailed_data;
      end
      model_data = detailed_data.model_data;
      max_err = max(max_err_seq);
      disp(['Detected maximum error prediction ',num2str(max_err),...
            ' for mu=[',...
            num2str(mu),']']);

      dmodel = rmodel.detailed_model;
      set_mu(dmodel, mu);
      U     = generate(this.generator, dmodel, model_data);

      cur_stop_Nmax      = get_field(detailed_data, 'stop_Nmax');
      cur_noof_pca_modes = get_field(detailed_data, 'noof_pca_modes');
      pca_modes = min(cur_noof_pca_modes, cur_stop_Nmax - size(detailed_data.RB,2));
      W = dmodel.get_inner_product_matrix(model_data);
      RBext     = PCA_fixspace(U, detailed_data.RB,...
                               W, pca_modes);
      if ~isempty(RBext)
        RB        = [detailed_data.RB, RBext];
        %        RBo       = model_orthonormalize_qr( dmodel.descr, model_data, RB, eps );
        %        K         = inner_product(dmodel, model_data, RBo, RB);
        %        [dummy,i] = max(abs(K));
        %        RB        = RBo(:,i);
        %        if model.debug
        %          K  = rb_engine.inner_product(model,detailed_data, RB, RB);
        %          if max(max(abs(K-eye(size(K)))))>1e-10
        %            error('error in orthonormal basis extension!');
        %          end
        %        end
        detailed_data.RB = RB;
        N = size(detailed_data.RB, 2);
        disp(['Extended RB to length ', num2str(N)]);
      else
        set_stop_flag(detailed_data, 'stopped_on_empty_extension');
        %! \todo Remember to catch this one in the container above when
        % combining eidetailed and pod-greedy!!!!
      end
    end
  end
end