porsche_get_Hessian_J.m

function [hessian, Delta_hessian]=porsche_get_Hessian_J(model, data)
%function [hessian, Delta_hessian]=porsche_get_Hessian_J(model, data)
%
% This function calculates the Hessian (w.r.t. the mus given in the model)
% of the objective-function via finite differences either in
% detailed or reduced mode, according to model.optimization.opt_mode. 
% The corresponding data (model_data or reduced_data) must be passed to
% this function.
%
% input:
% - model
% - data: model_data oder reduced_data
% 
% needed fields in model:
% model.optimization.opt_mode = 'reduced' or 'detailed'
%
% Oliver Zeeb, 31.05.2011

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

Delta_hessian = 0;
%disp('in porsche_get_Hessian_J: Fehlerschätzer noch nicht implementiert!');

h_range = 1000; % for constructing the difference quotient: how small is the variaton of the parameter
%note: h_range = 1000 gives quite good results when detailed and reduced
%calculation are compared. If h_range is too hig (e.g. 10.000), the results
%in the reduced case become worse!


model_original=model;
k=0;

%new: with flag model.optimization.opt_mode
%calculate derivative via small differences
func_mu = model.optimization.objective_function(model,data); %value of the objective function for the original parameter mu, calculated some lines above
hessian=zeros(length(model.mu_names)); 
for i=1:length(model.mu_names)
    if model.optimization.params_to_optimize(i) == 1
        h_i=(model.mu_ranges{i}(2)-model.mu_ranges{i}(1))/h_range;
        for j=i:length(model.mu_names)
            h_j=(model.mu_ranges{j}(2)-model.mu_ranges{j}(1))/h_range;
            if i==j %use other formula for d_ii than for d_ij
                %get the nedded values for the
                %hessian-finite-difference-formula
                model=model_original;
                model.(model.mu_names{i}) = model.(model.mu_names{i}) + 2*h_i;
                func_mu_p2hi = model.optimization.objective_function(model,data); %objective_function in mu + 2*h_i
                
                model=model_original;
                model.(model.mu_names{i}) = model.(model.mu_names{i}) + h_i;
                func_mu_phi = model.optimization.objective_function(model,data); %objective_function in mu + h_i
                
                model=model_original;
                model.(model.mu_names{i}) = model.(model.mu_names{i}) - h_i;
                func_mu_mhi = model.optimization.objective_function(model,data); %objective_function in mu - h_i
                
                model=model_original;
                model.(model.mu_names{i}) = model.(model.mu_names{i}) - 2*h_i;
                func_mu_m2hi = model.optimization.objective_function(model,data); %objective_function in mu - 2*h_i
                
                model=model_original;
                partial_der = (- func_mu_p2hi + 16*func_mu_phi - 30*func_mu + 16*func_mu_mhi - func_mu_m2hi) / (12*h_i^2);
                
                hessian(i,i) = partial_der;
                
                if(model.verbose>=8)
                    disp('entered porsche_get_Jacobian')
                    k=k+1;
                    disp(['values calculated: ' num2str(k) ' of ' num2str((length(model.mu_names)*(length(model.mu_names)+1))/2) ]);
                end
                
            else % i neq j --> calculate d_ij
                %get the nedded values for the
                %hessian-finite-difference-formula
                model=model_original;
                model.(model.mu_names{i}) = model.(model.mu_names{i}) + h_i;
                model.(model.mu_names{j}) = model.(model.mu_names{j}) + h_j;
                func_mu_phi_phj = model.optimization.objective_function(model,data); %objective_function in mu + h_i + h_j
                
                model=model_original;
                model.(model.mu_names{i}) = model.(model.mu_names{i}) + h_i;
                model.(model.mu_names{j}) = model.(model.mu_names{j}) - h_j;
                func_mu_phi_mhj = model.optimization.objective_function(model,data); %objective_function in mu + h_i - h_j
                
                model=model_original;
                model.(model.mu_names{i}) = model.(model.mu_names{i}) - h_i;
                model.(model.mu_names{j}) = model.(model.mu_names{j}) + h_j;
                func_mu_mhi_phj = model.optimization.objective_function(model,data); %objective_function in mu - h_i + h_j
                
                model=model_original;
                model.(model.mu_names{i}) = model.(model.mu_names{i}) - h_i;
                model.(model.mu_names{j}) = model.(model.mu_names{j}) - h_j;
                func_mu_mhi_mhj = model.optimization.objective_function(model,data); %objective_function in mu - h_i - h_j
                
                model=model_original; %setting the mus in the model back to their original values
                partial_der = (func_mu_phi_phj - func_mu_phi_mhj - func_mu_mhi_phj + func_mu_mhi_mhj)/(4*h_i*h_j); %calculating the partial derivative
                
                hessian(i,j) = partial_der;
                hessian(j,i) = partial_der; %symmetric Hessian --> H(i,j) = H(j,i)
                
                if(model.verbose>=8)
                    disp('entered porsche_get_Jacobian')
                    k=k+1;
                    disp(['values calculated: ' num2str(k) ' of ' num2str((length(model.mu_names)*(length(model.mu_names)+1))/2) ]);
                end
            end
        end %for j=i:length(model.mu_names)
    end %if model.optimization.params_to_optimize(i) == 1
end %for i=1:length(model.mu_names)


%old: without flag model.optimization.opt_mode
% %detailed case
% if strcmp(inputname(2),'model_data')
%     model_data = varargin{1};
%     %calculate derivative via small differences
%     func_mu = model.optimization.objective_function(model,model_data); %value of the objective function for the original parameter mu, calculated some lines above
%     hessian=zeros(length(model.mu_names)); %%%% ACHTUNG!!! Oli, 30.05.11
%     disp('in porsche_get_hessian_J: aufpassen, falls nicht alle params_to_optimimze = 1');
%     for i=1:length(model.mu_names)
%         if model.optimization.params_to_optimize(i) == 1
%             h_i=(model.mu_ranges{i}(2)-model.mu_ranges{i}(1))/h_range;
%             for j=i:length(model.mu_names)
%                 h_j=(model.mu_ranges{j}(2)-model.mu_ranges{j}(1))/h_range;
%                 if i==j %use other formula for d_ii than for d_ij
%                     %get the nedded values for the
%                     %hessian-finite-difference-formula
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) + 2*h_i;
%                     func_mu_p2hi = model.optimization.objective_function(model,model_data); %objective_function in mu + 2*h_i
%                     
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) + h_i;
%                     func_mu_phi = model.optimization.objective_function(model,model_data); %objective_function in mu + h_i
%                     
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) - h_i;
%                     func_mu_mhi = model.optimization.objective_function(model,model_data); %objective_function in mu - h_i
%                     
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) - 2*h_i;
%                     func_mu_m2hi = model.optimization.objective_function(model,model_data); %objective_function in mu - 2*h_i
%                     
%                     model=model_original;
%                     partial_der = (- func_mu_p2hi + 16*func_mu_phi - 30*func_mu + 16*func_mu_mhi - func_mu_m2hi) / (12*h_i^2);
%                     
%                     hessian(i,i) = partial_der;
%                     
%                     if(model.verbose>=8)
%                         disp('entered porsche_get_Jacobian')
%                         k=k+1;
%                         disp(['values calculated: ' num2str(k) ' of ' num2str((length(model.mu_names)*(length(model.mu_names)+1))/2) ]);
%                     end
%                     
%                 else % i neq j --> calculate d_ij
%                     %get the nedded values for the
%                     %hessian-finite-difference-formula
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) + h_i;
%                     model.(model.mu_names{j}) = model.(model.mu_names{j}) + h_j;
%                     func_mu_phi_phj = model.optimization.objective_function(model,model_data); %objective_function in mu + h_i + h_j
%                     
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) + h_i;
%                     model.(model.mu_names{j}) = model.(model.mu_names{j}) - h_j;
%                     func_mu_phi_mhj = model.optimization.objective_function(model,model_data); %objective_function in mu + h_i - h_j
%                     
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) - h_i;
%                     model.(model.mu_names{j}) = model.(model.mu_names{j}) + h_j;
%                     func_mu_mhi_phj = model.optimization.objective_function(model,model_data); %objective_function in mu - h_i + h_j
%                     
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) - h_i;
%                     model.(model.mu_names{j}) = model.(model.mu_names{j}) - h_j;
%                     func_mu_mhi_mhj = model.optimization.objective_function(model,model_data); %objective_function in mu - h_i - h_j
%                     
%                     model=model_original; %setting the mus in the model back to their original values
%                     partial_der = (func_mu_phi_phj - func_mu_phi_mhj - func_mu_mhi_phj + func_mu_mhi_mhj)/(4*h_i*h_j); %calculating the partial derivative
%                     
%                     hessian(i,j) = partial_der;
%                     hessian(j,i) = partial_der; %symmetric Hessian --> H(i,j) = H(j,i)
%                     
%                     if(model.verbose>=8)
%                         disp('entered porsche_get_Jacobian')
%                         k=k+1;
%                         disp(['values calculated: ' num2str(k) ' of ' num2str((length(model.mu_names)*(length(model.mu_names)+1))/2) ]);
%                     end
%                 end
%             end %for j=i:length(model.mu_names)
%         end %if model.optimization.params_to_optimize(i) == 1
%     end %for i=1:length(model.mu_names)
%     
%     
% %reduced case
% elseif strcmp(inputname(2), 'reduced_data')
%     reduced_data = varargin{1};
%     %calculate derivative via small differences
%     func_mu = model.optimization.objective_function(model,reduced_data); %value of the objective function for the original parameter mu, calculated some lines above
%     hessian=zeros(length(model.mu_names)); %%%% ACHTUNG!!! Oli, 30.05.11
%     disp('in porsche_get_hessian_J: aufpassen, falls nciht alle params_to_optimimze =1');
%     for i=1:length(model.mu_names)
%         if model.optimization.params_to_optimize(i) == 1
%             h_i=(model.mu_ranges{i}(2)-model.mu_ranges{i}(1))/h_range;
%             for j=i:length(model.mu_names)
%                 h_j=(model.mu_ranges{j}(2)-model.mu_ranges{j}(1))/h_range;
%                 if i==j %use other formula for d_ii than for d_ij
%                     %get the nedded values for the
%                     %hessian-finite-difference-formula
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) + 2*h_i;
%                     func_mu_p2hi = model.optimization.objective_function(model,reduced_data); %objective_function in mu + 2*h_i
%                     
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) + h_i;
%                     func_mu_phi = model.optimization.objective_function(model,reduced_data); %objective_function in mu + h_i
%                     
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) - h_i;
%                     func_mu_mhi = model.optimization.objective_function(model,reduced_data); %objective_function in mu - h_i
%                     
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) - 2*h_i;
%                     func_mu_m2hi = model.optimization.objective_function(model,reduced_data); %objective_function in mu - 2*h_i
%                     
%                     model=model_original;
%                     partial_der = (- func_mu_p2hi + 16*func_mu_phi - 30*func_mu + 16*func_mu_mhi - func_mu_m2hi) / (12*h_i^2);
%                     
%                     hessian(i,i) = partial_der;
%                     
%                     
%                 else % i neq j --> calculate d_ij
%                     %get the nedded values for the
%                     %hessian-finite-difference-formula
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) + h_i;
%                     model.(model.mu_names{j}) = model.(model.mu_names{j}) + h_j;
%                     func_mu_phi_phj = model.optimization.objective_function(model,reduced_data); %objective_function in mu + h_i + h_j
%                     
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) + h_i;
%                     model.(model.mu_names{j}) = model.(model.mu_names{j}) - h_j;
%                     func_mu_phi_mhj = model.optimization.objective_function(model,reduced_data); %objective_function in mu + h_i - h_j
%                     
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) - h_i;
%                     model.(model.mu_names{j}) = model.(model.mu_names{j}) + h_j;
%                     func_mu_mhi_phj = model.optimization.objective_function(model,reduced_data); %objective_function in mu - h_i + h_j
%                     
%                     model=model_original;
%                     model.(model.mu_names{i}) = model.(model.mu_names{i}) - h_i;
%                     model.(model.mu_names{j}) = model.(model.mu_names{j}) - h_j;
%                     func_mu_mhi_mhj = model.optimization.objective_function(model,reduced_data); %objective_function in mu - h_i - h_j
%                     
%                     model=model_original; %setting the mus in the model back to their original values
%                     partial_der = (func_mu_phi_phj - func_mu_phi_mhj - func_mu_mhi_phj + func_mu_mhi_mhj)/(4*h_i*h_j); %calculating the partial derivative
%                     
%                     hessian(i,j) = partial_der;
%                     hessian(j,i) = partial_der; %symmetric Hessian --> H(i,j) = H(j,i)
%                 end
%             end %for j=i:length(model.mu_names)
%         end %if model.optimization.params_to_optimize(i) == 1
%     end %for i=1:length(model.mu_names)
%     
%     
% else
%     warning('neither model_data nor reduced_data was given');
%     J=0;
% end


end