function distance2curve
% This practical example shows how to calculate the distance
% from the point b to some parametric curve in M-dimensional 
% space.
% The curve is described in the paper 
% Deb, K. et al., Scalable Test Problems for Evolutionary Multi-Objective Optimization
%
%
% http://strijov.com

% Initial constants 
N = 10000;    % number of points on the true curve (true POF)
M = 3;        % 3D objective space
k = 2;        % curve internal value

% Generate a random vector b
b = 2*rand(M,1);
%-------------------------------------
% a is a vector of the true pof curve 
% a1 is the vector in 1,...,M-1 subspace
a = get_DTLZ6_vector(M);
a1 = a(1:end-1)';
% b is a point of the given pof
b1 = b(1:end-1);
% p is projection b to the line by vector a
p1 = (a1'*b1)/(a1'*a1)*a1;
p = [p1; b(end)];   
% n is norm of p, lies on the true pof curve
n = p/sqrt(sumsqr(p));
% d is the distance from b to curve (to vector a)
d = sqrt(sumsqr(b - n));

% generate the curve as a set of points (partially random in the curve design space)
ds=[];
y = get_DTLZ6_true_pof(N, M, k);
for i = 1:N,
    ds = [ds sumsqr(y(i,:)'-b)];
end;
[ds, idx] = min(ds);
ds = sqrt(ds);
dsa = y(idx,:);
% the random generated nearest point to the vector b

% illustrate the example with the figre
if M == 3, 
    figure; hold on;
    view(72,14);
    axis tight;
    grid on;
    eps = 0.05;
    plot3(y(:,1),y(:,2),y(:,3),'b.'); 
    plot3(b(1),b(2),b(3),'b*'); 
    text(b(1)+eps,b(2)+eps,b(3),'b'); ;
    plot3(p(1),p(2),p(3),'r*'); 
    text(p(1)+eps,p(2)+eps,p(3),'p'); ;
    plot3(n(1),n(2),n(3),'g+'); 
    text(n(1)+eps,n(2)+eps,n(3),'n'); 
    
    plot3([0 b(1)],[0 b(2)],[0 b(3)],'b-'); 
    plot3([0 p(1)],[0 p(2)],[0 p(3)],'r-'); 
    plot3([0 n(1)],[0 n(2)],[0 n(3)],'r-');   
    
    plot3([b(1) p(1)],[b(2) p(2)],[b(3) p(3)],'r-');   
    plot3([b(1) n(1)],[b(2) n(2)],[b(3) n(3)],'g-');   
    
    plot3(dsa(1),dsa(2),dsa(3),'ro');    
end;

fprintf(1,'\n the distance to the actual nearest point on the curve is %f',d);
fprintf(1,'\n the distance to the random nearest point on the curve is %f',ds);
return
%===========================================================================
function y = get_DTLZ6_true_pof(N, M, k)
mode=1; %1 for true pof
x=[];
indim = k + M - 1;
for i=1:indim,
    if i < M,
        x=[x rand(N,1)];
    else,
        x=[x 0.0*ones(N,1)];
    end;
end;  


PP.M = M;
PP.k = k;
y = DTLZ6_obj(x,PP);

return

%=================================================================
function y = DTLZ6_obj(x,PP)
% MO test problem DTLZ2
% Source: Scalable Test Problems for Evolutionary Multi-Objective Optimization by K. Deb et al. July 17, 2001
% Rem: r = s the constant parameter of the algorithm 

indim = PP.k + PP.M - 1;
outdim = PP.M;

if (PP.k < 1) | (PP.M < 1), 
    error 'Invalid problem parameters';
end;
if indim ~= size(x,2), 
    error 'The input vector space dimension doesnt correspond to PP.indim. Consider PP.indim = size(x,2)';
end;

% definifion of the function value gr=g(r);
r=10; a=100;
%gr=r^10;            % see page 6, 9 top
%gr=x(1,1);         % see page 6, 5 bt.
%gr=x in [-1,1];    % see page 7, 6 bt.
%gr=0;            % see page 8, 9 top
% for the curve problem:
gr=1/r^a; % a>>1    % see page 8, 11 bt.

% The mapping procedure, t stands for theta
t=(pi/(4*(1+gr))*(1+2*gr*x));
t(:,1)=x(:,1);
t(:,PP.M:end)=x(:,PP.M:end);

for r=1:size(x,1),
%g = sum((x(r,PP.M:end) - 0.5).^2);      % is for DTLZ5 
g = sum(x(r,PP.M:end).^0.1);           % is for DTLZ6 (the only difference)

for i=1:PP.M,
    cos_nom = 1;
    for j=1:PP.M-i, 
        cos_nom=cos_nom*cos(t(r,j)*pi/2);   % change x for t
    end;
    
    if i==1,
        sin_nom=1;
    else,
        sin_nom=sin(t(r,PP.M-i+1)*pi/2);    % change x for t
    end;
    y(r,i) =(1 + g)*cos_nom*sin_nom;
  end;
end;    

%=================================================================

function y = get_DTLZ6_vector(M)
  for i=1:M,
    cos_nom = 1;
    for j=1:M-i, 
        cos_nom=cos_nom*cos(pi^2/8); 
    end;

    if i==1,
        sin_nom=1;
    else,
        sin_nom=sin(pi^2/8);
    end;
    y(i) = cos_nom*sin_nom;
  end;