%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%      Matlab scrip to plot the extendd oracle penalty function shape
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_oracle_penalty

clc

acc = 0.1;
oracle = 0;

x = zeros(1,1);
y = zeros(1,1);
z = zeros(1,1);

for s = 1:200
    for t = 1:100
        
        f   = -10 + s/10;
        res = t/10;
        
        x(s,t) = f;
        y(s,t) = res;
        
        z(s,t) = oracle_penalty(f,res,oracle,acc);
        
    end
end

mesh(x(1:200,1:100),y(1:200,1:100),z(1:200,1:100))
xlabel('Objective Function')
ylabel('Residual')
zlabel('Penalty Function')
saveas(gcf,'plot.eps', 'psc2')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     IMPLEMENTATION OF THE EXTENDED ORACLE PENALTY FUNCTION
%     ------------------------------------------------------
%
%     Input arguments:
%     ----------------
%     f           Objective function value
%     res         Residual function value
%     omega       Oracle parameter
%     acc         Accuracy for constraint violation 
%
%     Output arguments:
%     -----------------
%     p           Penalty function value
%     
%     ------------------------------------------------------
%
%     Reference:
%
%     Schlueter M., Gerdts M.: "The oracle penalty method",
%     Journal of Global Optimization (Springer), 
%     Vol. 47, Issue 2, Page 293-325, 2010.
%
%     ------------------------------------------------------
%
%     Author (C) :   Dr. Martin Schlueter
%                    Information Initiative Center (IIC),
%                    Division of Large Scale Computing Systems,
%                    Hokkaido University, JAPAN.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ p ] = oracle_penalty( f, res, omega, acc)    

if (f <= omega) && (res <= acc)

    p = f - omega;
    
else

    if f <= omega
        
        p = res;
        
    else
        
        if res < (f-omega)/3.0
            
            alpha = (abs(f-omega) * 0.80754991-res) / ((f-omega)-res);
            
        end
        
        if  (res <= (f-omega)) && (res >= (f-omega)/3.0)
            
            alpha = 1.0 - 0.5 / sqrt((f-omega)/res);
            
        end
        
        if  res > (f-omega)
            
            alpha = 0.5 * sqrt((f-omega)/res);
            
        end
        
        p = alpha * (f-omega) + (1.0-alpha) * res;

    end
end




