gradient_opt_non_iter_err.m

function [opt_data, model]=gradient_opt_non_iter_err(model,varargin)
%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_to_optimize(model);
%get the according boundaries 
[lower_bound, upper_bound] = get_bound_to_optimize(model);


if strcmp(model.optimization.opt_mode,'detailed')
    %if strcmp(inputname(2),'model_data')
        model_data  = varargin{1};
        detailed_data = [];
    %else
        %error('Its not called model data')
    %end
       % model.optimization.opt_mode = 'detailed';
elseif strcmp(model.optimization.opt_mode,'reduced')
    %if strcmp(inputname(2),'reduced_data')
    %    error('input argument is not called reduced data')
    %end
    if nargin >2
        reduced_data = varargin{1};
        model_data = varargin{2};
        detailed_data = varargin{3};
            model.optimization.opt_mode = 'reduced';
    else
        reduced_data = varargin{1};
            model.optimization.opt_mode = 'reduced';
    end
end


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

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


      
grad_f_x    = model.optimization.get_Jacobian(model,model_data,[],[])*min_or_max; nr_grad = nr_grad+1;
norm_grad_f_x = norm(grad_f_x);

% START OF OPTIMIZATION
while (norm(grad_f_x) > tol) && (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] = stepsize_armijo(model, model_data, output, x, d);
            nr_fct = nr_fct+nr_fct_stepsize;
           
        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;
            
            
        %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: ',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_to_optimize(model,x); % write the actual parameter setting in the model
    
    %additional break condition:
    %if there were no major changes in the last n steps --> finish!
    n=5;
    if (size(PM,1) > n) && (length(output) > n) %already more than n values to compare?
        PM_diff=[];
        output_diff=[];
        for k=1:n
            PM_diff = [PM_diff; abs(PM(end-k,:)-PM(end-k+1,:))]; %difference of the mus in the last n steps
            output_diff = [output_diff; abs(output(end-k)-output(end-k+1))]; %difference of the objective function in the last n steps
        end
        if (max(PM_diff) <= tol*1e-2) & (max(output_diff) <= tol*1e-2)
            disp('no major change in mu and objective function')
            disp('leaving gradient_opt')
            break
        end
    end
    
    [grad_f_x,dummy,output]=model.optimization.get_Jacobian(model,model_data,[],[]);
    norm_grad_f_x = [norm_grad_f_x,norm(grad_f_x)];
    if model.verbose >3
        disp(['actual output: ',num2str(output)]);
    end
    grad_f_x = grad_f_x * min_or_max;%now calculate the gradient with the new parameter setting
    nr_grad = nr_grad+1; 
    disp(['Actual gradient norm :',num2str(norm_grad_f_x)]);
end
% END OF OPTIMIZATION (detailed?)


%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.nr_opti_iter=opti_iter;
opt_data.norm_grad_f_x = norm_grad_f_x;
%MAL RAUSGENOMMEN!!!
% setting the mus values in the model back to the default values
%for k=1:length(model.mu_names)
%    model.(model.mu_names{k})=mu_original(k);
%end

else  %here: reduced optimization

    
[grad_f_x,Delta_grad,output,data] = 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);
nRB = data.nRB;
nRB_der = data.nRB_der;
[g_H,model] = model.gamma_H(model);
[L_H,model] = model.L_H(model);
break_value_vec = 2*(g_H)^2*(L_H)*(norm_grad_f_x);
Delta_mu_vec = 2*g_H * norm(norm_grad_f_x);

used_detailed_calculations = 0;
% PM          = x;        %for opt_data: array with all the mu_i
% output      = [];       %array containing the f(mu_k)
% max_iter    = 40;   %maximal number of iterations of the gradient method
% opti_iter   = 0;    %counter
% if isfield(model,'lipschitz_type')&&(strcmp(model.lipschitz_type,'Hessian'))
%     Delta_mu = zeros(length(grad_f_x),1);
% else
%     Delta_mu = 0;       %Error bound for reduced optimization optimal parameters
% end

% %Error estimation:
%     %fixing constants:
%     switch model.optimization.stepsize_rule
%         case {'qoutient'}
%             epsilon_m=zeros(max_iter+1,1);
%             C_alpha=1;%upper bound for the optimization stepsize
%             [C_L,model] = model.get_lipschitz_constant(model);
%             
%         otherwise
%             disp('unknown stepsize rule');
%             opt_data = [];
%             return;
%     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_to_optimize(model,x); % write the actual parameter setting in the model
    
    %additional break condition:
    %if there were no major changes in the last n steps --> finish!
%     n=5;
%     if (size(PM,1) > n) && (length(output) > n) %already more than n values to compare?
%         PM_diff=[];
%         output_diff=[];
%         for k=1:n
%             PM_diff = [PM_diff; abs(PM(end-k,:)-PM(end-k+1,:))]; %difference of the mus in the last n steps
%             output_diff = [output_diff; abs(output(end-k)-output(end-k+1))]; %difference of the objective function in the last n steps
%         end
%         if (max(max(PM_diff)) <= tol*1e-2) && (max(output_diff) <= tol*1e-2)
%             disp('no major change in mu and objective function')
%             disp('leaving gradient_opt')
%             break
%         end
%     end
%     
%     %Error estimation:
%     n_Delta = size(Delta_mu,2);
%     disp('next estimates')
%     disp(['augmenting constants:',num2str(1+C_alpha*norm(C_L))]);
%     disp(['norm of the gradient:',num2str(norm(grad_f_x))]);
%     disp(['Error estimator Gradient:',num2str(Delta_grad(:,end)')]);
%     disp(['actual adding:',num2str(C_alpha*Delta_grad(:,end)')]);
%     
%     if isfield(model,'lipschitz_type')&&(strcmp(model.lipschitz_type,'Hessian'))
%         Delta_mu_next = Delta_mu(:,n_Delta)+C_alpha.*(C_L*Delta_mu(:,n_Delta))+epsilon_m(n_Delta)*abs(grad_f_x(:,:))'+C_alpha*Delta_grad(:,end);
%     else
%         Delta_mu_next = Delta_mu(:,n_Delta)+C_alpha.*(C_L*Delta_mu(:,n_Delta))+epsilon_m(n_Delta)*norm(grad_f_x(:,:))'+C_alpha*norm(Delta_grad(:,end));
%     end
%     Delta_mu = [Delta_mu,Delta_mu_next];
%     
%     disp(['actual error estimate for parameters:',num2str(Delta_mu_next')]);
%     
    %Calculate gradient at the new parameter point:
    [grad_f_x, Delta_grad,output, data] = model.optimization.get_Jacobian(model,reduced_data);%now calculate the gradient with the new parameter setting
    norm_grad_f_x = norm(grad_f_x);
    %output=output_array{1};
    grad_f_x=grad_f_x*min_or_max;  
    nr_grad = nr_grad+1; 
    nRB = [nRB,data.nRB];
    nRB_der = [nRB_der, data.nRB_der];
    
    %decide to continue or not
    %break_value = 2*model.gamma_H^2*model.L_H*(norm_grad_f_x+norm(Delta_grad(:,end)));
    [g_H,model] = model.gamma_H(model);
    [L_H,model] = model.L_H(model);
    break_value = 2*(g_H)^2*(L_H)*(norm_grad_f_x+norm(Delta_grad(:,end)));
    break_value_vec = [break_value_vec, break_value];
    Delta_mu_vec = [Delta_mu_vec, 2*g_H * norm(norm_grad_f_x+norm(Delta_grad(:,end)))];
    disp(['Actual gradient norm :',num2str(norm_grad_f_x)]);
    disp(['Value of breaking criteria: ', num2str(break_value)]);
    %if (break_value<=1)&&(norm_grad_f_x<tol)
    if (norm_grad_f_x<tol)
        continue_optimization = 0;
    end
   
end
% END OF OPTIMIZATION

%noniterative error estimator
%Delta_mu = 2*model.gamma_H * norm(norm_grad_f_x+norm(Delta_grad(:,end)));
if break_value>1
    %...the requirement is not fullfilled
    if isfield(model.optimization,'allow_constant_calculation')&&model.optimization.allow_constant_calculation
        gamma = calculate_gamma_const(model.base_model);
        Lip_H = calculate_Lipschitz_const(model.base_model);
        break_value = 2*(gamma)^2*(Lip_H)*(norm_grad_f_x+norm(Delta_grad(:,end)));
        used_detailed_calculations = 1;
        g_H = gamma;
        L_H = Lip_H;
        if break_value<=1
            Delta_mu = 2*gamma * norm(norm_grad_f_x+norm(Delta_grad(:,end)));
        else
            disp('Error estimator not valid!!!!!!!!!!!!')
            Delta_mu = -1;
        end
    else
        %Error estimator not valid
        disp('Error estimator not valid!!!!!!!')
        Delta_mu = -1;
    end
else
    Delta_mu = 2*model.gamma_H(model) * norm(norm_grad_f_x+norm(Delta_grad(:,end)));
end
%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.nr_opti_iter=opti_iter;
opt_data.Delta_mu = Delta_mu;
opt_data.break_value = break_value;
opt_data.norm_grad_f_x = norm_grad_f_x;
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