porsche_get_Jacobian.m

function [jacobian,Delta_J] = porsche_get_Jacobian(model, data)
%function [Jacobian,Delta_J] = porsche_get_Jacobian(model, data)
%
%Function returning the Jacobian of the output functional
%
%For detailed optimization the Jacobian is calculated numerically via
%Finite Differences
%For reduced optimization the function checks if the @f$ \partial_{\mu_i}
%\Theta @f$ are available in the model. If not it is also calculated
%numerically via Finite Differences
%
%data can be either "model_data" in the detailed case or "reduced_data"
%in the reduced case.
%
% In the reduced case an error bound is given by Delta_J
%
% Markus Dihlmann 30.04.2010
% Oliver Zeeb
%

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

Delta_J = 0;
h_range = 100000; % for constructing the difference quotient: how small is the variaton of the parameter 

%new: with flag model.optimization.opt_mode


if (~model.optimization.derivatives_available||~model.compute_derivative_info)
    %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
    jacobian=[];
    for k=1:length(model.mu_names)
        h = (model.mu_ranges{k}(2)-model.mu_ranges{k}(1))/h_range;  % ranges must be intervals
        if model.optimization.params_to_optimize(k) == 1
            if model.mu_ranges{k}(2)-model.(model.mu_names{k}) >= h  % is it possible to get the righthanded difference quotient?
                old_mu = model.(model.mu_names{k});
                model.(model.mu_names{k}) = model.(model.mu_names{k}) + h; %getfield(model,model.mu_names{k})
                func = model.optimization.objective_function(model, data);           %value of the functional after changing one parameter
                grad_comp = (func - func_mu)/h;   %righthanded difference quotient
                model.(model.mu_names{k}) = old_mu; %set model to initial values for the next cycle run
            else
                old_mu = model.(model.mu_names{k});
                model.(model.mu_names{k}) = model.(model.mu_names{k}) - h;
                func = model.optimization.objective_function(model, data);
                grad_comp = (func_mu - func)/h;   % lefthanded difference quotient
                model.(model.mu_names{k}) = old_mu; %set model to initial values for the next cycle run
            end
            jacobian=[jacobian,grad_comp]; % add component to the end of gradient
        end
    end
    
    if(model.verbose>=8)
        disp(['calculated gradient is: ' num2str(J)])
    end
% else
%     %calculate derivative via derivative pde simulation
%     sim_data_der = detailed_simulation(model, model_data);
%     y_end=size(sim_data_der.y_der,1);
%     J=sim_data_der.y_der(y_end,:)';
end



%old: without flag model.optimization.opt_mode
% if strcmp(inputname(2),'model_data')
%         model_data = varargin{1}; 
%                
%         if (~model.optimization.derivatives_available||~model.compute_derivative_info)
%             %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
%             jacobian=[];
%             for k=1:length(model.mu_names)
%                 h = (model.mu_ranges{k}(2)-model.mu_ranges{k}(1))/h_range;  % ranges must be intervals
%                 if model.optimization.params_to_optimize(k) == 1
%                     if model.mu_ranges{k}(2)-model.(model.mu_names{k}) >= h  % is it possible to get the righthanded difference quotient?
%                         old_mu = model.(model.mu_names{k});
%                         model.(model.mu_names{k}) = model.(model.mu_names{k}) + h; %getfield(model,model.mu_names{k})
%                         func = model.optimization.objective_function(model, model_data);           %value of the functional after changing one parameter
%                         grad_comp = (func - func_mu)/h;   %righthanded difference quotient
%                         model.(model.mu_names{k}) = old_mu; %set model to initial values for the next cycle run
%                     else
%                         old_mu = model.(model.mu_names{k});
%                         model.(model.mu_names{k}) = model.(model.mu_names{k}) - h;
%                         func = model.optimization.objective_function(model, model_data);
%                         grad_comp = (func_mu - func)/h;   % lefthanded difference quotient
%                         model.(model.mu_names{k}) = old_mu; %set model to initial values for the next cycle run
%                     end
%                     jacobian=[jacobian,grad_comp]; % add component to the end of gradient
%                 end
%             end
%            %-------------------------------------------
%            %End of detailed calculation of the gradient
%            %-------------------------------------------        
%             if(model.verbose>=8)
%                 disp(['calculated gradient is: ' num2str(J)])
%             end
% %         else
% %             %calculate derivative via derivative pde simulation
% %             sim_data_der = detailed_simulation(model, model_data);
% %             y_end=size(sim_data_der.y_der,1);
% %             J=sim_data_der.y_der(y_end,:)';
%         end
%         
%         
% elseif strcmp(inputname(2), 'reduced_data')
%         reduced_data= varargin{1};
%                 
%         %are derivatives available?
%         %yes -> calculation of gradient by evolution equation
%         %no -> finite difference calculation of gradient using reduced
%         %simulations
%         
%         if (model.optimization.derivatives_available)
%             jacobian=[];
%             sim_data = rb_derivative_simulation(model, reduced_data);
%             %sim_data = model.rb_reconstruction(model, model_data, detailed_data, sim_data);
%             for i=1:size(sim_data.s_der,2)
%                 y=sim_data.s_der{i}(end);
%                 jacobian=[jacobian,y];
%             end
%             
%             Delta_J = sim_data.Delta_s_grad;
%         else
%            
%             func_mu = model.optimization.objective_function(model, reduced_data); %value of the functional at the actual parameter set, needed in every cycle run 
%             jacobian=[];
%             for k=1:length(model.mu_names)
%                 h = (model.mu_ranges{k}(2)-model.mu_ranges{k}(1))/h_range;  % ranges must be intervals
%                 if model.optimization.params_to_optimize(k) == 1
%                     if model.mu_ranges{k}(2)-model.(model.mu_names{k}) >= h     % is it possible to get the righthanded difference quotient?
%                         old_mu = model.(model.mu_names{k});
%                         model.(model.mu_names{k}) = model.(model.mu_names{k}) + h; %getfield(model,model.mu_names{k})
%                         func = model.optimization.objective_function(model, reduced_data);           %value of the functional after changing one parameter
%                         grad_comp = (func - func_mu)/h;   %righthanded difference quotient
%                         model.(model.mu_names{k}) = old_mu; %set model to initial values for the next cycle run
%                     else
%                         old_mu = model.(model.mu_names{k});
%                         model.(model.mu_names{k}) = model.(model.mu_names{k}) - h;
%                         func = model.optimization.objective_function(model, reduced_data);
%                         grad_comp = (func_mu - func)/h;  % lefthanded difference quotient
%                         model.(model.mu_names{k}) = old_mu; %set model to initial values for the next cycle run
%                     end
%                     jacobian=[jacobian,grad_comp]; % add component to the end of gradient
%                 end
%             end
% 
%          end
%         
%         
%         
%               
%         
% else
%     disp('neither model_data nor reduced_data was given');
%     J=0;
% end