function MO_Animate(varargin)
% This function generates objective space images showing why
% sum-weighted optimizer can not find all non-dominated
% solutions for non convex objective spaces in multi-ojective
% optimization
%
% Guillaume JACQUENOT


if nargin == 0
    Simu = 'Convex';
    % Simu = 'NonConvex';
    save_pictures = true;
    interpreter = 'none';
end


switch Simu
    case 'NonConvex'
        a = 0.1;
        b = 3;
        stepX = 1/200;
        stepY = 1/200;
    case 'Convex'
        a = 0.2;
        b = 1;
        stepX = 1/200;
        stepY = 1/200;
end

[X,Y] = meshgrid( 0:stepX:1,-2:stepY:2);

F1 = X;
F2 = 1+Y.^2-X-a*sin(b*pi*X);

figure;
grid on;
hold on;
box on;
axis square;
set(gca,'xtick',0:0.2:1);
set(gca,'ytick',0:0.2:1);

Ttr = get(gca,'XTickLabel');
Ttr(1,:)='0.0';
Ttr(end,:)='1.0';
set(gca,'XTickLabel',[repmat(' ',size(Ttr,1),1) Ttr]);

Ttr = get(gca,'YTickLabel');
Ttr(1,:)='0.0';
Ttr(end,:)='1.0';
set(gca,'YTickLabel',[repmat(' ',size(Ttr,1),1) Ttr]);

if strcmp(interpreter,'none')
    xlabel('f1','Interpreter','none');
    ylabel('f2','Interpreter','none','rotation',0);
else
    xlabel('f_1','Interpreter','Tex');
    ylabel('f_2','Interpreter','Tex','rotation',0);
end

set(gcf,'Units','centimeters')
set(gcf,'OuterPosition',[3 3 3+6 3+6])
set(gcf,'PaperPositionMode','auto')

[minF2,minF2_index] = min(F2);
minF2_index = minF2_index + (0:numel(minF2_index)-1)*size(X,1);

O1 = F1(minF2_index)';
O2 = minF2';

[pF,Pareto]=prtp([O1,O2]);

fill([O1( Pareto);1],[O2( Pareto);1],repmat(0.95,1,3));

text(0.45,0.75,'Objective space');
text(0.1,0.9,'\leftarrow Optimal Pareto front','Interpreter','TeX');

plot(O1( Pareto),O2( Pareto),'k-','LineWidth',2);
plot(O1(~Pareto),O2(~Pareto),'.','color',[1 1 1]*0.8);
V1 = O1( Pareto); V1 = V1(end:-1:1);
V2 = O2( Pareto); V2 = V2(end:-1:1);

O1P = O1( Pareto);
O2P = O2( Pareto);

O1PC = [O1P;max(O1P)];
O2PC = [O2P;max(O2P)];
ConvH = convhull(O1PC,O2PC);
ConvH(ConvH==numel(O2PC))=[];
c = setdiff(1:numel(O1P), ConvH);

% Non convex
O1PNC = O1PC(c);

[temp, I1] = min(O1PNC);
[temp, I2] = max(O1PNC);

if ~isempty(I1) && ~isempty(I2)
    plot(O1PC(c),O2PC(c),'-','color',[1 1 1]*0.7,'LineWidth',2);
end

p1 = (V2(1)-V2(2))/(V1(1)-V1(2));
hp = plot([0 1],[p1*(-V1(1))+V2(1) p1*(1-V1(1))+V2(1)]);
delete(hp);


Histo_X = [];
Histo_Y = [];
coeff = 0.02;
Sq1 = coeff *[0 1 1 0 0;0 0 1 1 0];
compt = 1;
for i = 2:1:length(V1)-1
    if ismember(i,ConvH)
        p1 = (V2(i+1)-V2(i-1))/(V1(i+1)-V1(i-1));
        x_inter = 1/(1+p1^2)*(p1^2*V1(i)-p1*V2(i));
        hp1 = plot([0 1],[p1*(-V1(i))+V2(i) p1*(1-V1(i))+V2(i)],'k');
        % hp2 = plot([x_inter],[-x_inter/p1],'k','Marker','.','MarkerSize',8)
        hp3 = plot([0 x_inter],[0 -x_inter/p1],'k-');
        hp4 = plot([x_inter 1],[-x_inter/p1 -1/p1],'k--');
        hp5 = plot(V1(i),V2(i),'ko','MarkerSize',10);

        % Plot the square for perpendicular lines
        alpha = atan(-1/p1);
        Mrot = [cos(alpha) -sin(alpha);sin(alpha) cos(alpha)];
        Sq_plot = repmat([x_inter;-x_inter/p1],1,5) + Mrot * Sq1;
        hp7 = plot(Sq_plot(1,:),Sq_plot(2,:),'k-');

        Histo_X = [Histo_X V1(i)];
        Histo_Y = [Histo_Y V2(i)];
        hp6 = plot(Histo_X,Histo_Y,'k.','MarkerSize',10);

        w1 = p1/(p1-1);
        w2 = 1-w1;
        Fweight_sum = V1(i)*w1+w2*V2(i);
        Fweight_sum = floor(1e3*Fweight_sum )/1e3;

        w1 = floor(1000*w1)/1e3;
        str1 = sprintf('%.3f',w1);
        str2 = sprintf('%.3f',1-w1);
        str3 = sprintf('%.3f',Fweight_sum);
        if (strcmp(str1,'0.500')||strcmp(str1,'0,500')) && strcmp(Simu,'NonConvex')
            disp('Two solutions');
        end
        title(['\omega_1 = ' str1 '  &  \omega_2 = ' str2 '  &  F = ' str3],'Interpreter','TeX');
        axis([0 1 0 1]);
        file = ['Frame' num2str(1000+compt)];
        if save_pictures
            saveas(gcf, file, 'epsc');
        end
        compt = compt +1;
        pause(0.001);
        delete(hp1);
        delete(hp3);
        delete(hp4);
        delete(hp5);
        delete(hp6);
        delete(hp7);
    end
end
disp(['Number of frames :' num2str(length(V1))]);
return;

function [A varargout]=prtp(B)
% Let Fi(X), i=1...n, are objective functions
% for minimization.
% A point X* is said to be Pareto optimal one
% if there is no X such that Fi(X)<=Fi(X*) for
% all i=1...n, with at least one strict inequality.
% A=prtp(B),
% B - m x n input matrix: B=
% [F1(X1) F2(X1) ... Fn(X1);
%  F1(X2) F2(X2) ... Fn(X2);
%  .......................
%  F1(Xm) F2(Xm) ... Fn(Xm)]
% A - an output matrix with rows which are Pareto
% points (rows) of input matrix B.
% [A,b]=prtp(B). b is a vector which contains serial
% numbers of matrix B Pareto points (rows).
% Example.
% B=[0 1 2; 1 2 3; 3 2 1; 4 0 2; 2 2 1;...
%    1 1 2; 2 1 1; 0 2 2];
% [A b]=prtp(B)
% A =
%      0     1     2
%      4     0     2
%      2     2     1
% b =
%      1     4     7
A=[]; varargout{1}=[];
sz1=size(B,1);
jj=0; kk(sz1)=0;
c(sz1,size(B,2))=0;
bb=c;
for k=1:sz1
    j=0;
    ak=B(k,:);
    for i=1:sz1
        if i~=k
            j=j+1;
            bb(j,:)=ak-B(i,:);
        end
    end
    if any(bb(1:j,:)'<0)
        jj=jj+1;
        c(jj,:)=ak;
        kk(jj)=k;
    end
end
if jj
  A=c(1:jj,:);
  varargout{1}=kk(1:jj);
else
  warning([mfilename ':w0'],...
    'There are no Pareto points. The result is an empty matrix.')
end
return;