vi_gen_detailed_data.m

function detailed_data = vi_gen_detailed_data(model,model_data)
%function detailed_data = vi_gen_detailed_data(model,model_data)
%
% function generating the required detailed data for the RB-Method
% for Parametrized Variational Inequalities. currently two
% generation-modes are available: 
% - pure_snapshots
% - enrichment_by_supremizers
% both modes work with a uniform parameter grid and
% qr-orthonormalization of the reduced primal space.
%
% Generated fields of detailed_data:
%  RB_U: reduced space for the primal solution
%  RB_L: reduced space for the Lagrange multiplier

% I.Maier 30.05.2011

detailed_data = model_data;

tic;

switch model.RB_generation_mode
  
 case 'pure_snapshots'
  
  % generate uniformly distributed parameter values
  par.numintervals = model.RB_numintervals;
  par.range = model.mu_ranges;
  
  mu_mesh = cubegrid(par);
  mu_train = get(mu_mesh,'vertex');
  
  K = model.get_inner_product_matrix(model,model_data);
  
  U = zeros(model_data.H,size(mu_train,1));
  L = zeros(model_data.H,size(mu_train,1));
  for m = 1:size(mu_train,1)
    
    model = set_mu(model,mu_train(m,:));  
    sim_data = detailed_simulation(model,model_data);
    U(:,m) = sim_data.U;
    L(:,m) = sim_data.L;
  end;
  
  detailed_data.RB_U = improve_conditioning(U,K);
  detailed_data.RB_L = L;
  detailed_data.RB_info.mu_train = mu_train;
  
 case 'enrichment_by_supremizers'
  
  % generate uniformly distributed parameter values
  par.numintervals = model.RB_numintervals;
  par.range = model.mu_ranges;
  
  mu_mesh = cubegrid(par);
  mu_train = get(mu_mesh,'vertex');
  
  K = model.get_inner_product_matrix(model,model_data);
  
  U = zeros(model_data.H,size(mu_train,1)*2);
  L = zeros(model_data.H,size(mu_train,1));
  for m = 1:size(mu_train,1)
    
    model = set_mu(model,mu_train(m,:));  
    sim_data = detailed_simulation(model,model_data);
    U(:,2*m-1) = sim_data.U;
    U(:,2*m) = K\sim_data.L;
    L(:,m) = sim_data.L;
  end;
  
  detailed_data.RB_U = orthonormalize(U,K,1e-16,'qr');
  detailed_data.RB_L = L;
  detailed_data.RB_info.mu_train = mu_train;

   case 'separate-greedy'
  % separate generation of U and lambda basis
  % by projection error minimization
  
  % generate uniformly distributed parameter values
  par.numintervals = model.RB_numintervals;
  par.range = model.mu_ranges;
  
  mu_mesh = cubegrid(par);
  mu_train = get(mu_mesh,'vertex');
  ntrain = size(mu_train,1);
  
  K = model.get_inner_product_matrix(model,model_data);
  
  disp('computing training snapshots');
  [Utrain, Ltrain ] = filecache_function(...
      @compute_training_snapshots,model,model_data,mu_train);

%  params.plot = @(data,grid,params) plot(data);
%  detailed_data.RB_L = Ltrain;
%  detailed_data.RB_U = Utrain;
%  plot_detailed_data(model,detailed_data);
%  keyboard;
  
  % initialize empty bases
  U = zeros(model_data.H,0);
  L = zeros(model_data.H,0);
  max_u_err_sequence = [];  
  max_lambda_err_sequence = [];  
  mu_u_sequence = zeros(length(model.mu_names),0);  
  mu_lambda_sequence = zeros(length(model.mu_names),0);  
  mu_u_index_sequence = [];
  mu_lambda_index_sequence = [];
  
  % greedy loop for u
  disp('Greedy loop for u');
  err = 1e10;
  N = 0;
  % initial mu: middle value:
  max_i = (length(mu_train)+1)/2;
  mu_max = mu_train(max_i,:);
  continue_loop =   (err > model.RB_stop_epsilon && N < ...
                     model.RB_stop_Nmax);
  
  while continue_loop
    model = set_mu(model,mu_max);  
    sim_data = detailed_simulation(model,model_data);
    U = [U,sim_data.U];
    mu_u_sequence = [mu_u_sequence,mu_max(:)];
    mu_u_index_sequence = [mu_u_index_sequence, max_i]
    
    % (the following should be done without orthonormalization in
    % the future...)
    % compute projection error:
    UO = orthonormalize_loop(U,K,2e-16);
    %      UO = improve_conditioning(U,K);
    PUtrain = UO * (UO' * K * Utrain);
    err_u = sqrt(sum((Utrain-PUtrain).* (K*(Utrain-PUtrain))));
    if max(err_u(mu_u_index_sequence))>1e-8
      disp('accuracy problem in orthogonalization, please check')
      keyboard;
    else % set to zero, such that not selected again
      err_u(mu_u_index_sequence);
    end;
    [err, max_i] = max(err_u);
    if ~isempty(find(mu_u_index_sequence==max_i))
      disp('selected mu twice, training set exhausted')
      N = N+1;
      continue_loop = 0;
    else
      mu_max = mu_train(max_i,:);
      max_u_err_sequence = [max_u_err_sequence, err];
      N = N+1;
      continue_loop = err > model.RB_stop_epsilon && N < model.RB_stop_Nmax;
    end;
    disp(['N = ',num2str(N),...
          ', err_u = ',num2str(err_u(max_i))]);
  end;
  
  % greedy loop for lambda
  disp('Greedy loop for lambda');
  err = 1e10;
  N = 0;
  % initial mu: middle value:
  mu_max = mu_train((length(mu_train)+1)/2,:);
  max_i = (length(mu_train)+1)/2;
  mu_lambda_index_sequence = []; 
  continue_loop = (err > model.RB_stop_epsilon && N < model.RB_stop_Nmax);
  while continue_loop
    model = set_mu(model,mu_max);  
    sim_data = detailed_simulation(model,model_data);
    L = [L,sim_data.L];
    Kinv_L = K\L;
    mu_lambda_sequence = [mu_lambda_sequence,mu_max(:)];
    mu_lambda_index_sequence = [mu_lambda_index_sequence, max_i];
    
    % (the following should be done without orthonormalization in
    % the future...)
    %    err_lambda = sqrt(sum((Ltrain-PLtrain).* (K*(Ltrain-PLtrain))));
    %LO = K * orthonormalize_loop(Kinv_L,K,2e-16);
    %    Kinv_Ltrain = K\Ltrain;
    %    PLtrain = LO * (LO' * Kinv_Ltrain);
    Kinv_L = K\L;
    PLtrain = Ltrain;
    % for each single snapshot compute projection by small QP....
    %%%  ordinary l2-projection:
    %%% PLtrain = L * pinv(L'*Kinv_L)* (Kinv_L' * Ltrain);
%    warning off
    H = 2*L'*Kinv_L;
    H = 0.5*(H+H');
    b = zeros(N,1);
    A = -eye(N,N);
    options = optimset('TolFun',0); % relative function change limit
    options.MaxIter = 1000;
    options.Display='off';
    for m = 1:ntrain;
      f = -2*Ltrain(:,m)'*Kinv_L;
      fprintf('*');
      [alphaN,Fval,Exitflag] =  quadprog(H,f,[],[],[],[],0,Inf,[],options);
%      keyboard;
      PLtrain(:,m) = L*alphaN;
%      if m == max_i % error should be minimum!
%       plot([L*alphaN,Ltrain(:,m)]);
%       keyboard;
%      end;
    end;
    fprintf('\n');
%    warning on
    Kinv_PLerr = K\(Ltrain-PLtrain);
    err_lambda = sqrt(sum((Ltrain-PLtrain).* (Kinv_PLerr)));
    % set earlier errors to 0...
    err_lambda(mu_lambda_index_sequence)= 0;
    [err, max_i] = max(err_lambda);
    if ~isempty(find(mu_lambda_index_sequence==max_i))
%      disp('selected mu twice, training set exhausted? Please check')
%      keyboard;
      %      m = max_i;
      %      f = -2*Ltrain(:,m)'*Kinv_L;
      %      [alphaN,Fval,Exitflag] =  quadprog(H,f,[],[],[],[],0,Inf,[],options);
      %      %      keyboard;
      %      plot([L*alphaN,Ltrain(:,m)]);
      %      disp('error from optimization functional:')
      %      sqrt(Fval + Ltrain(:,m)'*(K\Ltrain(:,m)))
      %      disp('error from norm:')
      %      Ldiff = L*alphaN-Ltrain(:,m);
      %      sqrt(Ldiff'*(K\Ldiff))
      N = N+1;
      continue_loop = 0;
    else
      mu_max = mu_train(max_i,:);
      max_lambda_err_sequence = [max_lambda_err_sequence, err];
      N = N+1;
      continue_loop = (err > model.RB_stop_epsilon && N < model.RB_stop_Nmax);
    end;
    disp(['N = ',num2str(N),...
          ', err_lambda = ',num2str(err_lambda(max_i))]);
  end;
   
  if model.supremizer_enrichment
    % include supremizers into primal basis
    U = [U, K\L];  
    % U_part = [U(:,1:n_part), K\(L(:,1:n_part))];  
    detailed_data.RB_U_with_supr = U;
  end;  

  detailed_data.RB_U = improve_conditioning(U,K);
  detailed_data.RB_L = L;
  detailed_data.RB_info.max_u_err_sequence = max_u_err_sequence;
  detailed_data.RB_info.max_lambda_err_sequence = max_lambda_err_sequence;
  detailed_data.RB_info.mu_u_sequence = mu_u_sequence;
  detailed_data.RB_info.mu_lambda_sequence = mu_lambda_sequence;
  detailed_data.RB_info.mu_lambda_index_sequence = ...
      mu_lambda_index_sequence; 

 case {'greedy-true-projection-error','greedy-true-rb-error',...
       'greedy-rb-estimator'};
  
  if ~isfield(model,'dual_lin_independence_mode')
    model.dual_lin_independence_mode = 0;
  end;
  
  if isequal(model.RB_generation_mode, ...
             'greedy-true-projection-error')
    use_projection_error = 1;
  else
    use_projection_error = 0;
  end;
  
  if isequal(model.RB_generation_mode, ...
             'greedy-true-rb-error')
    use_true_rb_error = 1;
  else
    use_true_rb_error = 0;
  end;
  
  % generate uniformly distributed parameter values
  par.numintervals = model.RB_numintervals;
  par.range = model.mu_ranges;
  
  mu_mesh = cubegrid(par);
  mu_train = get(mu_mesh,'vertex');
  ntrain = size(mu_train,2);
  
  K = model.get_inner_product_matrix(model,model_data);
  
  if (use_true_rb_error ~= 0) | (use_projection_error ~= 0)    
    disp('computing training snapshots');
    [Utrain, Ltrain ] = filecache_function(...
        @compute_training_snapshots,model,model_data,mu_train);
  end;
  
  % initialize empty bases
  U = zeros(model_data.H,0);
  L = zeros(model_data.H,0);
  max_err_sequence = [];  
  mu_sequence = zeros(length(model.mu_names),0);  
  mu_index_sequence = [];
  max_u_err_sequence = [];  
  max_lambda_err_sequence = [];  
  
  % greedy loop 
  err = 1e10;
  N = 0;
  % initial mu: middle value:
  mu_max = mu_train((length(mu_train)+1)/2,:);
  max_i = (length(mu_train)+1)/2;
  
  while err > model.RB_stop_epsilon && N < model.RB_stop_Nmax;
    model = set_mu(model,mu_max);  
    sim_data = filecache_function(@detailed_simulation,model,model_data);
%    sim_data = detailed_simulation(model,model_data);
    U = [U,sim_data.U];
    L = [L,sim_data.L];
    mu_sequence = [mu_sequence,mu_max(:)];
    mu_index_sequence = [mu_index_sequence,max_i];
    if use_projection_error
      % (the following should be done without orthonormalization in
      % the future...)
      % compute projection error:
      UO = orthonormalize_loop(U,K,2e-16);
      %      UO = improve_conditioning(U,K);
      Kinv_L = K\L;
      PUtrain = UO * (UO' * K * Utrain);
      
      LO = K * orthonormalize_loop(Kinv_L,K,2e-16);
      %    Kinv_Ltrain = K\Ltrain;
      %    PLtrain = LO * (LO' * Kinv_Ltrain);
      Kinv_L = K\L;
      PLtrain = L * pinv(L'*Kinv_L)* (Kinv_L' * Ltrain);
      
      Kinv_PLerr = K\(Ltrain-PLtrain);
      
      err_u = sqrt(sum((Utrain-PUtrain).* (K*(Utrain-PUtrain))));
      %    err_lambda = sqrt(sum((Ltrain-PLtrain).* (K*(Ltrain-PLtrain))));
      err_lambda = sqrt(sum((Ltrain-PLtrain).* (Kinv_PLerr)));
    else % use true rb error or estimator
      
      % construct reduced basis (with supr.) for rb simulation:
      detailed_data = model_data;
%      detailed_data = [];
      if model.dual_lin_independence_mode == 0 %
        detailed_data.RB_L = L;
      elseif model.dual_lin_independence_mode == 1
        dual_speye = speye(size(L,1));
        L_dofs = find(max(abs(L),[],2)>1e-11);
%       keyboard;
        detailed_data.RB_L = full(dual_speye(:,L_dofs));
        disp(['Found Linear independent dual basis of size ' ,...
              num2str(length(L_dofs)),' for ',...
              num2str(size(L,2)),' dual snapshots']);
      end;
      if model.supremizer_enrichment
        %       detailed_data.RB_U = [U,K\L];
        %       detailed_data.RB_U = orthonormalize_loop([U, K\L] ,K,2e-16);
        detailed_data.RB_U = improve_conditioning(...
            [U, K\detailed_data.RB_L],K);
      else
        %       detailed_data.RB_U = U;
        %       detailed_data.RB_U = orthonormalize_loop([U] ,K,2e-16);
        detailed_data.RB_U = improve_conditioning([U],K);
      end;
%      detailed_data.X = model_data.X;
%      detailed_data.dx = model_data.dx;
      disp(['Size of U without enrichment: ',num2str(size(U,2)),...
            ', Size of orth(U + possible supr.) : ',num2str(size(detailed_data.RB_U,2))])
      
      if use_true_rb_error      
        UNtrain = zeros(model_data.H,ntrain);
        LNtrain = zeros(model_data.H,ntrain);
      else % use rb error estimator
        err_u= zeros(1,size(mu_train,1));
        err_lambda= zeros(1,size(mu_train,1));
        model.enable_error_estimator = 1;
      end;
      reduced_data = gen_reduced_data(model, detailed_data);
      for m = 1:size(mu_train,1)
        model = set_mu(model,mu_train(m,:));
        sim_data = rb_simulation(model,reduced_data);
        if use_true_rb_error    
          UNtrain(:,m) = detailed_data.RB_U * sim_data.UN;
          LNtrain(:,m) = detailed_data.RB_L * sim_data.LN;
        else % use rb_error_estimator
          err_u(m) = sim_data.Delta_u;
          err_lambda(m) = sim_data.Delta_lambda;
        end;
      end;
      if use_true_rb_error      
        err_u = sqrt(sum((Utrain-UNtrain).* (K*(Utrain-UNtrain))));
        %    err_lambda = sqrt(sum((Ltrain-PLtrain).* (K*(Ltrain-PLtrain))));
        Kinv_Lerr = K\(Ltrain-LNtrain);
        err_lambda = sqrt(sum((Ltrain-LNtrain).* (Kinv_Lerr)));
      end;
    end;
    
    % set errors of recent values to 0 in order to prevent "doubles"
    err_u(mu_index_sequence) = 0;
    err_lambda(mu_index_sequence) = 0;
    [err, max_i] = max(model.w_u * err_u + model.w_lambda * ...
                       err_lambda);
    disp(['N = ',num2str(N),...
          ', err = ',num2str(err),...
          ', err_u = ',num2str(err_u(max_i)),...
          ', err_lambda = ',num2str(err_lambda(max_i)) ]);
    mu_max = mu_train(max_i,:);
    max_err_sequence = [max_err_sequence, err];
    max_u_err_sequence = [max_u_err_sequence, err_u(max_i)];
    max_lambda_err_sequence = [max_lambda_err_sequence, err_lambda(max_i)];
    N = N+1;
  end;

  % return results
  detailed_data = model_data;
 
  detailed_data.RB_U_no_supr = U;
  %n_part = ceil(size(U,2)/2);
  %U_part = [U(:,1:n_part)];
  
  if model.dual_lin_independence_mode == 0 %
    detailed_data.RB_L = L;
  elseif model.dual_lin_independence_mode == 1
    dual_speye = speye(size(L,1));
    L_dofs = find(max(abs(L),[],2)>1e-11);
    detailed_data.RB_L = full(dual_speye(:,L_dofs));
    detailed_data.L_no_post_processing = L; 
  end;
  
  if model.supremizer_enrichment
    % include supremizers into primal basis
    U = [U, K\detailed_data.RB_L];  
    % U_part = [U(:,1:n_part), K\(L(:,1:n_part))];  
    detailed_data.RB_U_with_supr = U;
  end;  
  
  % no orthonormalization, works for basis up to 50:
  %detailed_data.RB_U = U;
  
  % numerical orthonormalization may lead to additional approximation error
  %  detailed_data.RB_U = orthonormalize_loop(U,K,2e-16);
  detailed_data.RB_U = improve_conditioning(U,K);
  % return part of the basis to see, whether order is maintained
  %detailed_data.RB_U_part = orthonormalize_loop(U_part,K,2e-16);
  
  %  detailed_data.RB_L = L;
  detailed_data.RB_info.max_err_sequence = max_err_sequence;
  detailed_data.RB_info.max_u_err_sequence = max_u_err_sequence;
  detailed_data.RB_info.max_lambda_err_sequence = max_lambda_err_sequence;
  detailed_data.RB_info.mu_sequence = mu_sequence;
  
 otherwise
  error('RB generation mode unknwon');
end;

t = toc;
detailed_data.RB_info.computation_time = t;

function [Utrain,Ltrain] = compute_training_snapshots(model,model_data,mu_train)
disp('computing training snapshots');
% compute training set
Utrain = zeros(model_data.H,size(mu_train,1));
Ltrain = zeros(model_data.H,size(mu_train,1));
for m = 1:size(mu_train,1)    
  model = set_mu(model,mu_train(m,:));  
  sim_data = detailed_simulation(model,model_data);
  Utrain(:,m) = sim_data.U;
  Ltrain(:,m) = sim_data.L;
  fprintf('.');
end;
fprintf('\n');

function UO = orthonormalize_loop(U,K,eps);

%if (cond(U'*K*U)>1e6);
  disp('orthogonalizing!');
  
  % perform not exact orthongonalization, but 
  % conditioning improvement.
%  keyboard;
  %[u,s,v] = svd(U,0);
  %UO = u;
  %unorminv = (sqrt(sum((K * UO).*UO,1))).^(-1);
  %UO = UO * diag(unorminv(:));
  
  %% the following performs repeated qr orthogonalization
  UO = orthonormalize_qr(U,K,eps);
  KN = UO' * K * UO;
  e = (max(max(abs(KN-eye(size(KN))))));
  while e > 1e-12
    disp('Detected orthonormalization inaccuracy.');
    disp('PERFORMING REPEATED ORTHONORMALIZATION');
    UO = orthonormalize_qr(UO,K,2e-16);
    KN = UO' * K * UO;
    e = (max(max(abs(KN-eye(size(KN))))));
    %    error('error in orthonormalization, please check!');
  end;

%else
%  UO = U;
%end;