gradient_opt.m

function [opt_data, model]=gradient_opt(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
    

nr_fct = 0;     %counter for number of function calls
nr_grad = 0;    %counter for number of gradient calculations

%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    = 50;   %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)];
    %output=output_array{1};
    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
% 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] = 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);
% 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 {'fixed_stepsize'}
            C_alpha = model.optimization.stepsize_fix;
            [C_L,model] = model.get_lipschitz_constant(model);
            epsilon_m=zeros(max_iter+1,1);
        case {'quotient'}
            epsilon_m=zeros(max_iter+1,1);
            C_alpha=1;%upper bound for the optimization stepsize
            [C_L,model] = model.get_lipschitz_constant(model);
            
        case {'armijo'}
            epsilon_m = ones(max_iter+1,1).*model.optimization.armijo_t_min*model.optimization.stepsize_factor;
            C_alpha = 1;
            [C_L,model] = model.get_lipschitz_constant(model);
        otherwise
            disp('unknown stepsize rule');
            opt_data = [];
            return;
    end
    

% 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, reduced_data, output, x, d);
            nr_fct = nr_fct+nr_fct_stepsize;
           [C_L,model] = model.get_lipschitz_constant(model); %update Lipschitz constant
        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] = 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_f_x)];
    %output=output_array{1};
    grad_f_x=grad_f_x*min_or_max;  
    nr_grad = nr_grad+1; 
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.nr_opti_iter=opti_iter;
opt_data.Delta_mu = Delta_mu;
opt_data.norm_grad_f_x = norm_grad_f_x;

end