gamm2013_exp.m

function gamm2013_exp(step, exp_params)
%function gamm2013_exp(step)
%
% experiments for the GAMM2013 paper: optimization with 2x2 thermal block

% B. Haasdonk, M. Dihlmann, 22.5.2013

if nargin<1
     step = 1;
    exp_params =[];
elseif nargin<2
    exp_params = [];
end;

params.B1 = 2; 
params.B2 = 2;
params.numintervals_per_block = 50;
model = thermalblock_model_struct(params);
model.RB_stop_epsilon = 1e-6;
model.RB_stop_Nmax = 30;
model.RB_train_size = 5000;
model_data = gen_model_data(model);
plot_params = [];
plot_params.axis_tight = 1;
plot_params.yscale_uicontrols = 0.5;

switch step
 case 1 % 2x2 thermal block, check detailed simulation by gui
  
  %  model = model.set_mu(model,[0.1,0.8,0.1,0.1]);
  %  sim_data = detailed_simulation(model,model_data);
  demo_detailed_gui(model, model_data,[],plot_params);
  %  detailed_data = gen_detailed_data(model,model_data);
  %  reduced_data = gen_reduced_data(model, detailed_data);
  %  demo_rb_gui(model,detailed_data,reduced_data);

 case 2 % reference solution for optimization
  mu_star = [2,1,1,2];
  model = model.set_mu(model,mu_star);
  sim_data = detailed_simulation(model,model_data);
  plot_sim_data(model,model_data,sim_data,plot_params);
  
 case 3 %detailed data generation, mu3, mu4 fixed
  mu3 = 1; mu4 = 2;
  % keep parameters mu3 and mu4 fixed in greedy basis generation  
  %model.RB_mu_list = 
  model.mu_ranges = {[0.5,2.5],[0.5,2.5],[mu3,mu3],[mu4,mu4]};
  %id = 'MATLAB:nearlySingularMatrix';
  %warning('off',id)
  % set random number generator for deterministic result:
  RandStream.setDefaultStream(RandStream('mt19937ar','seed',1010101));
  detailed_data = gen_detailed_data(model, model_data);
  %warning('on',id)  
  % 15 Basisvektoren erzeugt

  save('gamm2013_detailed_data.mat');

 case 4 % inspect reduced basis and rb_demo_gui
  
  load('gamm2013_detailed_data');  
  params.plot = @plot_vertex_data;
  plot_sequence(detailed_data.RB, model_data.grid, params)
  figure;
  
  demo_rb_gui(model,detailed_data,reduced_data,plot_params);  
  
  figure;
  
 case 5 % add optimization functionals to detailed_data
  load('gamm2013_detailed_data');
  
  %define optimization functionals by error to reference solution
  mu_star = [2,1];
  model = model.set_mu(model,mu_star);
  ref_sim_data = detailed_simulation(model,model_data);
  ref_dofs = ref_sim_data.uh.dofs;
  model.set_mu = @(model,mu)my_pfusch_set_mu(model,mu,mu3,mu4);
  model.get_mu = @(model)model.mus(1:2);
  J = @(dofs) J_functional(dofs, ...
                           model_data.df_info.l2_inner_product_matrix, ...
                           ref_dofs);
  detailed_optimization_functional = @(mu1_mu2) ... 
      my_detailed_optimization_functional(mu1_mu2, J, model, model_data);

  reduced_data = gen_reduced_data(model,detailed_data);
  %add functinal fields in reduced data
  RB = detailed_data.RB;
  M = model_data.df_info.l2_inner_product_matrix;
  reduced_data.J1 = RB'*M*RB;
  reduced_data.J2 = ref_dofs'*M*RB;
  reduced_data.J3 = ref_dofs'*M*ref_dofs;
  %JN = 0.5*(a'*J1*a-2*J2*a+J3);
  %gradient J = a'*gJ1*pmui_a - gJ2 pmui_a
  reduced_data.gJ1 = RB'*M*RB;
  reduced_data.gJ2 = ref_dofs'*M*RB;
  
  reduced_optimization_functional = @(mu1_mu2) ... 
      my_reduced_optimization_functional(mu1_mu2, model,...
                                         reduced_data);
    
  model_data.ref_dofs = ref_dofs;  
                 
  save('gamm2013_detailed_data.mat');
  
 case 6 % plot reduced output landscape and time measurement
  load('gamm2013_detailed_data');

  if isfield(exp_params, 'nmu')
      n_mu1 = exp_params.nmu;
      n_mu2 = exp_params.nmu;
  else
    n_mu1 = 21; n_mu2 = 21;
  end
%  mu1_range = [0.5,2.5];
%  mu1_range = [1.9,2.1];
  mu1_range = model.mu_ranges{1};
  mu1s = linspace(mu1_range(1),mu1_range(2),n_mu1);
  mu2_range = model.mu_ranges{2};
%  mu2_range = [0.5,2.5];
%  mu2_range = [0.9,1.1];
  mu2s = linspace(mu2_range(1),mu2_range(2),n_mu2);
  Js = zeros(n_mu1,n_mu2);
  tic;
  for n1 = 1:n_mu1
    for n2= 1:n_mu2  
      mu1_mu2=[mu1s(n1),mu2s(n2)];
      Js(n1,n2)= reduced_optimization_functional(mu1_mu2);
    end;
  end;
  
  t = toc;
  disp(['computation time: ',num2str(t),' seconds.']);
  
  
  cmaplog = create_log_colormap();
  colormap(cmaplog);
  surf(mu1s,mu2s,Js','FaceColor','interp');  
  xlabel('{\mu}_1');
  ylabel('{\mu}_2');

  % check reduced model and reduced optimization functional at
  % reference config:
  model = model.set_mu(model,mu_star); 
  sim_data = rb_simulation(model,reduced_data);  
  sim_data = rb_reconstruction(model, detailed_data, sim_data);
  disp(['H10 error to reference solution:']);
  err = sim_data.uh.dofs - ref_dofs;
  h10_error = sqrt(err'*model_data.df_info.h10_inner_product_matrix * err)
  disp(['Delta error estimator:']);
  sim_data.Delta
  disp(['L2 error to reference solution:']);
  l2_error = sqrt(err'*model_data.df_info.l2_inner_product_matrix * err)
  
  JN_ref = reduced_optimization_functional(mu_star(1:2));
  disp(['reduced optimization functional at reference parameter: ',...
        num2str(JN_ref)]);
  
  keyboard;
  
 case 7 % run detailed optimizer and time-measurement
  load('gamm2013_detailed_data');
  mu0 = [0.5,0.5];
  lb = [0.5,0.5];
  ub = [2.5,2.5];
%  options = [];
  options = optimset();
  options.MaxFunEvals = 1000;
  options.TolFun = 1e-14;
  options.TolX = 1e-14;
  fun = detailed_optimization_functional;
  tic;
  [mumin,fval,exitflag,output,lambda,grad] = fmincon(fun,mu0,[],[], ...
                                                  [],[],lb,ub,[],options);
  t = toc;  
  disp('found optimum mu_min:')
  mumin
  disp('J(mu_min): ');
  fun(mumin)
  disp('detailed optimization time:')
  t
  disp(['functional gradient at optimum: ',num2str(grad')]);
  
 case 8 % run reduced optimizer and time-measurement and error estimation
  load('gamm2013_detailed_data');
  mu0 = [0.5,0.5];
  lb = [0.5,0.5];
  ub = [2.5,2.5];
%  options = [];
  options = optimset();
  options.MaxFunEvals = 1000;
  options.TolFun = 1e-14;
  options.TolX = 1e-14;
  fun = reduced_optimization_functional;
  tic;
  [mumin,fval,exitflag,output,lambda,grad] = fmincon(fun,mu0,[],[], ...
                                                  [],[],lb,ub,[],options);
  t = toc;  
  format long;
  disp('found optimum mu_min:')
  mumin
  disp('J(mu_min): ');
  fun(mumin)
  disp('reduced optimization time:')
  t

  disp('l2_error in mu_min:') 
  norm(mumin-mu_star(1:2))

  disp('computing a-posteriori error estimator (expensive):')
  
  % infinity norm * sqrt2 is upper bound for 2-norm of gradient
  %epsilon_J = norm(grad);%output.firstorderopt*sqrt(2);
  grad_J = my_reduced_functional_gradient(model,...
                    reduced_data);
  epsilon_J = norm(grad_J);
  
  disp(['norm functional gradient at optimum: ',num2str(epsilon_J)]);
  

  % error estimator computation (expensive!!!) 
  
  [Delta_mu, val_indicator] = my_expensive_error_estimator(mumin,model, ...
                                          detailed_data,reduced_data,...
                                          detailed_optimization_functional,...
                                          epsilon_J, ...
                                          ref_dofs);
  
  disp('error bound:')
  Delta_mu
  
  disp('validity indicator, should be less than 1 for validity: ')
  val_indicator
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555
  % do a variation of number of basis vectors N in RB
  % track online simulation time
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 9 
        
        load('gamm2013_detailed_data');
        detailed_data_full = detailed_data;
        %J = @(dofs) J_functional(dofs, ...
        %    model_data.df_info.l2_inner_product_matrix, ...
        %    ref_dofs);
        %
        %detailed_optimization_functional = @(mu1_mu2) ...
        %    my_detailed_optimization_functional(mu1_mu2, J, model, model_data, ...
        %    mu3,mu4);
        
        %optimization settings
        mu0 = [0.5,0.5];
        lb = [0.5,0.5];
        ub = [2.5,2.5];
        %  options = [];
        options = optimset();
        options.MaxFunEvals = 1000;
        options.TolFun = 1e-14;
        options.TolX = 1e-14;
        
        %detailed_optimization
        fun = detailed_optimization_functional;
        tic;
        [mumin_det,fval_det,exitflag_det,output_det,lambda,grad_det] = fmincon(fun,mu0,[],[], ...
                                                  [],[],lb,ub,[],options);
        t_det = toc;  
        
        
        %N - vector
        Nmax = size(detailed_data.RB,2);
        N_vec = 14:1:Nmax;
        %initialization 
        Delta_vec = zeros(1,length(N_vec));
        Delta_pmui_vec = zeros(2,length(N_vec));
        Delta_grad_J_vec = zeros(2,length(N_vec));
        s_vec = zeros(1,length(N_vec));
        t_sim = zeros(1,length(N_vec));
        t_sim_ext = zeros(1,length(N_vec));
        t_J = zeros(1,length(N_vec));
        t_opt = zeros(1,length(N_vec));
        mu_min_vec = zeros(2,length(N_vec));
        val_ind_vec = zeros(1,length(N_vec));
        Delta_mu_vec = zeros(1,length(N_vec));
        %loop over several basis sizes
        for n=1:length(N_vec)
            disp(n)
            N=N_vec(n);
            detailed_data.RB = detailed_data_full.RB(:,1:N);
            reduced_data = gen_reduced_data(model, detailed_data);
            RB = detailed_data.RB;
              M = model_data.df_info.l2_inner_product_matrix;
              reduced_data.J1 = RB'*M*RB;
              reduced_data.J2 = ref_dofs'*M*RB;
              reduced_data.J3 = ref_dofs'*M*ref_dofs;
              %JN = 0.5*(a'*J1*a-2*J2*a+J3);
              %gradient J = a'*gJ1*pmui_a - gJ2 pmui_a
              reduced_data.gJ1 = RB'*M*RB;
              reduced_data.gJ2 = ref_dofs'*M*RB;
            reduced_optimization_functional = @(mu1_mu2) ... 
                  my_reduced_optimization_functional(mu1_mu2, model,...
                                 reduced_data);
  
            
            tic;sim_data = rb_simulation(model, reduced_data);t_s = toc;
            sim_data = rb_reconstruction(model, detailed_data,sim_data);
            
            tic;
            sim_data_ext = extended_rb_simulation(model,reduced_data,detailed_data);
            t_se = toc;
            
            %Delta grad
            Delta_grad_J_vec(:,n) = my_Delta_partial_Js(sim_data_ext,detailed_data,ref_dofs);
            
            
            fun = reduced_optimization_functional;
            tic;
            [mumin,fval,exitflag,output,lambda,grad] = fmincon(fun,mu0,[],[], ...
                [],[],lb,ub,[],options);
            t_opt(n) = toc;
            %epsilon_J = norm(grad);%output.firstorderopt*sqrt(2);
            model = set_mu(model,mumin);
             grad_J = my_reduced_functional_gradient(model,reduced_data);
            epsilon_J = norm(grad_J);
            
            
            % error estimator computation (expensive!!!)
            
            [Delta_mu_vec(n), val_ind_vec(n)] = my_expensive_error_estimator(mumin,model, ...
                detailed_data,reduced_data,...
                detailed_optimization_functional,...
                epsilon_J, ...
                ref_dofs);
            
            
            
            %track data
            Delta_vec(n) = sim_data.Delta;
            Delta_pmui_vec(:,n) = sim_data_ext.Delta_partial_mui_us';
            tic;
            s_vec(n) = J(sim_data.uh.dofs);t_J(n) = toc;
            t_sim(n) = t_s;
            t_sim_ext(n) = t_se;
            mu_min_vec(:,n) = mumin';
        end
            
        t_J
            
        figure(1); clf
        semilogy(N_vec,Delta_vec)
        title('Error estimator for solution')
        xlabel('N')
        
        figure(2);clf
        semilogy(N_vec, Delta_pmui_vec(1,:))
        hold on
        semilogy(N_vec, Delta_pmui_vec(2,:),'r')
        xlabel('N')
        legend('pmu1','pmu2')
        title('Fehlersch?tzer Ableitungen')
        
        figure(3); clf;
        semilogy(N_vec, s_vec)
        title('Output')
        xlabel('N');
        ylabel('s');
        
        figure(4); clf;
        semilogy(N_vec,Delta_grad_J_vec(1,:));
        hold on
        semilogy(N_vec, Delta_grad_J_vec(2,:),'r');
        title('Gradient J Error estimator')
        xlabel('N');
        legend('pmu1','pmu2')
        
        figure(5); clf;
        plot(N_vec,t_sim)
        hold on; plot(N_vec,t_sim_ext,'r')
        title('online simulation times');
        legend('rb sim','rb sim expensive')
        xlabel('N')
        
        
        mu_err = mu_min_vec-repmat(mu_star(1:2)',1,length(N_vec));
        for i=1:length(N_vec)
            mu_err_norm(i) = norm(mu_err(:,i));
        end
        figure(6);clf;
        semilogy(N_vec,mu_err_norm)
        hold on
        semilogy(N_vec,Delta_mu_vec,'r');
        title('Error in optimal parameters')
        legend('true','est')
        xlabel('N');
        
        %Aufteilung in g?ltige und ung?ltige Fehlersch?tzer
        Delta_mu_v = [];
        Delta_mu_iv = [];
        i=1;
        while val_ind_vec(i)>1
            Delta_mu_iv = [Delta_mu_iv,Delta_mu_vec(i)];
            i = i+1;
        end
        Delta_mu_iv = [Delta_mu_iv,Delta_mu_vec(i)];
        for j=i:length(N_vec)
            Delta_mu_v = [Delta_mu_v, Delta_mu_vec(j)];
        end
        
        figure(7); clf;
        semilogy(N_vec,mu_err_norm,'LineWidth',2)
        hold on
        semilogy(N_vec(1:length(Delta_mu_iv)),Delta_mu_iv,'r-.','LineWidth',2);
        semilogy(N_vec(length(Delta_mu_iv):end),Delta_mu_v,'r-','LineWidth',2);
        %title('Error in optimal parameters')
        legend('true error','error est. (invalid)','error est. (valid)')
        xlabel('number of basis vectors N');
        ylabel('error in optimal parameter')
        
         disp('----------------------------------------')
        disp('Simplex optimization (Matlab fmincon):')
        disp('======================================')
        disp(['Detailed optimization time :',num2str(t_det)]);
        disp(['Reduced optimization time  :', num2str(t_opt(end))]);
        disp(' ')
        disp(['Last true error in optimization: ',num2str(mu_err_norm(end)),' and Delta_mu: ',num2str(Delta_mu_vec(end))])
        disp(['Effectivity of error estimator: ',num2str(Delta_mu_vec(end)/mu_err_norm(end))])
        disp(' ')
        disp(['Reduced optimal parameters: ',num2str(mumin(:)')])
        disp(['Number of functional calculations : ',num2str(output.funcCount)])
        disp(['Functional value at optimum: ',num2str(fval)]);
        disp(['Gradient at optimum: ',num2str(norm(grad(:)))])
        
        keyboard;
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% implement and test reduced functional gradient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    case 10
        
        load('gamm2013_detailed_data');
        mu1_2 = [0.5;0.5];
        mu = mu1_2;
        model = set_mu(model,mu);
        tic;
        grad_J = my_reduced_functional_gradient(model, reduced_data);
        t = toc
        %check numerically with fd
        h = 0.0001;
        hv = eye(2,2).*h;
        gJn=zeros(2,1);
        J0 = reduced_optimization_functional(mu);
        for i=1:2
            gJn(i) = (detailed_optimization_functional(mu+hv(:,i))-J0);
            gJn(i) = gJn(i)/h;
        end
        
        disp(grad_J)
        disp(gJn)
        
        %detailed gradient
        Jdet = my_detailed_functional_gradient(model, model_data);
        disp(Jdet);
        
        
        keyboard;
        
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
 % gradient optimization
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
    case 11
        
        load('gamm2013_detailed_data');
        model.optimization = [];
        model.optimization.init_params = [0.5,0.5];
        model.optimization.tol = 1e-7;
        model.optimization.stepsize_rule = 'armijo';
        model.optimization.stepsize_fix = 1;
        model.optimization.get_Jacobian = @my_reduced_functional_gradient;
        model.optimization.objective_function = @(model,reduced_data)my_reduced_optimization_functional([],model,reduced_data);
        model.optimization.derivatives_available = 0;
        model.optimization.params_to_optimize = [1,1,0,0];
        model.optimization.opt_mode = 'reduced';
        model.optimization.initial_stepsize = 1000;
        model.optimization.sigma = 0.6;
        
        model.compute_err_estimator = 0;
        
        model.expensive_err_est = @(mu,model,grad_J)my_expensive_error_estimator(mu,model, ...
                                                  detailed_data,reduced_data,...
                                                  detailed_optimization_functional,...
                                                  grad_J, ...
                                                  ref_dofs);
        
        %reduced
        tic;
        [opt_data_red, model]=my_gradient_opt(model,reduced_data);
        t_optred = toc
        
        model.compute_err_estimator = 1;
        [opt_data_red, model]=my_gradient_opt(model,reduced_data);
        
        
        %detailed
        model.optimization.get_Jacobian = @my_detailed_functional_gradient;
        model.optimization.objective_function = @(model,model_data)my_detailed_optimization_functional([],J,model,model_data);
        tic;
        [opt_data_det] = my_gradient_opt(model, model_data);
        t_optdet = toc
        
        
        opt_data_red.t_opt = t_optred;
        opt_data_det.t_opt = t_optdet;
        
        save('Optimization_data_gamm.mat','opt_data_red','opt_data_det');
        
        keyboard;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% plot figures for GAMM proceedings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%first conduct cases 3 and 5

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
% plot functional landscape
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 100
        
        exp_params.nmu = 50; %wieviele Gitterpunkte
        
        gamm2013_exp(6,exp_params);
        
        keyboard;
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% plot reference solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 101
        
        load('gamm2013_detailed_data');
        
        model = model.set_mu(model,[2,1]);
        sim_data = detailed_simulation(model,model_data);
        plot_sim_data(model,model_data,sim_data,plot_params);
        
        
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot gradient steps and error estimator decay
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case 102
  
        %gamm2013_exp(11);
        
        mu_ref = [2;1];
        load('gamm2013_detailed_data.mat')
        load('Optimization_data_gamm.mat');
        
        %sortieren in g?ltige und ung?ltige Fehlersch?tzerwerte
        Delta_mu_vec = opt_data_red.Delta_mu_vec;
        norm_grad = opt_data_red.norm_grad_vec;
        
        ind_v = find(opt_data_red.break_value_vec<1);
        ind_iv = find(opt_data_red.break_value_vec>=1);
        Delta_mu_iv = Delta_mu_vec(ind_iv(1):ind_iv(end)+1);
        norm_grad_iv = norm_grad(ind_iv(1):ind_iv(end)+1);
        Delta_mu_v = Delta_mu_vec(ind_v);
        norm_grad_v = norm_grad(ind_v);
        
        %sort
        [norm_grad,idx1] = sort(norm_grad,'descend');
        [norm_grad_iv,idx2] = sort(norm_grad_iv,'descend');
        Delta_mu_iv = Delta_mu_iv(idx2);
        [norm_grad_v,idx3] = sort(norm_grad_v,'descend');
        Delta_mu_v = Delta_mu_v(idx3);
        
        %true error
        dist = opt_data_red.parameter_sets-repmat(mu_ref,1,size(opt_data_red.parameter_sets,2));
        err = zeros(1,size(dist,2));
        for i=1:size(dist,2)
            err(i) = norm(dist(:,i));
        end
        err = err(idx1);
        
        figure
        loglog(norm_grad,err,'LineWidth',2);
        hold on
        loglog(norm_grad_iv,Delta_mu_iv,'r-.','LineWidth',2);
        loglog(norm_grad_v,Delta_mu_v,'r','LineWidth',2);
        set(gca, 'XDir', 'reverse')
        xlabel('gradient norm')
        ylabel('error in parameters')
        legend('true error','error est. (invalid)','error est. (valid)')
        
        %show some data on screen
        disp('----------------------------------------')
        disp('Gradient optimization:')
        disp('=======================')
        disp(['Detailed optimization time :',num2str(opt_data_det.t_opt)]);
        disp(['Reduced optimization time  :', num2str(opt_data_red.t_opt)]);
        disp(' ')
        disp(['Last true error in optimization: ',num2str(err(end)),' and Delta_mu: ',num2str(Delta_mu_vec(end))])
        disp(['Effectivity of error estimator: ',num2str(Delta_mu_vec(end)/err(end))])
        disp(' ')
        disp(['Reduced optimal parameters: ',num2str(opt_data_red.optimal_params(:)')])
        disp(['Number of functional calculations (armijo): ',num2str(opt_data_red.nr_fct_calls)])
        disp(['Number o gradient calculations : ',num2str(opt_data_red.nr_grad_calc)])
        disp(' ')
        disp(['Functional value at optimum: ',num2str(my_reduced_optimization_functional(opt_data_red.optimal_params,model,reduced_data))]);
        disp(['Functional value at optimum(det): ',num2str(my_detailed_optimization_functional(opt_data_red.optimal_params,J,model,model_data))]);
        disp(['Gradient at optimum: ',num2str(opt_data_red.norm_grad_f_x)])
        keyboard;
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5        
%        plot_variation of number of basis vectors
%   (simplex optimization)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55

    case 103
        
        gamm2013_exp(9);
        
  
 otherwise
  error('step number unknown')
end;

% functional J, here l2 error functional J(u) := \int_Omega |u-u_star|^2
function res = J_functional(dofs, l2_inner_product_matrix, ref_dofs)
err =  dofs-ref_dofs;
res = 0.5* err'*l2_inner_product_matrix* err;

% Implementation of the detailed optimization functional
function [res,dummy] = my_detailed_optimization_functional(mu1_mu2, J, model, model_data)
dummy = [];
if ~isempty(mu1_mu2)
    model = model.set_mu(model,[mu1_mu2(:)]);
end
sim_data = detailed_simulation(model,model_data);
fprintf('.');
res = J(sim_data.uh.dofs);

% reduced optimization functional: non offline/online decomposed!
function [res,dummy] = my_reduced_optimization_functional(mu1_mu2, model, ...
                reduced_data)
    dummy = [];
if ~isempty(mu1_mu2)                    
    model = model.set_mu(model,mu1_mu2(:));
end
sim_data = rb_simulation(model,reduced_data);
uN = sim_data.uN;
res = 0.5*(uN'*reduced_data.J1*uN-2*reduced_data.J2*uN+reduced_data.J3);%J(detailed_data.RB*sim_data.uN);
fprintf('.');

function [Delta_mu, val_indicator] = my_expensive_error_estimator(mumin,model, ...
                                                  detailed_data,reduced_data,...
                                                  optimization_functional,...
                                                  epsilon_J, ...
                                                  ref_dofs);
% compute parameter error estimator Delta_mu and 
% validity indicator val_indicator = 4 gamma^2 L epsilon
% should be less than 1 in order the bound to be valid.

% gamma: Induced norm of inverse detailed hessian
%  => approximate Hessian by finite difference +-h
h = 1e-6;
H = FD_Hessian(mumin,optimization_functional, h)
gamma = norm(inv(H))

% Compute Delta_Grad_J by extended RB-simulation (for sensitivities)
sim_data = extended_rb_simulation(model,reduced_data,detailed_data);
Delta_partial_Js = my_Delta_partial_Js(sim_data,detailed_data,ref_dofs);
% sum up all partial derivative estimators
Delta_Grad_J = norm(Delta_partial_Js);

epsilon = Delta_Grad_J + epsilon_J;

Delta_mu = 2*gamma*epsilon;

% L: maximum Hessian Lipschitz constant in 2*gamma*epsilon  ball around mumin
% for now: simple maximum search over few points
L = 0;
mus = [0,1; 1,0; sqrt(2)/2,sqrt(2)/2];
nmus = size(mus,1);
for n = 1:nmus
  H2 = FD_Hessian(mumin(:)'+mus(n,:)*Delta_mu,optimization_functional, h);
  L2 = norm(H-H2)/Delta_mu; % Delta_mu is difference between mus
  if L2>L
    L = L2;
  end;
end;

val_indicator = 4 * gamma^2 * L * epsilon;

% compute approximate Hessian by 6 point evaluations
%            f12         f22
%   f01      f11 = f(x)  f21 = f(x+[h,0])
%            f10  
function H = FD_Hessian(x,f, h);
x = x(:);
H = zeros(2,2);
f11 = f(x); 
f01 = f(x-[h;0]); 
f21 = f(x+[h;0]); 
f10 = f(x-[0;h]); 
f12 = f(x+[0;h]); 
f22 = f(x+[h;h]); 
% d1d1f:
H(1,1) = (f21-2*f11+f01)/h^2;
% d2d2f:
H(2,2) = (f12-2*f11+f10)/h^2;
% d1d2f: (a bit inaccurate...)
H(1,2)= ((f22-f21)/h -(f12-f11)/h)/h;
H(2,1)=H(1,2);

function grad_J = my_detailed_functional_gradient(model, model_data)

grad_J = zeros(2,1);

sim_data = extended_detailed_simulation(model, model_data);

M = model_data.df_info.l2_inner_product_matrix;

for i = 1:2
    grad_J(i) = (sim_data.uh.dofs-model_data.ref_dofs)'*M*sim_data.u_der(:,i);
end


% return gradient of the functional J (reduced)
function grad_J = my_reduced_functional_gradient(model, reduced_data)
%perform a reduced simulation
grad_J = zeros(2,1);

model.compute_err_estimator = 0;
sim_data = extended_rb_simulation(model, reduced_data,[]);

%M = detailed_data.df_info.l2_inner_product_matrix;
%RB = detailed_data.RB;
%for i=1:2
%    grad_J(i) = sim_data.uN'*RB'*M*RB*sim_data.partial_mui_uNs(:,i) ...
%                - ref_dofs'*M*RB*sim_data.partial_mui_uNs(:,i);
%end
for i=1:2
    grad_J(i) = sim_data.uN'*reduced_data.gJ1*sim_data.partial_mui_uNs(:,i) ...
                - reduced_data.gJ2*sim_data.partial_mui_uNs(:,i);
end


% error estimators for gradient of J for quadratic functional
% |partial_mui J(u) - partial_mui J(uN)|
function Delta_partial_Js = my_Delta_partial_Js(sim_data,detailed_data,ref_dofs)
Delta_partial_Js = [0,0]; % using Delta_partial_us
err = detailed_data.RB * sim_data.uN - ref_dofs;
norm_uN_minus_ud = sqrt(err' * ...
                        detailed_data.df_info.h10_inner_product_matrix ...
                        * err);
for i = 1:2
  partial_mui_uN_dofs = detailed_data.RB*sim_data.partial_mui_uNs(:,i);
  norm_partial_uNi = sqrt(partial_mui_uN_dofs' * ...
                          detailed_data.df_info.h10_inner_product_matrix * ...
                          partial_mui_uN_dofs);                 
  Delta_partial_Js(i) = sim_data.Delta* (sim_data.Delta_partial_mui_us(i) ...
                                         + norm_partial_uNi) ...
      + norm_uN_minus_ud*sim_data.Delta_partial_mui_us(i);
end;

% RB-simulation including sensitivities and error estimators
% but expensive computation by high-dimensional operations!
function sim_data = extended_rb_simulation(model,reduced_data,detailed_data)

if ~isfield(model,'compute_err_estimator')
    model.compute_err_estimator = 1;
end
sim_data = rb_simulation(model,reduced_data);
sim_data.Delta_s = []; %not available
sim_data.s = []; %not available here
% solve sensitivity problems:

%old_mode = model.decomp_mode;
model.decomp_mode = 2; % coefficients
[A_coeff,f_coeff] = ...
    model.operators(model,[]);
AN = lincomb_sequence(reduced_data.AN_comp, A_coeff);
fN = lincomb_sequence(reduced_data.fN_comp, f_coeff);
uN = AN\fN;

% solve sensitivity problems, coefficient derivatives by finite difference:

% matrix of reduced sensitivity coefficients, each column one reduced solution:
sim_data.partial_mui_uNs = zeros(length(uN),2);
% vector of error estimators:
sim_data.Delta_partial_mui_us = zeros(1,2);
mu = model.get_mu(model);
h = 1e-6;
E = eye(2);
if model.compute_err_estimator
    K = model.get_inner_product_matrix(detailed_data);
    model.decomp_mode = 1;
    [A_comp,f_comp] = ...
    model.operators(model,detailed_data);
end
model.decomp_mode = 2; % coefficients
alpha = model.coercivity_alpha(model);


for i = 1:2
  % 1. solve for sensitivity solutions:
  model = model.set_mu(model,mu+h*E(:,i));
  [A_coeff2,f_coeff2] = ...
      model.operators(model,[]);
  partial_mui_A_coeff = (A_coeff2-A_coeff)/h;
  partial_mui_f_coeff = (f_coeff2-f_coeff)/h;  
  partial_mui_AN = lincomb_sequence(reduced_data.AN_comp, partial_mui_A_coeff);
  partial_mui_fN = lincomb_sequence(reduced_data.fN_comp, partial_mui_f_coeff);
  partial_mui_uN = AN \ (partial_mui_fN - partial_mui_AN*uN);
  sim_data.partial_mui_uNs(:,i) = partial_mui_uN; 

  if model.compute_err_estimator
      %2. compute residuals and error estimators
      % Delta_partial_mui_u = gamma_partial_mui_a/alpha * Delta_u + |r|/alpha
      % with (r,v ) = partial_mui_i f(v)-a(partial_mui_uN,v) -partial_mui_a(u_N,v) 
      A = lincomb_sequence(A_comp, A_coeff);
      partial_mui_A = lincomb_sequence(A_comp, partial_mui_A_coeff);
      partial_mui_f = lincomb_sequence(f_comp, partial_mui_f_coeff);

      res_partial_mui_A = partial_mui_A * (detailed_data.RB * uN);
      res_A = A * (detailed_data.RB * partial_mui_uN);
      res = partial_mui_f - res_partial_mui_A - res_A;
      %%% residuum functional is res * v 
      %%% riesz representant: v_r' K v = (v_r,v) = res*v
      %%% so res = K * v_r;
      v_r = K \ res;
      %%% res_norm_sqr = (v_r,v_r) = v_r' K v_r = v_r' * res;
      res_norm_sqr = v_r' * res;  
      res_norm = sqrt(max(res_norm_sqr,0));

      % continuity constant of bilinear form partial_mui_a:
      % keyboard;
      gamma_partial_mui_a = get_continuity_constant(partial_mui_A,K);

      sim_data.Delta_partial_mui_us(i) = ...
          gamma_partial_mui_a / alpha * sim_data.Delta ...
          + res_norm / alpha;
  end

end;

%detailed simulation returning sensitivities
function sim_data = extended_detailed_simulation(model, model_data)

sim_data = detailed_simulation(model,model_data);

% solve sensitivity problems:
model.decomp_mode = 2; % coefficients
[A_coeff,f_coeff] = ...
    model.operators(model,[]);

% solve sensitivity problems, coefficient derivatives by finite difference:

% matrix of reduced sensitivity coefficients, each column one reduced solution:
u = sim_data.uh.dofs;
u_der = zeros(length(u),2);


mu = model.get_mu(model);
h = 1e-6;
E = eye(2);
%K = model.get_inner_product_matrix(detailed_data);
model.decomp_mode = 1;
[A_comp,f_comp] = ...
 model.operators(model,model_data);

model.decomp_mode = 2; % coefficients



for i = 1:2
  % 1. solve for sensitivity solutions:
  model = model.set_mu(model,mu+h*E(:,i));
  [A_coeff2,f_coeff2] = ...
      model.operators(model,[]);
  partial_mui_A_coeff = (A_coeff2-A_coeff)/h;
  partial_mui_f_coeff = (f_coeff2-f_coeff)/h;  
  partial_mui_A = lincomb_sequence(A_comp, partial_mui_A_coeff);
  partial_mui_f = lincomb_sequence(f_comp, partial_mui_f_coeff);
  A = lincomb_sequence(A_comp, A_coeff);
  u_der(:,i) = A \ (partial_mui_f - partial_mui_A*u);
  sim_data.u_der = u_der;

end;


% function computing the continuity constant of a sparse matrix A
% with respect to inner product matrix K
% formally: c = largest singular value of matrix (K^-0.5 * A * K^-0.5)
% where K^-0.5 are inverse cholesky factors. These are not
% explicitly computed as full, but linear systems are solved.
% largest singular value of M = inv(C')*A*inv(C) is computed by
% largest eigenvalue of (M'* M)
function c = get_continuity_constant(A,K)
c=1;
%C = chol(K); % => K = C'*C
%f = @(x) my_matrix_prod(x,A,C);
%lambda_max = eigs(f,size(A,1),1);
%c = sqrt(lambda_max);

function res = my_matrix_prod(x,A,C)
y = C\x; %y = invC*x;
y2 = (C') \ (A*y);% y2 = inv(CT)*A*inv(C)*x
y3 = C\y2;
y4 = (C') \ (A'*y3); % y4 = inv(CT)*AT*inv(C)*inv(CT)*A*inv(C)*x
res = y4;

function cmaplog = create_log_colormap()

%usual colormap
cmap = colormap;

len = size(cmap,1);

x = (len-1)/log10(len-1);
cmaplog = zeros(size(cmap));

for n=1:len
    ind = x*log10(n)+1;
    a = floor(ind);
    b = ceil(ind);
    if b>len
        b = len;
    end
    m = ind-a;
    v1 = cmap(a,:);
    v2 = cmap(b,:);
    
    cmaplog(n,:) = (1-m) * v1 + m*v2;
end

cmap = cmaplog;


function model = my_pfusch_set_mu(model,mu1_2,mu3,mu4)
model.mus = [mu1_2(:);mu3;mu4];

%function [grad_f_x,Delta_grad,output,data] = my_return_jacobian(model,reduced_data)
%grad_f_x = my_reduced_functional_gradient(model, reduced_data);
%Delta_grad = [];
%output = [];
%data = [];
%if model.compute_err_estimator
%    Delta_grad = 
%end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%gradient optimization

function [opt_data, model]=my_gradient_opt(model,reduced_data)
%function [x_opt,nr_fct_total,nr_grad_total]=gradient_opt(model, model_data, reduced_data)
%
% Function performing an optimization either reduced oder detailed using gradient method
% with either Armijo-stepsize rule, Wolfe-Powell-stepsize rule or
% dichotomie-algorithm
%
%Required fields of model:
% model.optimization.init_params: initial value
% model.optimization.tol: tolerance how close to 0 shall the norm be?
% model.optimization.get_Jacobian: gradient function
% model.optimization.objective_function: not needed in gradient_opt but
% needed in armijo_stepsize.m and wolfe_powell_stepsize.m
% model.optimization.stepsize_rule: either 'armijo' or 'wolfepowell' or
% 'dichotomie'
%
%Generated fields of opt_data:
% opt_data.x_optimal: the optimal value
% opt_data.nr_fct_calls: Number of function calls
% opt_data.nr_grad_calc: Number of needed calculations of the gradient
%
%
% Oliver Zeeb 25.05.2010%



if(model.verbose>=8)
    disp('entered gradient_opt')
end

%adding some fields to the model, if not existent
if ~isfield(model.optimization,'tol')
    model.optimization.tol=1e-3;
end
if ~isfield(model.optimization,'stepsize_rule')
    model.optimization.stepsize_rule='quotient';
end
if ~isfield(model.optimization,'min_or_max')
    model.optimization.min_or_max = 'min';
end

if model.optimization.derivatives_available
    model.compute_derivative_info = 1;
end

%setting signum for maximizing or minimizing,
% i.e. changes f(x) to -f(x) if function f(x) shall be maximized
if model.optimization.min_or_max == 'min'
    min_or_max=1;
elseif model.optimization.min_or_max == 'max'
    min_or_max=-1;
else disp('maximize or minimize? model.optimization.min_or_max not properly defined')
    return
end
    
%for statistics
nr_fct = 0;     %counter for number of function calls
nr_grad = 0;    %counter for number of gradient calculations
nRB= []; %basis sizes used for reduced simulations
nRB_der = []; %basis sizes used for reduced derivative simulations

%save current mu values of the model

mu_original   = get_mu(model);


%set the mus in the model to the mus given in init_params
%model.(model.mu_names{k})=model.optimization.init_params(k);
model = set_mu(model, model.optimization.init_params);


%get initial parameters that should be optzimized
x = get_mu(model);
%get the according boundaries 
[lower_bound, upper_bound] = get_bound_to_optimize(model);




tol         = model.optimization.tol;    
PM          = x;        %for opt_data: array with all the mu_i
output      = [];       %array containing the f(mu_k)
max_iter    = 100;   %maximal number of iterations of the gradient method
opti_iter   = 0;    %counter
%output_array={};

if strcmp(model.optimization.opt_mode,'detailed')    

else  %here: reduced optimization

    
[grad_f_x] = model.optimization.get_Jacobian(model,reduced_data);
grad_f_x    = grad_f_x*min_or_max; nr_grad = nr_grad+1;
norm_grad_f_x = norm(grad_f_x);
norm_grad_vec = norm_grad_f_x;
if model.compute_err_estimator
    [Delta_mu_vec,break_value_vec] =  model.expensive_err_est(x,model,norm_grad_f_x);
end



continue_optimization = 1;

% START OF OPTIMIZATION
while (continue_optimization)&& (opti_iter <= max_iter)
    d = -grad_f_x;
    
    %calculation of the stepsize t:
    switch model.optimization.stepsize_rule
        case {'fixed_stepsize'}
            t = model.optimization.stepsize_fix;
        case {'armijo'} % Armijo Regel
            [model,t,nr_fct_stepsize,output, data] = stepsize_armijo(model, reduced_data, output, x, d);
            nr_fct = nr_fct+nr_fct_stepsize;
            nRB = [nRB,data.nRB];
           
        case {'dichotomie'} 
            [model,t,nr_fct_stepsize,output] = stepsize_dichotomie(model, model_data, output, x, d);
            nr_fct = nr_fct+nr_fct_stepsize;     
            
        case {'exponential'}
            if ~exist('exp_step')
                exp_step=0;
            end
            [model,t,nr_fct_stepsize,output] = stepsize_exponential(model, model_data, output, x, d, exp_step);
            exp_step=exp_step+1;
            nr_fct = nr_fct+nr_fct_stepsize;
            
        case {'quotient'}
            if ~exist('quot_step')
                quot_step=1;
            end
            [model,t,nr_fct_stepsize] = stepsize_quotient(model, x, d, quot_step);
            quot_step=quot_step+1;
            nr_fct = nr_fct+nr_fct_stepsize;
            %C_alpha = t; %update of the maximum stepsize in this step for error estimation
            %[C_L,model] = model.get_lipschitz_constant(model); %update Lipschitz constant
            
        %WOLFE POWELL NOT WORKING AT THE MOMENT!!!
        case {'wolfepowell'} 
            [t, nr_fct_stepsize, nr_grad_stepsize] = stepsize_wolfe_powell(model, x, d);
            nr_fct = nr_fct+nr_fct_stepsize;
            nr_grad = nr_grad+nr_grad_stepsize;

        otherwise
            fprintf('unknown  stepsize rule')
            opt_data=[];
            return
    end
    
    opti_iter = opti_iter+1;
    x = x+t*d;
    if model.verbose >3
        disp(['Next parameter in reduced gradient optimization: ',num2str(x(:)')]);
        disp('------------------------------------')
    end
    %check, whether the new x is allowed by the according mu_range or whether it exceeds the
    %boundaries. If so: set this component to the boundary.
    for k=1:length(x)
        if x(k) > upper_bound(k)
            x(k) = upper_bound(k);
        elseif x(k) < lower_bound(k)
            x(k) = lower_bound(k);
        end
    end
        
    PM = [PM, x(:)];
    model = set_mu(model,x); % write the actual parameter setting in the model
    
    
    %Calculate gradient at the new parameter point:
    [grad_f_x] = model.optimization.get_Jacobian(model,reduced_data);%now calculate the gradient with the new parameter setting
    norm_grad_f_x = norm(grad_f_x);
    norm_grad_vec = [norm_grad_vec,norm_grad_f_x];
    %output=output_array{1};
    grad_f_x=grad_f_x*min_or_max;  
    nr_grad = nr_grad+1; 
    
    %decide to continue or not
    %break_value = 2*model.gamma_H^2*model.L_H*(norm_grad_f_x+norm(Delta_grad(:,end)));
    if model.compute_err_estimator
        [Delta_mu,break_value] =  model.expensive_err_est(x,model,norm_grad_f_x);
        break_value_vec = [break_value_vec, break_value];
        Delta_mu_vec = [Delta_mu_vec, Delta_mu];
        disp(['Actual gradient norm :',num2str(norm_grad_f_x(:)')]);
        disp(['Value of breaking criteria: ', num2str(break_value)]);
    end
    %if (break_value<=1)&&(norm_grad_f_x<tol)
    if (norm_grad_f_x<tol)
        continue_optimization = 0;
    end
   
end
% END OF OPTIMIZATION


%if function was maximized instead of minimized, the values must be
%multiplied by -1:
%output=output*min_or_max;

%setting opt_data
opt_data = [];
opt_data.optimal_params = x;
opt_data.nr_fct_calls   = nr_fct;
opt_data.nr_grad_calc   = nr_grad;
opt_data.parameter_sets = PM;
%opt_data.output = output;
%opt_data.max_output = max(output);
opt_data.max_output_paramsets = PM(end,:);
 opt_data.norm_grad_f_x = norm_grad_f_x;
 opt_data.norm_grad_vec = norm_grad_vec;
if model.compute_err_estimator
    opt_data.nr_opti_iter=opti_iter;
    opt_data.Delta_mu = Delta_mu;
    opt_data.break_value = break_value;
    %opt_data.nRB = nRB;
    %opt_data.nRB_der = nRB_der;
    %opt_data.used_detailed_calculations = used_detailed_calculations;
    %opt_data.gamma = g_H;
    %opt_data.Lipschitz_H = L_H;
    opt_data.break_value_vec = break_value_vec;
    opt_data.Delta_mu_vec = Delta_mu_vec;
end
end






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Improvement possibilities:
%  - offline-online decomposition of reduced optimization
%    functional
%  - extended RB simulation: offline online decomposition
%  - Make RB basisgeneration more stable (error estimator below 1e-6)?
%    currently about 20 basisvectors, then oscillations start
%  - determine number of function evaluations for detailed/reduced
%  - use gradient scheme for optimization and view improvements and
%    validity of error estimator
%  - nicer Hessian

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Findings:
% - full optimization: 35 sec. Reduced: 0.35 sec factor 100 acceleration