% Software for the computation of a data matrix from given data points
%
% Copyright (c) 2024 Tomass Andersons, Martina Beese, University of Rostock, Department of Mathematics.
% 
% The software can be used for academic, research and other similar noncommercial uses. The user
% acknowledges that this software is still in the development stage and that it is provided by the
% copyright holders and contributors "as is" and any express or implied warranties, fitness for a
% particular purpose are disclaimed. In no event shall the copyright owner or contributors be liable
% for any direct, indirect, incidental, special, exemplary, or consequential damages.
% 
% The copyright holders provide no reassurances that the source code provided does not infringe any
% patent, copyright, or any other intellectual property rights of third parties. The copyright
% holders disclaim any liability to any recipient for claims brought against recipient by any third
% party for infringement of that parties intellectual property rights.

%% Info
% The GUI can be started by running this script.

%%
function recover_matrix_PAA

fig = figure(222);
pause(.1);

scrsz = get(0,'ScreenSize');

set(fig,'Position',[1/6*scrsz(3:4) 4/5*scrsz(3) 1/2*scrsz(4)]);

pause(.5);
figpos = get(fig,'Position');

set(fig,'Name', 'Polytope adaptation algorithm','NumberTitle','off','MenuBar','none','Toolbar','none','color','white');

if ispc
    basicfsz = 5;
    basichi = 0.05;
else
    basicfsz = 6.5;
    basichi = 0.05;
end
Styles = {'bitstream charter';'clean';'Courier';'courier 10 pitch';'frank ruehl';'Helvetica';'Helvetica-Narrow';'itc bookman';'urw bookman l'};
sty = Styles{6};
bgc = 0.95*ones(3,1);
% lcs = [0 0 1;0 0.5 0;1 0 0;0 0.5 0.5;0.5 0 0.5;0.5 0.5 0];
lcs = [0 0.4470 0.7410; 0.8500 0.3250 0.0980; 0.9290 0.6940 0.1250; 0.4940 0.1840 0.5560; 0.4660 0.6740 0.1880; 0.3010 0.7450 0.9330];
LCS = [0.85 0.9 1; 0.995 0.8250 0.75; 1 0.9 0.75; 0.9 0.6 1; 0.7 1 0.5; 0.7 0.9 1];

% title names
titles(1) = uicontrol('Style','text','String','target I_S','HorizontalAlignment','center','Fontsize',3*basicfsz,'BackgroundColor','w',...
    'Units','normalized','FontName',sty,'Position',[0.02, 0.89, 0.3, basichi]);
titles(2) = uicontrol('Style','text','String','target I_C','HorizontalAlignment','center','Fontsize',3*basicfsz,'BackgroundColor','w',...
    'Units','normalized','FontName',sty,'Position',[0.35, 0.89, 0.3, basichi]);
titles(3) = uicontrol('Style','text','String','D','HorizontalAlignment','center','Fontsize',3*basicfsz,'BackgroundColor','w',...
    'Units','normalized','FontName',sty,'Position',[0.68, 0.89, 0.3, basichi]);


% subplots
AxesIPS = axes('Fontsize',2*basicfsz,'Units','normalized','FontName',sty,'Position',[0.02, 0.2, 0.3, 0.65]);
plot(AxesIPS,1,0,'w');

AxesIPC = axes('Fontsize',2*basicfsz,'Units','normalized','FontName',sty,'Position',[0.35, 0.2, 0.3, 0.65]);
plot(AxesIPC,1,0,'w');

AxesD = axes('Fontsize',2*basicfsz,'Units','normalized','FontName',sty,'Position',[0.68, 0.2, 0.3, 0.65]);
plot(AxesD,1,0,'w');

% Note
note=uicontrol('Style','text','String','Please draw target I_S','HorizontalAlignment','left','Fontsize',2*basicfsz,'BackgroundColor','w',...
    'Units','normalized','FontName',sty,'Position',[0.01, 0.0, 0.5, basichi]);

% The buttons with all UserData in emptybutton
emptybutton = uicontrol('Style','pushbutton','Position',[0.02, 0.8, 0.14, basichi],'Visible','off','Callback',{@emptybutton_CB});


Button(3) = uicontrol('Style','pushbutton','UserData',0,'String','Save Results','HorizontalAlignment','center','Fontsize',2*basicfsz,...
    'Units','normalized','FontName',sty,'Position',[0.78, 0.05, 0.1, basichi],'Enable','off',...
    'Callback', {@buttonsave, emptybutton});

Button(2) = uicontrol('Style','pushbutton','UserData',0,'String','Confirm target I_C','HorizontalAlignment','center','Fontsize',2*basicfsz,...
    'Units','normalized','FontName',sty,'Position',[0.45, 0.05, 0.1, basichi],'Enable','off',...
    'Callback', {@buttonD2, emptybutton,Button,AxesIPS,AxesIPC,AxesD,note});

Button(1) = uicontrol('Style','pushbutton','UserData',0,'String','Confirm target I_S','HorizontalAlignment','center','Fontsize',2*basicfsz,...
    'Units','normalized','FontName',sty,'Position',[0.12, 0.05, 0.1, basichi],'Enable','off',...
    'Callback', {@buttonD1, emptybutton, Button,AxesIPS, AxesIPC,note});

onoffButton(Button,1,[])

% Button(2) = uicontrol('Style','pushbutton','UserData',0,'String','Compute polytopes','HorizontalAlignment','center','Fontsize',2*basicfsz,...
%     'Units','normalized','FontName',sty,'Position',[0.02, 0.4, 0.14, basichi],'Enable','off',...
%     'Callback', {@buttonEDP, emptybutton, Button, options, optionsray, AxesU, AxesV,AxesD1, AxesD2, lcs});


emptybutton_CB(emptybutton)
drawIPS(emptybutton,AxesIPS);

set(gcf,'WindowButtonDownFcn', {@fcnbd, fig, emptybutton, Button, AxesIPS, AxesIPC})

end

function emptybutton_CB(hObject, eventdata)
data.ismatrix = 0;
data.y=[];
data.z=[];
data.drIPS=1;
data.drIPC=0;

N=3;
data.D0 = 3*diag(ones(1,N)) + diag(ones(1,N-1),1) + diag(ones(1,N-1),-1);
set(hObject,'UserData',data);
end

function buttonsave(hObject, eventdata, Empty_Button)
data = get(Empty_Button,'UserData');
if isfield(data,'D') == 0
    return;
end
D=data.D;

% data.PathName = actpwd(1:end-10);
% cd(data.PathName);
% [FileName,PathName,FilterIndex] = uigetfile('*.mat','Save Matrix D');

uisave({'D'},'results.mat') %[pwd,'PAA_',data.FileName]
end

function buttonD1(hObject, eventdata, Empty_Button, Button,AxesIPS,AxesIPC,note)
data = get(Empty_Button,'UserData');

if isempty(data.y)
    return;
end

% check if in FP
if length(convhull(data.FPS))~=length(convhull([data.FPS;data.y]))
    warning('The target I_S is not included in the outer polygon')
    data.y=[];
    drawIPS(Empty_Button,AxesIPS);
    set(Empty_Button,'UserData',data);
    return;
end

% check if 0 enthalten
if length(convhull(data.y))~=length(convhull([data.y; 0 0]))
    % neu malen lassen, Bild mit 0 initialisieren
    warning('The origin is not included in target I_S')
    data.y=[];
    drawIPS( Empty_Button,AxesIPS);
    set(Empty_Button,'UserData',data);
    return;
end

% beide Bilder neu
data.D1=[ones(size(data.y,1),1),data.y]*data.V0';
set(Empty_Button,'UserData',data);

% function draw after D1
data=drawD1(Empty_Button, AxesIPS,AxesIPC);

data.drIPS=0;
data.drIPC=1;

% onoffButton(Button,2,1)
onoffButton(Button,2,[])
% onoffButton(Button,[],1)

set(hObject,'Enable','off')
set(note,'String','Please draw target I_C')

set(Empty_Button,'UserData',data);
end

function buttonD2(hObject, eventdata, Empty_Button,Button,AxesIPS,AxesIPC,AxesD,note)
data = get(Empty_Button,'UserData');


if isempty(data.z)
    return;
end

% check if in FP
if length(convhull(data.FPC))~=length(convhull([data.FPC;data.z]))
    warning('The target I_C is not included in the outer polygon')
    data.z=[];
    data=drawD1(Empty_Button,AxesIPS,AxesIPC);
    set(Empty_Button,'UserData',data);
    return;
end

% check if 0 enthalten
if length(convhull(data.z))~=length(convhull([data.z; 0 0]))
    % neu malen lassen, Bild mit 0 initialisieren
    warning('The origin is not included in target I_C')
    data.z=[];
    data=drawD1(Empty_Button,AxesIPS,AxesIPC);
    set(Empty_Button,'UserData',data);
    return;
end

D2=[ones(size(data.z,1),1),data.z]*data.U1';
data.D=D2';
set(Empty_Button,'UserData',data);
% beide Bilder neu

data=drawD2(Empty_Button, AxesIPS,AxesIPC);

AxesD.XLimitMethod='tight';
AxesD.YLimitMethod='tight';

% D plotten
% plot(AxesD,data.D','k')
surf(AxesD,data.D)
axtoolbar(AxesD,{'rotate','restoreview'})

data.drIPC=0;
onoffButton(Button,3,[])
set(hObject,'Enable','off')
set(note,'String','')

set(Empty_Button,'UserData',data);
end



function onoffButton(handle, on, off)
for i=1:length(on)
    set(handle(on(i)),'Enable','on');
end
for i=1:length(off)
    set(handle(off(i)),'Enable','off');
end
end

function fcnbd(hObject, eventdata, parentfig, Empty_Button, Buttons, AxesIPS, AxesIPC)
data = get(Empty_Button,'UserData');
% data.y=[];
% data.z=[];
SelTyp = get(gcbf, 'SelectionType');

if strcmp(SelTyp,'alt') == 0
    pt = get(hObject, 'CurrentPoint');
    sz = get(hObject, 'Position');
    x = pt(1,1:2)./sz(1,3:4);

    posU = get(AxesIPS, 'Position');
    posV = get(AxesIPC, 'Position');

    if (posU(1) < x(1)) && (x(1) < posU(1)+posU(3)) && (posU(2) < x(2)) && (x(2) < posU(2)+posU(4)) && data.drIPS
        xl = get(AxesIPS,'xlim');
        yl = get(AxesIPS,'ylim');
        y(1)=(x(1)-posU(1))/posU(3)*abs(xl(2)-xl(1))+xl(1);
        y(2)=(x(2)-posU(2))/posU(4)*abs(yl(2)-yl(1))+yl(1);
        
        hold(AxesIPS,"on")
        plot(AxesIPS,y(1),y(2),'*','Color',[ 0.4940 0.1840 0.5560])
        hold(AxesIPS,"off")

        data.y=[data.y;y];

    elseif (posV(1) < x(1)) && (x(1) < posV(1)+posV(3)) && (posV(2) < x(2)) && (x(2) < posV(2)+posV(4)) && data.drIPC
        xl = get(AxesIPC,'xlim');
        yl = get(AxesIPC,'ylim');
        y(1)=(x(1)-posV(1))/posV(3)*abs(xl(2)-xl(1))+xl(1);
        y(2)=(x(2)-posV(2))/posV(4)*abs(yl(2)-yl(1))+yl(1);
        
        hold(AxesIPC,"on")
        plot(AxesIPC,y(1),y(2),'*','Color',[ 0.4940 0.1840 0.5560])
        hold(AxesIPC,"off")

        data.z=[data.z;y];

    end
elseif strcmp(SelTyp,'alt') == 1
    pt = get(hObject, 'CurrentPoint');
    sz = get(hObject, 'Position');
    x = pt(1,1:2)./sz(1,3:4);

    posU = get(AxesIPS, 'Position');
    posV = get(AxesIPC, 'Position');

    if (posU(1) < x(1)) && (x(1) < posU(1)+posU(3)) && (posU(2) < x(2)) && (x(2) < posU(2)+posU(4))
        tmpc = figure;
        copyobj(AxesIPS,tmpc);
        tmpa = get(tmpc,'CurrentAxes');
        set(tmpa,'position',[0.15 0.15 0.75 0.75]);
        title(tmpa, 'target I\_S');
%         xlabel(tmpa, '$T_{1,2}$','interpreter','latex');
%         ylabel(tmpa, '$T_{1,3}$','interpreter','latex');

    elseif (posV(1) < x(1)) && (x(1) < posV(1)+posV(3)) && (posV(2) < x(2)) && (x(2) < posV(2)+posV(4))
        tmpc = figure;
        copyobj(AxesIPC,tmpc);
        tmpa = get(tmpc,'CurrentAxes');
        set(tmpa,'position',[0.15 0.15 0.75 0.75]);
        title(tmpa, 'target I\_C');
%         xlabel(tmpa, '$(T^{-1})_{1,2}$','interpreter','latex');
%         ylabel(tmpa, '$(T^{-1})_{1,3}$','interpreter','latex');

    end
end

set(Empty_Button,'UserData',data);
end

%%
function data=drawIPS( Empty_Button,AxesIPS)
data = get(Empty_Button,'UserData');
D=data.D0;

[U,S,V] = svds(D,3);
for i=1:3
    if min(V(:,i)) < -max(V(:,i))
        V(:,i) = -V(:,i);
        U(:,i) = -U(:,i);
    end
end
a = diag(1./U(:,1)/S(1,1))*U(:,2:3)*S(2:3,2:3);
b = diag(1./V(:,1))*V(:,2:3);

% computation of the indices of a resp. b that define the vertices of the inner polygons
cips = convhull(a);
cipc = convhull(b);

% computation of the inner polygons
IPS = a(cips,:);
IPC = b(cipc,:);

FPS = [];
m = size(IPC,1);
for i=1:m-1
    X = [IPC(i,:); IPC(i+1,:)];
    FPS(end+1,:) = X\[-1;-1];
end
FPS(end+1,:) = FPS(1,:);

xls = [min(FPS(:,1)) max(FPS(:,1))];
yls = [min(FPS(:,2)) max(FPS(:,2))];

col = [0     0.4470     0.7410;
    0.8500    0.3250    0.0980];

plot(AxesIPS,FPS(:,1), FPS(:,2),'color',col(1,:), 'linewidth', 1)
line(AxesIPS,IPS(:,1), IPS(:,2),'color',col(2,:), 'linewidth', 1)

% hold(AxesIPS,"on")
% plot(AxesIPS,a(:,1),a(:,2),'o','color',[.7 .7 .7])
% hold(AxesIPS,"off")

set(AxesIPS,'xlim',1.05*xls,'ylim',1.05*yls)

data.FPS=FPS;
data.V0=V;
set(Empty_Button,'UserData',data);
end

%%
function data=drawD1(Empty_Button, AxesIPS,AxesIPC)
data = get(Empty_Button,'UserData');


D=data.D1';

col = [0     0.4470     0.7410;
    0.8500    0.3250    0.0980];

s = 3;
[U,S,V] = svds(D,s);
for i=1:s
    if min(V(:,i)) < -max(V(:,i))
        V(:,i) = -V(:,i);
        U(:,i) = -U(:,i);
    end
end
a = diag(1./U(:,1)/S(1,1))*U(:,2:3)*S(2:3,2:3);
b = diag(1./V(:,1))*V(:,2:3);

% computation of the indices of a resp. b that define the vertices of the inner polygons
cips = convhull(a);
cipc = convhull(b);

% computation of the inner polygons
IPS = a(cips,:);
IPC = b(cipc,:);

% computation of outer polygon (FIRPOL) for S using duality and IPC
FPS = [];
m = size(IPC,1);
for i=1:m-1
    X = [IPC(i,:); IPC(i+1,:)];
    FPS(end+1,:) = X\[-1;-1];
end
FPS(end+1,:) = FPS(1,:);

% computation of outer polygon (FIRPOL) for C using duality and IPS
FPC = [];
m = size(IPS,1);
for i=1:m-1
    X = [IPS(i,:); IPS(i+1,:)];
    FPC(end+1,:) = X\[-1;-1];
end
FPC(end+1,:) = FPC(1,:);

xls = [min(FPS(:,1)) max(FPS(:,1))];
yls = [min(FPS(:,2)) max(FPS(:,2))];
plot(AxesIPC,FPS(:,1), FPS(:,2),'color',col(2,:), 'linewidth', 1)
line(AxesIPC,IPS(:,1), IPS(:,2),'color',col(1,:), 'linewidth', 1)
set(AxesIPC,'xlim',1.05*xls,'ylim',1.05*yls)
hold(AxesIPC,"on")
plot(AxesIPC,a(:,1),a(:,2),'o','color',[.7 .7 .7])
hold(AxesIPC,"off")

xlc = [min(FPC(:,1)) max(FPC(:,1))];
ylc = [min(FPC(:,2)) max(FPC(:,2))];

plot(AxesIPS,FPC(:,1), FPC(:,2),'color',col(1,:), 'linewidth', 1)
line(AxesIPS,IPC(:,1), IPC(:,2),'color',col(2,:), 'linewidth', 1)
set(AxesIPS,'xlim',1.05*xlc,'ylim',1.05*ylc)
hold(AxesIPS,"on")
plot(AxesIPS,b(:,1),b(:,2),'o','color',[.7 .7 .7])
hold(AxesIPS,"off")

data.U1=V;
data.FPC=FPS;
set(Empty_Button,'UserData',data);
end

%%
function data=drawD2(Empty_Button, AxesIPS,AxesIPC)
data = get(Empty_Button,'UserData');


D=data.D;

col = [0     0.4470     0.7410;
    0.8500    0.3250    0.0980];

s = 3;
[U,S,V] = svds(D,s);
for i=1:s
    if min(V(:,i)) < -max(V(:,i))
        V(:,i) = -V(:,i);
        U(:,i) = -U(:,i);
    end
end
a = diag(1./U(:,1)/S(1,1))*U(:,2:3)*S(2:3,2:3);
b = diag(1./V(:,1))*V(:,2:3);

% computation of the indices of a resp. b that define the vertices of the inner polygons
cips = convhull(a);
cipc = convhull(b);

% computation of the inner polygons
IPS = a(cips,:);
IPC = b(cipc,:);

% computation of outer polygon (FIRPOL) for S using duality and IPC
FPS = [];
m = size(IPC,1);
for i=1:m-1
    X = [IPC(i,:); IPC(i+1,:)];
    FPS(end+1,:) = X\[-1;-1];
end
FPS(end+1,:) = FPS(1,:);

% computation of outer polygon (FIRPOL) for C using duality and IPS
FPC = [];
m = size(IPS,1);
for i=1:m-1
    X = [IPS(i,:); IPS(i+1,:)];
    FPC(end+1,:) = X\[-1;-1];
end
FPC(end+1,:) = FPC(1,:);

xls = [min(FPS(:,1)) max(FPS(:,1))];
yls = [min(FPS(:,2)) max(FPS(:,2))];
plot(AxesIPS,FPS(:,1), FPS(:,2),'color',col(1,:), 'linewidth', 1)
line(AxesIPS,IPS(:,1), IPS(:,2),'color',col(2,:), 'linewidth', 1)
set(AxesIPS,'xlim',1.05*xls,'ylim',1.05*yls)
hold(AxesIPS,"on")
plot(AxesIPS,a(:,1),a(:,2),'o','color',[.7 .7 .7])
hold(AxesIPS,"off")

xlc = [min(FPC(:,1)) max(FPC(:,1))];
ylc = [min(FPC(:,2)) max(FPC(:,2))];

plot(AxesIPC,FPC(:,1), FPC(:,2),'color',col(1,:), 'linewidth', 1)
line(AxesIPC,IPC(:,1), IPC(:,2),'color',col(2,:), 'linewidth', 1)
set(AxesIPC,'xlim',1.05*xlc,'ylim',1.05*ylc)
hold(AxesIPC,"on")
plot(AxesIPC,b(:,1),b(:,2),'o','color',[.7 .7 .7])
hold(AxesIPC,"off")

set(Empty_Button,'UserData',data);
end
