calculate_gamma_const.m

function [gamma_H,model] = calculate_gamma_const(model, varargin)
%function gamma_H = calculate_gamma_const(model, varargin)
%
% function caluclates the gamma constant using detailed simulations, which is needed for the
% non-iterative error estimator in optimization.
%
% Optional inputs:
% model_data
% mu
%

% Markus Dihlmann 03.01.2013

compute_derivative_info_old = model.compute_derivative_info;

%see if there's model_data
if nargin>1
    model_data = varargin{1};
    if isempty(model_data)
        model_data = gen_model_data(model);
    end
else
    model_data = gen_model_data(model);
end

%see if there's a mu given
if nargin>2
    mu = varargin{2};
    if ~isempty(mu)
        model = set_mu(model, mu);
    end
end

%claculate the Hessian
disp('performing detailed simulations to calculate the gamma constant')
H = lin_evol_opt_fd_Hessian(model, model_data,[]);

model = setfield_in_model_and_base_model(model,'compute_derivative_info', 0);
HJ = zeros(size(H));
for k=1:size(H,1)
    for l=1:size(H,2)
        sim_data_dummy.U = H{k,l};
        sim_data_dummy = model.get_output(model, sim_data_dummy, model_data);
        HJ(k,l) = sim_data_dummy.y(end);
    end
end

%reset model
model= setfield_in_model_and_base_model(model,'compute_derivative_info', compute_derivative_info_old);

%return Gamma
gamma_H = norm(inv(HJ));