% 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.
% The optimization problem is located in the function obj.calcMatrix


%%
classdef recover_matrix_opt < handle
    properties
        % Notes:
        % Generally the following numbers are used for the figures
        % 1 left side, 2 middle, 3 right side
        % and for the plots
        % 1 dual F, 2 dual I, 3 primal F, 4 primal I,
        % 5 delete in dual, 6 delete in primal
        
        basicfsz = 10
        sty = 'Helvetica'
        colorMat
        
        D       % Nonnegative Matrix
        F1      % Primal FIRPOL = drawn points
        kF1     % Primal FIRPOL coefficients
        I1      % Primal INNPOL = primal data points
        kI1     % Primal INNPOL coefficients
        F2      % Dual FIRPOL = drawn points
        kF2     % Dual FIRPOL coefficients
        I2      % Dual INNPOL = dual data points
        kI2     % Dual INNPOL coefficients
        Res     % Here the checked results are saved 
        
        hRootFigure     matlab.ui.Figure
        
        hLeftPanel
        hLeftAxes
        hButtonsLeft
        hLeftPlots
        
        hMiddlePanel
        hMiddleAxes
        hButtonsMiddle
        hMiddlePlots
        
        hRightPanel
        hRightAxes
        hButtonsRight
    end
    
     %% Ordinary methods
    methods
        %% drawData: Constructor
        function obj = recover_matrix_opt(version,figNR)
            if nargin == 0 % Started as standalone
                % Create the root figure
                obj.hRootFigure = figure('Visible','off','outerposition',[50 50 1820 810],...
                    'NumberTitle','off','MenuBar','none','Toolbar','none','color','white',...
                    'Name','Geometric data construction');
            else % Started in Facpack
                % Create the root figure
                obj.hRootFigure = figure(figNR);
                set(obj.hRootFigure,'Visible','off','outerposition',[50 50 1820 810],...
                    'NumberTitle','off','MenuBar','none','Toolbar','none','color','white',...
                    'Name',['FACPACK (Geometric data construction, version ',num2str(version,'%2.2f'), '-alpha)']);
            end
            createComp(obj)
            obj.F1 = [];      % Primal FIRPOL = drawn points
            obj.kF1 = [];     % Primal FIRPOL coefficients
            obj.I1 = [];      % Primal INNPOL = primal data points
            obj.kI1 = [];     % Primal INNPOL coefficients
            obj.F2 = [];      % Dual FIRPOL = drawn points
            obj.kF2 = [];     % Dual FIRPOL coefficients
            obj.I2 = [];      % Dual INNPOL = dual data points
            obj.kI2 = [];     % Dual INNPOL coefficients
        end
        %% createComp: Create GUI components
        function createComp(obj)
            % Left side
            obj.hLeftPanel = uipanel('Title','U-Space',...
                'BackgroundColor','w','FontUnits','pixel','Units','pixel','FontSize',obj.basicfsz*1.2,...
                'FontName',obj.sty,'Position',[5 5 600 750],'Visible','on','Units','normalized');
            obj.hLeftAxes = axes('Parent',obj.hLeftPanel,'FontSize',obj.basicfsz,...
                'FontName',obj.sty,'units','pixel','Position',[50 50 500 500],'Units','normalized');
            xlim(obj.hLeftAxes,[-2 2]);
            ylim(obj.hLeftAxes,[-2 2]);
            box(obj.hLeftAxes,'on');
            hold(obj.hLeftAxes,'on');
            axtoolbar(obj.hLeftAxes,'default');
            obj.hLeftPlots.f = plot(obj.hLeftAxes,[0],[0],'w');
            obj.hLeftPlots.i = plot(obj.hLeftAxes,[0],[0],'w');
            obj.hLeftPlots.dp1 = plot(obj.hLeftAxes,[0],[0],'w');
            obj.hLeftPlots.dp2 = plot(obj.hLeftAxes,[0],[0],'w');
            obj.hLeftPlots.fnew = plot(obj.hLeftAxes,[0],[0],'w');
            obj.hLeftPlots.inew = plot(obj.hLeftAxes,[0],[0],'w');
            plot(obj.hLeftAxes,[0],[0],'Color',[.7 .7 .7],'Marker','+');
            obj.hButtonsLeft.hButton1 = uicontrol('Style','togglebutton','string','Draw F','Parent',obj.hLeftPanel,'Position',...
                [50 630 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gDrawButton(hObj,event,1)},...
                'Tooltip','F: outer polygon, FIRPOL. Only convex hull is used.');
            obj.hButtonsLeft.hButton2 = uicontrol('Style','pushbutton','string','Load F','Parent',obj.hLeftPanel,'Position',...
                [200 630 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gLoad(hObj,event,1)},...
                'Tooltip','(optional) F: outer polygon, FIRPOL. Only convex hull is used.');
            obj.hButtonsLeft.hButton3 = uicontrol('Style','pushbutton','string','Reset F','Parent',obj.hLeftPanel,'Position',...
                [350 630 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gReset(hObj,event,1)},...
                'Tooltip','Also resets the dual I with all its points');
            obj.hButtonsLeft.hButton4 = uicontrol('Style','togglebutton','string','Draw I','Parent',obj.hLeftPanel,'Position',...
                [50 580 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gDrawButton(hObj,event,2)},...
                'Tooltip','I: inner polygon, INNPOL. All points are used.');
            obj.hButtonsLeft.hButton5 = uicontrol('Style','pushbutton','string','Load I','Parent',obj.hLeftPanel,'Position',...
                [200 580 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gLoad(hObj,event,2)},...
                'Tooltip','(optional) I: inner polygon, INNPOL. All points are used.');
            obj.hButtonsLeft.hButton6 = uicontrol('Style','togglebutton','string','Delete Pts','Parent',obj.hLeftPanel,'Position',...
                [350 580 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gDrawButton(hObj,event,5)},...
                'Tooltip','Click on points to delete them.');
            % Middle
            obj.hMiddlePanel = uipanel('Title','V-Space',...
                'BackgroundColor','w','FontUnits','pixel','Units','pixel','FontSize',obj.basicfsz*1.2,...
                'FontName',obj.sty,'Position',[610 5 600 750],'Visible','on','Units','normalized');
            obj.hMiddleAxes = axes('Parent',obj.hMiddlePanel,'FontSize',obj.basicfsz,...
                'FontName',obj.sty,'units','pixel','Position',[50 50 500 500],'Units','normalized');
            xlim(obj.hMiddleAxes,[-2 2]);
            ylim(obj.hMiddleAxes,[-2 2]);
            box(obj.hMiddleAxes,'on');
            hold(obj.hMiddleAxes,'on');
            axtoolbar(obj.hMiddleAxes,'default');
            obj.hMiddlePlots.f = plot(obj.hMiddleAxes,[0],[0],'w');
            obj.hMiddlePlots.i = plot(obj.hMiddleAxes,[0],[0],'w');
            obj.hMiddlePlots.dp3 = plot(obj.hMiddleAxes,[0],[0],'w');
            obj.hMiddlePlots.dp4 = plot(obj.hMiddleAxes,[0],[0],'w');
            obj.hMiddlePlots.fnew = plot(obj.hMiddleAxes,[0],[0],'w');
            obj.hMiddlePlots.inew = plot(obj.hMiddleAxes,[0],[0],'w');
            plot(obj.hMiddleAxes,[0],[0],'Color',[.7 .7 .7],'Marker','+');
            obj.hButtonsMiddle.hButton1 = uicontrol('Style','togglebutton','string','Draw F','Parent',obj.hMiddlePanel,'Position',...
                [50 630 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gDrawButton(hObj,event,3)},...
                'Tooltip','F: outer polygon, FIRPOL. Only convex hull is used.');
            obj.hButtonsMiddle.hButton2 = uicontrol('Style','pushbutton','string','Load F','Parent',obj.hMiddlePanel,'Position',...
                [200 630 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gLoad(hObj,event,3)},...
                'Tooltip','(optional) F: outer polygon, FIRPOL. Only convex hull is used.');
            obj.hButtonsMiddle.hButton3 = uicontrol('Style','pushbutton','string','Reset F','Parent',obj.hMiddlePanel,'Position',...
                [350 630 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gReset(hObj,event,3)},...
                'Tooltip','Also resets the dual I with all its points');
            obj.hButtonsMiddle.hButton4 = uicontrol('Style','togglebutton','string','Draw I','Parent',obj.hMiddlePanel,'Position',...
                [50 580 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gDrawButton(hObj,event,4)},...
                'Tooltip','I: inner polygon, INNPOL. All points are used.');
            obj.hButtonsMiddle.hButton5 = uicontrol('Style','pushbutton','string','Load I','Parent',obj.hMiddlePanel,'Position',...
                [200 580 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gLoad(hObj,event,4)},...
                'Tooltip','(optional) I: inner polygon, INNPOL. All points are used.');
            obj.hButtonsMiddle.hButton6 = uicontrol('Style','togglebutton','string','Delete Pts','Parent',obj.hMiddlePanel,'Position',...
                [350 580 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gDrawButton(hObj,event,6)},...
                'Tooltip','Click on points to delete them.');
            % Right side
            obj.hRightPanel = uipanel('Title','Data matrix',...
                'BackgroundColor','w','FontUnits','pixel','Units','pixel','FontSize',obj.basicfsz*1.2,...
                'FontName',obj.sty,'Position',[1215 5 600 750],'Visible','on','Units','normalized');
            obj.hRightAxes = axes('Parent',obj.hRightPanel,'FontSize',obj.basicfsz,...
                'FontName',obj.sty,'units','pixel','Position',[50 50 500 500],'Units','normalized');
            obj.hButtonsRight.hButton71 = uicontrol('Style','text','string','D: small values','Parent',obj.hRightPanel,'Position',...
                [175 580 60 40],'Units','normalized',...
                'Tooltip','Increase to possibly avoid lots of near 0 entries in D. 0 (default), >1 (no optimal solution possible). Coefficient for the minimal values of first singular vectors.');
            obj.hButtonsRight.hButton7 = uicontrol('Style','edit','string','0','Parent',obj.hRightPanel,'Position',...
                [235 580 60 40],'Units','normalized',...
                'Tooltip','Coefficient for the minimal values of first singular vectors: 0 (default), >1 (no optimal solution possible).');
            obj.hButtonsRight.hButton8 = uicontrol('Style','popupmenu','string',...
                {'fminunc','fminsearch','fsolve','lsqnonlin','patternsearch','simulannealbnd','fminmax','particleswarm','ga'},...
                'Parent',obj.hRightPanel,'Position',[175 630 120 40],'Units','normalized',...
                'Tooltip','Optimization routine for the problem (feel free to adjust options).');
            obj.hButtonsRight.hButton3 = uicontrol('Style','popupmenu','string',{'Const start vect.','Rand start vect.'},...
                'Parent',obj.hRightPanel,'Position',[175 680 120 40],'Units','normalized',...
                'Tooltip','Starting vector of the optimization problem as normed ones or normed random vector.');
            obj.hButtonsRight.hButton1 = uicontrol('Style','pushbutton','string','Accept results','Parent',obj.hRightPanel,'Position',...
                [425 580 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gAcceptLoaded(hObj,event)},...
                'Tooltip','The original F and I is overwritten with the calculated points.');
            obj.hButtonsRight.hButton2 = uicontrol('Style','pushbutton','string','Calculate','FontWeight','bold','Parent',...
                obj.hRightPanel,'Position',[300 680 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gCalc(hObj,event)},...
                'Tooltip','Calculate a matrix D for the given F and I.');
            obj.hButtonsRight.hButton9 = uicontrol('Style','pushbutton','string','Load D','Parent',obj.hRightPanel,'Position',...
                [50 680 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gLoadD(hObj,event)},...
                'Tooltip','(optional) Load a matrix D and use its F and I. Use this to adjust I and F of a given dataset.');
            obj.hButtonsRight.hButton4 = uicontrol('Style','pushbutton','string','Export','Parent',obj.hRightPanel,'Position',...
                [50 580 120 40],'Units','normalized','Callback',{@(hObj,event) obj.gSave(hObj,event)},...
                'Tooltip','Save D and last accepted F and I.');
            obj.hButtonsRight.hButton5 = uicontrol('Style','text','string','Residual: ','Parent',obj.hRightPanel,'Position',...
                [300 580 120 40],'Units','normalized','HorizontalAlignment','left',...
                'Tooltip','Residual of the chosen solver. Be cautious if the value is large.');
            obj.hButtonsRight.hButton6 = uicontrol('Style','text','string','Min(D): ','Parent',obj.hRightPanel,'Position',...
                [300 630 120 40],'Units','normalized','HorizontalAlignment','left',...
                'Tooltip','Smallest value in D. Be cautious if the value is negative.');
            % Figure callback
            obj.hLeftAxes.ButtonDownFcn = {@(hObj,event) obj.gAxesCallback(hObj,event,1)};
            obj.hMiddleAxes.ButtonDownFcn = {@(hObj,event) obj.gAxesCallback(hObj,event,2)};
            obj.hRightAxes.ButtonDownFcn = {@(hObj,event) obj.gAxesCallback(hObj,event,3)};
            % Plot callback
            obj.hLeftPlots.dp1.ButtonDownFcn = {@(hObj,event) obj.gPlotsCallback(hObj,event,1)};
            obj.hLeftPlots.dp2.ButtonDownFcn = {@(hObj,event) obj.gPlotsCallback(hObj,event,2)};
            obj.hMiddlePlots.dp3.ButtonDownFcn = {@(hObj,event) obj.gPlotsCallback(hObj,event,3)};
            obj.hMiddlePlots.dp4.ButtonDownFcn = {@(hObj,event) obj.gPlotsCallback(hObj,event,4)};
            % Set position for the current screen size (comment this if
            % debugging)
            scrsz = get(0,'ScreenSize');
            set(obj.hRootFigure,'Position',[1/20*scrsz(3:4) 9/10*scrsz(3) min(1/2.5*9/10*scrsz(3),9/10*scrsz(4))]);
            % Show everything
            set(obj.hRootFigure,'Visible','on');
        end
        %% gCalc: Start calculation
        function gCalc(obj, hObj, event)
            % Turn buttons off
            obj.gTurnButtonsOff([1 1 1 1 1 1]);
            % Check whether there are enough points for a calculation
            if size(obj.F1,1)<3 || size(obj.I1,1)<3
                errordlg('Use inner and outer boundary with at least 3 points!','Error');
            else
                % Calculate matrix
                optimError = obj.calcMatrix(obj.hButtonsRight.hButton3.Value,obj.hButtonsRight.hButton8.Value,...
                    str2double(obj.hButtonsRight.hButton7.String));
                % Normalize matrix
                obj.D = obj.D./max(max(obj.D));
                % Show results
                obj.hButtonsRight.hButton5.String = ['Residual: ',num2str(optimError)];
                obj.hButtonsRight.hButton6.String = ['min(D): ',num2str(min(min(obj.D)))];
                surf(obj.hRightAxes,obj.D');
                obj.hRightAxes.ButtonDownFcn = {@(hObj,event) obj.gAxesCallback(hObj,event,3)};
                axtoolbar(obj.hRightAxes,'default');
                % Recalculate boundaries from D, show them
                obj.checkResults()
                % Plot new results for comparison
                xlim(obj.hLeftAxes,'auto');
                ylim(obj.hLeftAxes,'auto');
                obj.hLeftPlots.fnew.XData=obj.Res.F2(obj.Res.kF2,1);
                obj.hLeftPlots.fnew.YData=obj.Res.F2(obj.Res.kF2,2);
                obj.hLeftPlots.fnew.Color = [1 0 1];
                obj.hLeftPlots.fnew.LineStyle='-';
                obj.hLeftPlots.inew.XData=obj.Res.I2(obj.Res.kI2,1);
                obj.hLeftPlots.inew.YData=obj.Res.I2(obj.Res.kI2,2);
                obj.hLeftPlots.inew.Color = [0 1 1];
                obj.hLeftPlots.inew.LineStyle='-';
                xlim(obj.hMiddleAxes,'auto');
                ylim(obj.hMiddleAxes,'auto');
                obj.hMiddlePlots.fnew.XData=obj.Res.F1(obj.Res.kF1,1);
                obj.hMiddlePlots.fnew.YData=obj.Res.F1(obj.Res.kF1,2);
                obj.hMiddlePlots.fnew.Color = [0 1 1];
                obj.hMiddlePlots.fnew.LineStyle='-';
                obj.hMiddlePlots.inew.XData=obj.Res.I1(obj.Res.kI1,1);
                obj.hMiddlePlots.inew.YData=obj.Res.I1(obj.Res.kI1,2);
                obj.hMiddlePlots.inew.Color = [1 0 1];
                obj.hMiddlePlots.inew.LineStyle='-';
            end
        end
        %% gAxesCallback: Add points, plot them, open plots in new window
        function gAxesCallback(obj, hObj, event, axNr)
            if event.Button==1
                if axNr==1
                    if obj.hButtonsLeft.hButton1.Value==1
                        obj.F2 = [obj.F2;obj.hLeftAxes.CurrentPoint(1, 1:2)];
                        obj.hLeftPlots.dp1.XData=obj.F2(:,1);
                        obj.hLeftPlots.dp1.YData=obj.F2(:,2);
                        obj.hLeftPlots.dp1.Color = [1 0 0];
                        obj.hLeftPlots.dp1.LineStyle='none';
                        obj.hLeftPlots.dp1.Marker='*';
                    elseif obj.hButtonsLeft.hButton4.Value==1
                        obj.I2 = [obj.I2;obj.hLeftAxes.CurrentPoint(1, 1:2)];
                        obj.hLeftPlots.dp2.XData=obj.I2(:,1);
                        obj.hLeftPlots.dp2.YData=obj.I2(:,2);
                        obj.hLeftPlots.dp2.Color = [0 0 1];
                        obj.hLeftPlots.dp2.LineStyle='none';
                        obj.hLeftPlots.dp2.Marker='*';
                    end
                elseif axNr==2
                    if obj.hButtonsMiddle.hButton1.Value==1
                        obj.F1 = [obj.F1;obj.hMiddleAxes.CurrentPoint(1, 1:2)];
                        obj.hMiddlePlots.dp3.XData=obj.F1(:,1);
                        obj.hMiddlePlots.dp3.YData=obj.F1(:,2);
                        obj.hMiddlePlots.dp3.Color = [0 0 1];
                        obj.hMiddlePlots.dp3.LineStyle='none';
                        obj.hMiddlePlots.dp3.Marker='*';
                    elseif obj.hButtonsMiddle.hButton4.Value==1
                        obj.I1 = [obj.I1;obj.hMiddleAxes.CurrentPoint(1, 1:2)];
                        obj.hMiddlePlots.dp4.XData=obj.I1(:,1);
                        obj.hMiddlePlots.dp4.YData=obj.I1(:,2);
                        obj.hMiddlePlots.dp4.Color = [1 0 0];
                        obj.hMiddlePlots.dp4.LineStyle='none';
                        obj.hMiddlePlots.dp4.Marker='*';
                    end
                end
            % This opens the clicked plot in a new window
            elseif event.Button==3
                fignew = figure();
                if axNr==1
                    temp_plot = copyobj(obj.hLeftAxes,fignew);
                    temp_plot.Position = [0.1300 0.1100 0.7750 0.8150];
                elseif axNr==2
                    temp_plot = copyobj(obj.hMiddleAxes,fignew);
                    temp_plot.Position = [0.1300 0.1100 0.7750 0.8150];
                elseif axNr==3
                    temp_plot = copyobj(obj.hRightAxes,fignew);
                    temp_plot.Position = [0.1300 0.1100 0.7750 0.8150];

                end
            end
        end
        %% gPlotsCallback: To remove points
        function gPlotsCallback(obj, hObj, event, boundary)
            if event.Button==1
                Pt = event.IntersectionPoint(1:2);
                switch boundary
                    case 1
                        id = ismembertol(obj.F2(:,1),Pt(1),1e-5)&ismembertol(obj.F2(:,2),Pt(2),1e-5);
                        obj.F2 = obj.F2(~id,:);
                        obj.hLeftPlots.dp1.XData=obj.F2(:,1);
                        obj.hLeftPlots.dp1.YData=obj.F2(:,2);
                    case 2
                        id = ismembertol(obj.I2(:,1),Pt(1),1e-5)&ismembertol(obj.I2(:,2),Pt(2),1e-5);
                        obj.I2 = obj.I2(~id,:);
                        obj.hLeftPlots.dp2.XData=obj.I2(:,1);
                        obj.hLeftPlots.dp2.YData=obj.I2(:,2);
                    case 3
                        id = ismembertol(obj.F1(:,1),Pt(1),1e-5)&ismembertol(obj.F1(:,2),Pt(2),1e-5);
                        obj.F1 = obj.F1(~id,:);
                        obj.hMiddlePlots.dp3.XData=obj.F1(:,1);
                        obj.hMiddlePlots.dp3.YData=obj.F1(:,2);
                    case 4
                        id = ismembertol(obj.I1(:,1),Pt(1),1e-5)&ismembertol(obj.I1(:,2),Pt(2),1e-5);
                        obj.I1 = obj.I1(~id,:);
                        obj.hMiddlePlots.dp4.XData=obj.I1(:,1);
                        obj.hMiddlePlots.dp4.YData=obj.I1(:,2);
                end
            end
        end
        %% gDrawButton: Turn draw buttons off on click
        % and stop the scaling of axes to allow for drawing
        function gDrawButton(obj, hObj, event, boundary)
            if hObj.Value==1
                switch boundary
                    case 1
                        obj.gTurnButtonsOff([0 1 1 1 1 1])
                        xlim(obj.hLeftAxes);
                        xlim(obj.hLeftAxes,'manual');
                        ylim(obj.hLeftAxes,'manual');
                    case 2
                        obj.gTurnButtonsOff([1 0 1 1 1 1])
                        xlim(obj.hLeftAxes);
                        xlim(obj.hLeftAxes,'manual');
                        ylim(obj.hLeftAxes,'manual');
                    case 3
                        obj.gTurnButtonsOff([1 1 0 1 1 1])
                        xlim(obj.hMiddleAxes);
                        xlim(obj.hMiddleAxes,'manual');
                        ylim(obj.hMiddleAxes,'manual');
                    case 4
                        obj.gTurnButtonsOff([1 1 1 0 1 1])
                        xlim(obj.hMiddleAxes);
                        xlim(obj.hMiddleAxes,'manual');
                        ylim(obj.hMiddleAxes,'manual');
                    case 5
                        obj.gTurnButtonsOff([1 1 1 1 0 1])
                        xlim(obj.hLeftAxes);
                        xlim(obj.hLeftAxes,'manual');
                        ylim(obj.hLeftAxes,'manual');
                    case 6
                        obj.gTurnButtonsOff([1 1 1 1 1 0])
                        xlim(obj.hMiddleAxes);
                        xlim(obj.hMiddleAxes,'manual');
                        ylim(obj.hMiddleAxes,'manual');
                end
            end
            if hObj.Value==0
                switch boundary
                    case 1
                        obj.hButtonsLeft.hButton1.Value = 1;
                        obj.gTurnButtonsOff([1 0 0 0 0 0]);
                        xlim(obj.hLeftAxes);
                        xlim(obj.hLeftAxes,'manual');
                        ylim(obj.hLeftAxes,'manual');
                    case 2
                        obj.hButtonsLeft.hButton4.Value = 1;
                        obj.gTurnButtonsOff([0 1 0 0 0 0]);
                        xlim(obj.hLeftAxes);
                        xlim(obj.hLeftAxes,'manual');
                        ylim(obj.hLeftAxes,'manual');
                    case 3
                        obj.hButtonsMiddle.hButton1.Value = 1;
                        obj.gTurnButtonsOff([0 0 1 0 0 0]);
                        xlim(obj.hMiddleAxes);
                        xlim(obj.hMiddleAxes,'manual');
                        ylim(obj.hMiddleAxes,'manual');
                    case 4
                        obj.hButtonsMiddle.hButton4.Value = 1;
                        obj.gTurnButtonsOff([0 0 0 1 0 0]);
                        xlim(obj.hMiddleAxes);
                        xlim(obj.hMiddleAxes,'manual');
                        ylim(obj.hMiddleAxes,'manual');
                    case 5
                        obj.hButtonsLeft.hButton6.Value = 1;
                        obj.gTurnButtonsOff([0 0 0 0 1 0]);
                        xlim(obj.hLeftAxes);
                        xlim(obj.hLeftAxes,'manual');
                        ylim(obj.hLeftAxes,'manual');
                    case 6
                        obj.hButtonsMiddle.hButton6.Value = 1;
                        obj.gTurnButtonsOff([0 0 0 0 0 1]);
                        xlim(obj.hMiddleAxes);
                        xlim(obj.hMiddleAxes,'manual');
                        ylim(obj.hMiddleAxes,'manual');
                end
            end
        end
        %% gTurnButtonsOff: Turn buttons off, if they were on, 
        % then calculate convex hull and dual plot again
        function gTurnButtonsOff(obj,ButtonNr)
            if ButtonNr(1)==1
                if obj.hButtonsLeft.hButton1.Value==1
                    obj.hButtonsLeft.hButton1.Value = 0;
                    if size(obj.F2,1)>2
                        obj.kF2 = convhull(obj.F2);
                        if inpolygon(0,0,obj.F2(obj.kF2,1),obj.F2(obj.kF2,2))==0
                            errordlg('Origin (0,0) not included in the boundary!','Error');
                        end
                        obj.hLeftPlots.f.XData=obj.F2(obj.kF2,1);
                        obj.hLeftPlots.f.YData=obj.F2(obj.kF2,2);
                        obj.hLeftPlots.f.Color = [1 0 0];
                        obj.hLeftPlots.f.LineStyle='-';
                        % Dual plot
                        xlim(obj.hMiddleAxes,'auto');
                        ylim(obj.hMiddleAxes,'auto');
                        [obj.I1,obj.kI1] = obj.calcDual(obj.F2,obj.kF2,obj.I1,obj.kI1);
                        obj.hMiddlePlots.dp4.XData=obj.I1(:,1);
                        obj.hMiddlePlots.dp4.YData=obj.I1(:,2);
                        obj.hMiddlePlots.dp4.Color = [1 0 0];
                        obj.hMiddlePlots.dp4.LineStyle='none';
                        obj.hMiddlePlots.dp4.Marker='*';
                        obj.hMiddlePlots.i.XData=obj.I1(obj.kI1,1);
                        obj.hMiddlePlots.i.YData=obj.I1(obj.kI1,2);
                        obj.hMiddlePlots.i.Color = [1 0 0];
                        obj.hMiddlePlots.i.LineStyle='-';
                    end
                end
            end
            if ButtonNr(2)==1
                if obj.hButtonsLeft.hButton4.Value==1
                    obj.hButtonsLeft.hButton4.Value = 0;
                    if size(obj.I2,1)>2
                        obj.kI2 = convhull(obj.I2);
                        if inpolygon(0,0,obj.I2(obj.kI2,1),obj.I2(obj.kI2,2))==0
                            errordlg('Origin (0,0) not included in the boundary!','Error');
                        end
                        obj.hLeftPlots.i.XData=obj.I2(obj.kI2,1);
                        obj.hLeftPlots.i.YData=obj.I2(obj.kI2,2);
                        obj.hLeftPlots.i.Color = [0 0 1];
                        obj.hLeftPlots.i.LineStyle='-';
                        % Dual plot
                        xlim(obj.hMiddleAxes,'auto');
                        ylim(obj.hMiddleAxes,'auto');
                        [obj.F1,obj.kF1] = obj.calcDual(obj.I2,obj.kI2);
                        obj.hMiddlePlots.dp3.XData=obj.F1(:,1);
                        obj.hMiddlePlots.dp3.YData=obj.F1(:,2);
                        obj.hMiddlePlots.dp3.Color = [0 0 1];
                        obj.hMiddlePlots.dp3.LineStyle='none';
                        obj.hMiddlePlots.dp3.Marker='*';
                        obj.hMiddlePlots.f.XData=obj.F1(obj.kF1,1);
                        obj.hMiddlePlots.f.YData=obj.F1(obj.kF1,2);
                        obj.hMiddlePlots.f.Color = [0 0 1];
                        obj.hMiddlePlots.f.LineStyle='-';
                    end
                end
            end
            if ButtonNr(3)==1
                if obj.hButtonsMiddle.hButton1.Value==1
                    obj.hButtonsMiddle.hButton1.Value = 0;
                    if size(obj.F1,1)>2
                        obj.kF1 = convhull(obj.F1);
                        if inpolygon(0,0,obj.F1(obj.kF1,1),obj.F1(obj.kF1,2))==0
                            errordlg('Origin (0,0) not included in the boundary!','Error');
                        end
                        obj.hMiddlePlots.f.XData=obj.F1(obj.kF1,1);
                        obj.hMiddlePlots.f.YData=obj.F1(obj.kF1,2);
                        obj.hMiddlePlots.f.Color = [0 0 1];
                        obj.hMiddlePlots.f.LineStyle='-';
                        % Dual plot
                        xlim(obj.hLeftAxes,'auto');
                        ylim(obj.hLeftAxes,'auto');
                        [obj.I2,obj.kI2] = obj.calcDual(obj.F1,obj.kF1,obj.I2,obj.kI2);
                        obj.hLeftPlots.dp2.XData=obj.I2(:,1);
                        obj.hLeftPlots.dp2.YData=obj.I2(:,2);
                        obj.hLeftPlots.dp2.Color = [0 0 1];
                        obj.hLeftPlots.dp2.LineStyle='none';
                        obj.hLeftPlots.dp2.Marker='*';
                        obj.hLeftPlots.i.XData=obj.I2(obj.kI2,1);
                        obj.hLeftPlots.i.YData=obj.I2(obj.kI2,2);
                        obj.hLeftPlots.i.Color = [0 0 1];
                        obj.hLeftPlots.i.LineStyle='-';
                    end
                end
            end
            if ButtonNr(4)==1
                if obj.hButtonsMiddle.hButton4.Value==1
                    obj.hButtonsMiddle.hButton4.Value = 0;
                    if size(obj.I1,1)>2
                        obj.kI1 = convhull(obj.I1);
                        if inpolygon(0,0,obj.I1(obj.kI1,1),obj.I1(obj.kI1,2))==0
                            errordlg('Origin (0,0) not included in the boundary!','Error');
                        end
                        obj.hMiddlePlots.i.XData=obj.I1(obj.kI1,1);
                        obj.hMiddlePlots.i.YData=obj.I1(obj.kI1,2);
                        obj.hMiddlePlots.i.Color = [1 0 0];
                        obj.hMiddlePlots.i.LineStyle='-';
                        % Dual plot
                        xlim(obj.hLeftAxes,'auto');
                        ylim(obj.hLeftAxes,'auto');
                        [obj.F2,obj.kF2] = obj.calcDual(obj.I1,obj.kI1);
                        obj.hLeftPlots.dp1.XData=obj.F2(:,1);
                        obj.hLeftPlots.dp1.YData=obj.F2(:,2);
                        obj.hLeftPlots.dp1.Color = [1 0 0];
                        obj.hLeftPlots.dp1.LineStyle='none';
                        obj.hLeftPlots.dp1.Marker='*';
                        obj.hLeftPlots.f.XData=obj.F2(obj.kF2,1);
                        obj.hLeftPlots.f.YData=obj.F2(obj.kF2,2);
                        obj.hLeftPlots.f.Color = [1 0 0];
                        obj.hLeftPlots.f.LineStyle='-';
                    end
                end
            end
            if ButtonNr(5)==1
                if obj.hButtonsLeft.hButton6.Value==1
                    obj.hButtonsLeft.hButton6.Value = 0;
                    % F2
                    if size(obj.F2,1)>2
                        obj.kF2 = convhull(obj.F2);
                        if inpolygon(0,0,obj.F2(obj.kF2,1),obj.F2(obj.kF2,2))==0
                            errordlg('Origin (0,0) not included in the boundary!','Error');
                        end
                        obj.hLeftPlots.f.XData=obj.F2(obj.kF2,1);
                        obj.hLeftPlots.f.YData=obj.F2(obj.kF2,2);
                        obj.hLeftPlots.f.Color = [1 0 0];
                        obj.hLeftPlots.f.LineStyle='-';
                        % Dual plot
                        xlim(obj.hMiddleAxes,'auto');
                        ylim(obj.hMiddleAxes,'auto');
                        [obj.I1,obj.kI1] = obj.calcDual(obj.F2,obj.kF2,obj.I1,obj.kI1);
                        obj.hMiddlePlots.dp4.XData=obj.I1(:,1);
                        obj.hMiddlePlots.dp4.YData=obj.I1(:,2);
                        obj.hMiddlePlots.dp4.Color = [1 0 0];
                        obj.hMiddlePlots.dp4.LineStyle='none';
                        obj.hMiddlePlots.dp4.Marker='*';
                        obj.hMiddlePlots.i.XData=obj.I1(obj.kI1,1);
                        obj.hMiddlePlots.i.YData=obj.I1(obj.kI1,2);
                        obj.hMiddlePlots.i.Color = [1 0 0];
                        obj.hMiddlePlots.i.LineStyle='-';
                    end
                    % I2
                    if size(obj.I2,1)>2
                        obj.kI2 = convhull(obj.I2);
                        if inpolygon(0,0,obj.I2(obj.kI2,1),obj.I2(obj.kI2,2))==0
                            errordlg('Origin (0,0) not included in the boundary!','Error');
                        end
                        obj.hLeftPlots.i.XData=obj.I2(obj.kI2,1);
                        obj.hLeftPlots.i.YData=obj.I2(obj.kI2,2);
                        obj.hLeftPlots.i.Color = [0 0 1];
                        obj.hLeftPlots.i.LineStyle='-';
                        % Dual plot
                        xlim(obj.hMiddleAxes,'auto');
                        ylim(obj.hMiddleAxes,'auto');
                        [obj.F1,obj.kF1] = obj.calcDual(obj.I2,obj.kI2);
                        obj.hMiddlePlots.dp3.XData=obj.F1(:,1);
                        obj.hMiddlePlots.dp3.YData=obj.F1(:,2);
                        obj.hMiddlePlots.dp3.Color = [0 0 1];
                        obj.hMiddlePlots.dp3.LineStyle='none';
                        obj.hMiddlePlots.dp3.Marker='*';
                        obj.hMiddlePlots.f.XData=obj.F1(obj.kF1,1);
                        obj.hMiddlePlots.f.YData=obj.F1(obj.kF1,2);
                        obj.hMiddlePlots.f.Color = [0 0 1];
                        obj.hMiddlePlots.f.LineStyle='-';
                    end
                end
            end
            if ButtonNr(6)==1
                if obj.hButtonsMiddle.hButton6.Value==1
                    obj.hButtonsMiddle.hButton6.Value = 0;
                    % F1
                    if size(obj.F1,1)>2
                        obj.kF1 = convhull(obj.F1);
                        if inpolygon(0,0,obj.F1(obj.kF1,1),obj.F1(obj.kF1,2))==0
                            errordlg('Origin (0,0) not included in the boundary!','Error');
                        end
                        obj.hMiddlePlots.f.XData=obj.F1(obj.kF1,1);
                        obj.hMiddlePlots.f.YData=obj.F1(obj.kF1,2);
                        obj.hMiddlePlots.f.Color = [0 0 1];
                        obj.hMiddlePlots.f.LineStyle='-';
                        % Dual plot
                        xlim(obj.hLeftAxes,'auto');
                        ylim(obj.hLeftAxes,'auto');
                        [obj.I2,obj.kI2] = obj.calcDual(obj.F1,obj.kF1,obj.I2,obj.kI2);
                        obj.hLeftPlots.dp2.XData=obj.I2(:,1);
                        obj.hLeftPlots.dp2.YData=obj.I2(:,2);
                        obj.hLeftPlots.dp2.Color = [0 0 1];
                        obj.hLeftPlots.dp2.LineStyle='none';
                        obj.hLeftPlots.dp2.Marker='*';
                        obj.hLeftPlots.i.XData=obj.I2(obj.kI2,1);
                        obj.hLeftPlots.i.YData=obj.I2(obj.kI2,2);
                        obj.hLeftPlots.i.Color = [0 0 1];
                        obj.hLeftPlots.i.LineStyle='-';
                    end
                    % I1
                    if size(obj.I1,1)>2
                        obj.kI1 = convhull(obj.I1);
                        if inpolygon(0,0,obj.I1(obj.kI1,1),obj.I1(obj.kI1,2))==0
                            errordlg('Origin (0,0) not included in the boundary!','Error');
                        end
                        obj.hMiddlePlots.i.XData=obj.I1(obj.kI1,1);
                        obj.hMiddlePlots.i.YData=obj.I1(obj.kI1,2);
                        obj.hMiddlePlots.i.Color = [1 0 0];
                        obj.hMiddlePlots.i.LineStyle='-';
                        % Dual plot
                        xlim(obj.hLeftAxes,'auto');
                        ylim(obj.hLeftAxes,'auto');
                        [obj.F2,obj.kF2] = obj.calcDual(obj.I1,obj.kI1);
                        obj.hLeftPlots.dp1.XData=obj.F2(:,1);
                        obj.hLeftPlots.dp1.YData=obj.F2(:,2);
                        obj.hLeftPlots.dp1.Color = [1 0 0];
                        obj.hLeftPlots.dp1.LineStyle='none';
                        obj.hLeftPlots.dp1.Marker='*';
                        obj.hLeftPlots.f.XData=obj.F2(obj.kF2,1);
                        obj.hLeftPlots.f.YData=obj.F2(obj.kF2,2);
                        obj.hLeftPlots.f.Color = [1 0 0];
                        obj.hLeftPlots.f.LineStyle='-';
                    end
                end
            end
        end
        %% gReset: Reset a boundary
        function gReset(obj, hObj, event, boundary)
            % Remove previous results            
            obj.hLeftPlots.fnew.XData=[0];
            obj.hLeftPlots.fnew.YData=[0];
            obj.hLeftPlots.fnew.Color = [0 0 0];
            obj.hLeftPlots.inew.XData=[0];
            obj.hLeftPlots.inew.YData=[0];
            obj.hLeftPlots.inew.Color = [0 0 0];
            obj.hMiddlePlots.fnew.XData=[0];
            obj.hMiddlePlots.fnew.YData=[0];
            obj.hMiddlePlots.fnew.Color = [0 0 0];
            obj.hMiddlePlots.inew.XData=[0];
            obj.hMiddlePlots.inew.YData=[0];
            obj.hMiddlePlots.inew.Color = [0 0 0];
            % Turn buttons off
            obj.gTurnButtonsOff([1 1 1 1 1 1]);
            % Remove the values
            switch boundary
                case 1
                    obj.F2 = [];
                    obj.kF2 = [];
                    obj.hLeftPlots.dp1.XData=[0];
                    obj.hLeftPlots.dp1.YData=[0];
                    obj.hLeftPlots.dp1.Color = [1 1 1];
                    obj.hLeftPlots.f.XData=[0];
                    obj.hLeftPlots.f.YData=[0];
                    obj.hLeftPlots.f.Color = [1 1 1];
                    obj.I1 = [];
                    obj.kI1 = [];
                    obj.hMiddlePlots.dp4.XData=[0];
                    obj.hMiddlePlots.dp4.YData=[0];
                    obj.hMiddlePlots.dp4.Color = [1 1 1];
                    obj.hMiddlePlots.i.XData=[0];
                    obj.hMiddlePlots.i.YData=[0];
                    obj.hMiddlePlots.i.Color = [1 1 1];
                case 3
                    obj.F1 = [];
                    obj.kF1 = [];
                    obj.hMiddlePlots.dp3.XData=[0];
                    obj.hMiddlePlots.dp3.YData=[0];
                    obj.hMiddlePlots.dp3.Color = [1 1 1];
                    obj.hMiddlePlots.f.XData=[0];
                    obj.hMiddlePlots.f.YData=[0];
                    obj.hMiddlePlots.f.Color = [1 1 1];
                    obj.I2 = [];
                    obj.kI2 = [];
                    obj.hLeftPlots.dp2.XData=[0];
                    obj.hLeftPlots.dp2.YData=[0];
                    obj.hLeftPlots.dp2.Color = [1 1 1];
                    obj.hLeftPlots.i.XData=[0];
                    obj.hLeftPlots.i.YData=[0];
                    obj.hLeftPlots.i.Color = [1 1 1];
            end
        end
        %% gLoad: Load F/I
        function gLoad(obj, hObj, event,boundary)
            % Turn buttons off
            obj.gTurnButtonsOff([1 1 1 1 1 1]);
            % Load file (no safeguards for all error cases)
            [fname, fpath] = uigetfile({'*.mat'},...
                'Select *.mat file with only one matrix (point coordinates in rows)');
            if max(fname)>0 && max(fpath)>0
                S=load([fpath fname]);
                names=fieldnames(S);
                if not(isempty(names))
                    temp = getfield(S,names{1});
                else
                    disp('ERROR: empty mat file')
                end
            end
            % Stop if nothing is loaded
            if exist('temp','var')==0
                return
            end
            % Add the boundary
            switch boundary
                case 1
                    obj.F2 = temp;
                    obj.I1 = [];
                    obj.kI1 = [];
                    obj.hLeftPlots.dp1.XData=obj.F2(:,1);
                    obj.hLeftPlots.dp1.YData=obj.F2(:,2);
                    obj.hLeftPlots.dp1.Color = [1 0 0];
                    obj.hLeftPlots.dp1.LineStyle='none';
                    obj.hLeftPlots.dp1.Marker='*';
                    obj.hButtonsLeft.hButton1.Value = 1;
                    obj.gTurnButtonsOff([1 0 0 0 0 0])
                case 2
                    obj.I2 = temp;
                    obj.hLeftPlots.dp2.XData=obj.I2(:,1);
                    obj.hLeftPlots.dp2.YData=obj.I2(:,2);
                    obj.hLeftPlots.dp2.Color = [0 0 1];
                    obj.hLeftPlots.dp2.LineStyle='none';
                    obj.hLeftPlots.dp2.Marker='*';
                    obj.hButtonsLeft.hButton4.Value = 1;
                    obj.gTurnButtonsOff([0 1 0 0 0 0]);
                case 3
                    obj.F1 = temp;
                    obj.I2 = [];
                    obj.kI2 = [];
                    obj.hMiddlePlots.dp3.XData=obj.F1(:,1);
                    obj.hMiddlePlots.dp3.YData=obj.F1(:,2);
                    obj.hMiddlePlots.dp3.Color = [0 0 1];
                    obj.hMiddlePlots.dp3.LineStyle='none';
                    obj.hMiddlePlots.dp3.Marker='*';
                    obj.hButtonsMiddle.hButton1.Value = 1;
                    obj.gTurnButtonsOff([0 0 1 0 0 0]);
                case 4
                    obj.I1 = temp;
                    obj.hMiddlePlots.dp4.XData=obj.I1(:,1);
                    obj.hMiddlePlots.dp4.YData=obj.I1(:,2);
                    obj.hMiddlePlots.dp4.Color = [1 0 0];
                    obj.hMiddlePlots.dp4.LineStyle='none';
                    obj.hMiddlePlots.dp4.Marker='*';
                    obj.hButtonsMiddle.hButton4.Value = 1;
                    obj.gTurnButtonsOff([0 0 0 1 0 0]);
            end
        end
        %% gLoadD: Load data matrix
        function gLoadD(obj, hObj, event,boundary)
            % Turn buttons off
            obj.gTurnButtonsOff([1 1 1 1 1 1]);
            obj.gReset([],[],1);
            obj.gReset([],[],3);
            % Load file (no safeguards for all error cases)
            [fname, fpath] = uigetfile({'*.mat'},...
                'Select *.mat file with matrix D');
            if max(fname)>0 && max(fpath)>0
                S=load([fpath fname]);
                names=fieldnames(S);
                if not(isempty(names))
                    temp = getfield(S,'D');
                else
                    disp('ERROR: empty mat file')
                end
            end
            % Stop if nothing is loaded
            if exist('temp','var')==0
                return
            end
            obj.D = temp;
            obj.checkResults();
            surf(obj.hRightAxes,obj.D');
            obj.hRightAxes.ButtonDownFcn = {@(hObj,event) obj.gAxesCallback(hObj,event,3)};
            axtoolbar(obj.hRightAxes,'default');
            obj.gAcceptLoaded([],[]);
        end
        %% gAcceptLoaded: Reset after loading or results
        function gAcceptLoaded(obj, hObj, event)
            % Remove previous results            
            obj.hLeftPlots.fnew.XData=[0];
            obj.hLeftPlots.fnew.YData=[0];
            obj.hLeftPlots.fnew.Color = [0 0 0];
            obj.hLeftPlots.inew.XData=[0];
            obj.hLeftPlots.inew.YData=[0];
            obj.hLeftPlots.inew.Color = [0 0 0];
            obj.hMiddlePlots.fnew.XData=[0];
            obj.hMiddlePlots.fnew.YData=[0];
            obj.hMiddlePlots.fnew.Color = [0 0 0];
            obj.hMiddlePlots.inew.XData=[0];
            obj.hMiddlePlots.inew.YData=[0];
            obj.hMiddlePlots.inew.Color = [0 0 0];
            % Turn buttons off
            obj.gTurnButtonsOff([1 1 1 1 1 1]);
            xlim(obj.hLeftAxes);
            xlim(obj.hLeftAxes,'manual');
            ylim(obj.hLeftAxes,'manual');
            xlim(obj.hMiddleAxes);
            xlim(obj.hMiddleAxes,'manual');
            ylim(obj.hMiddleAxes,'manual');
            % I2 (and F1 with turnButtonsOff)
            obj.I2 = obj.Res.I2;
            obj.hLeftPlots.dp2.XData=obj.I2(:,1);
            obj.hLeftPlots.dp2.YData=obj.I2(:,2);
            obj.hLeftPlots.dp2.Color = [0 0 1];
            obj.hLeftPlots.dp2.LineStyle='none';
            obj.hLeftPlots.dp2.Marker='*';
            obj.hButtonsLeft.hButton4.Value = 1;
            obj.gTurnButtonsOff([0 1 0 0 0 0]);
            % I1 (and F2 with turnButtonsOff)
            obj.I1 = obj.Res.I1;
            obj.hMiddlePlots.dp4.XData=obj.I1(:,1);
            obj.hMiddlePlots.dp4.YData=obj.I1(:,2);
            obj.hMiddlePlots.dp4.Color = [1 0 0];
            obj.hMiddlePlots.dp4.LineStyle='none';
            obj.hMiddlePlots.dp4.Marker='*';
            obj.hButtonsMiddle.hButton4.Value = 1;
            obj.gTurnButtonsOff([0 0 0 1 0 0]);
        end
        %% gSave: Save results
        function gSave(obj, hObj, event)
            % Turn buttons off
            obj.gTurnButtonsOff([1 1 1 1 1 1]);
            % Prepare file names
            D = obj.D;
            F1 = obj.F1;
            I1 = obj.I1;
            F2 = obj.F2;
            I2 = obj.I2;
            % Inner and outer boundaries
            FIRPOL_1 = obj.F1(obj.kF1,:);
            INNPOL_1 = obj.I1(obj.kI1,:);
            FIRPOL_2 = obj.F2(obj.kF2,:);
            INNPOL_2 = obj.I2(obj.kI2,:);
            % Save
            uisave({'D','F1','I1','F2','I2','FIRPOL_1','INNPOL_1','FIRPOL_2','INNPOL_2'},'results.mat')
        end
        %% calcDual: calculate the dual plot from points
        function [pd,kd] = calcDual(obj,p,k,iold,ikold)
            sizeTemp = size(k,1)-1;
            pd = zeros(sizeTemp,2);
            pd(:,1) = (-p(k(1:sizeTemp),2)+p(k(2:sizeTemp+1),2))./...
                (p(k(1:sizeTemp),2).*p(k(2:sizeTemp+1),1)-p(k(2:sizeTemp+1),2).*p(k(1:sizeTemp),1));
            pd(:,2) = (-p(k(1:sizeTemp),1)+p(k(2:sizeTemp+1),1))./...
                (p(k(1:sizeTemp),1).*p(k(2:sizeTemp+1),2)-p(k(2:sizeTemp+1),1).*p(k(1:sizeTemp),2));
            kd = convhull(pd);
            % All of the old points inside of new INNPOL remain.
            if nargin==5
                % The old points from I, which are inside of the new conv
                % hull are appended at the end of the matrix. Note, that
                % this changes the order of the points in the trace.
                if ~isempty(iold)
                    [in,on] = inpolygon(iold(:,1),iold(:,2),pd(kd,1),pd(kd,2));
                    pd = [pd; iold(in&~on,1) iold(in&~on,2)];
                end
            end
        end
        %% calcMatrix: Here the optimization is started
        function err = calcMatrix(obj,startType,optimType,bndCoeff)
            % The options can be set as needed
            % e.g. 'Display','iter' or 'final' or 'none'
            % Starting vectors
            if startType==1
                startVectV = (ones(1,size(obj.I2,1)))./sqrt(size(obj.I2,1));
                startVectUS =  (ones(1,size(obj.I1,1)+3))./sqrt(size(obj.I1,1)+3);
            else
                startVectV = (rand(1,size(obj.I2,1)))./sqrt(size(obj.I2,1));
                startVectUS = (rand(1,size(obj.I1,1)+3))./sqrt(size(obj.I1,1)+3);
            end
            % Bounds for the minimal positive value in v1, u1
            bndV = bndCoeff/sqrt(size(obj.I2,1));
            bndUS = bndCoeff/sqrt(size(obj.I1,1)+3);
            % Various optimizations
            switch optimType
                case 1
                    options = optimoptions('fminunc','Display','final',...
                        'StepTolerance',1E-8,'OptimalityTolerance',1E-8,'MaxFunctionEvaluations',100000);
                    [v1,resnorm] = fminunc(@(v1) obj.recV(v1,obj.I2,bndV),...
                        startVectV,options);
                    [us1,resnorm2] = fminunc(@(us1) obj.recUS(us1,obj.I1,bndUS),...
                        startVectUS,options);
                case 2
                    options = optimset('Display','final',...
                        'TolX',1E-8,'TolFun',1E-8,'MaxIter',100000,'MaxFunEvals',100000);
                    [v1,resnorm] = fminsearch(@(v1) obj.recV(v1,obj.I2,bndV),...
                        startVectV,options);
                    [us1,resnorm2] = fminsearch(@(us1) obj.recUS(us1,obj.I1,bndUS),...
                        startVectUS,options);
                case 3
                    options = optimoptions('fsolve','Display','final',...
                        'MaxIterations',100000,'MaxFunctionEvaluations',1000000);
                    [v1,resnorm] = fsolve(@(v1) obj.recV(v1,obj.I2,bndV),...
                        startVectV,options);
                    [us1,resnorm2] = fsolve(@(us1) obj.recUS(us1,obj.I1,bndUS),...
                        startVectUS,options);
                case 4 % Not optimized: not vectorized, therefore similar to fsolve.
                    [v1,resnorm] = lsqnonlin(@(v1) obj.recV(v1,obj.I2,bndV),...
                        startVectV);
                    [us1,resnorm2] = lsqnonlin(@(us1) obj.recUS(us1,obj.I1,bndUS),...
                        startVectUS);
                case 5
                    [v1,resnorm] = patternsearch(@(v1) obj.recV(v1,obj.I2,bndV),...
                        startVectV);
                    [us1,resnorm2] = patternsearch(@(us1) obj.recUS(us1,obj.I1,bndUS),...
                        startVectUS);
                case 6
                    [v1,resnorm] = simulannealbnd(@(v1) obj.recV(v1,obj.I2,bndV),...
                        startVectV);
                    [us1,resnorm2] = simulannealbnd(@(us1) obj.recUS(us1,obj.I1,bndUS),...
                        startVectUS);
                case 7
                    [v1,resnorm] = fminimax(@(v1) obj.recV(v1,obj.I2,bndV),...
                        startVectV);
                    [us1,resnorm2] = fminimax(@(us1) obj.recUS(us1,obj.I1,bndUS),...
                        startVectUS);
                case 8
                    [v1,resnorm] = particleswarm(@(v1) obj.recV(v1,obj.I2,bndV),...
                        size(obj.I2,1));
                    [us1,resnorm2] = particleswarm(@(us1) obj.recUS(us1,obj.I1,bndUS),...
                        size(obj.I1,1)+3);
                case 9
                    [v1,resnorm] = ga(@(v1) obj.recV(v1,obj.I2,bndV),...
                        size(obj.I2,1));
                    [us1,resnorm2] = ga(@(us1) obj.recUS(us1,obj.I1,bndUS),...
                        size(obj.I1,1)+3);
            end
            % Calculate the SVD
            V(:,1) = v1';
            V(:,2) = V(:,1).*obj.I2(:,1);
            V(:,3) = V(:,1).*obj.I2(:,2);
            Sig = abs(diag(us1(1:3)));
            U(:,1) = (us1(4:end))';
            U(:,2) = Sig(1,1)/Sig(2,2)*U(:,1).*obj.I1(:,1);
            U(:,3) = Sig(1,1)/Sig(3,3)*U(:,1).*obj.I1(:,2);
            % Calculate the matrix and the total residual
            obj.D = U*Sig*V';
            err = resnorm+resnorm2;
        end
        %% checkResults: Calculate INNPOL and FIRPOL for both factors for a given D
        function checkResults(obj)
            if isempty(obj.F1)
                % Recalculate the SVD with the standard orientation
                [U,Sig,V]=svds(obj.D,3);
                p = zeros(3);
                for i = 1:3
                    if max(V(:,i))>=max(-V(:,i))
                        p(i,i) = 1;
                    else
                        p(i,i) = -1;
                    end
                end
                U = U(:,1:3)*p;
                Sig = Sig(1:3,1:3);
                V = V(:,1:3)*p;
                % Calculate the boundaries
                obj.calcBounds(U,Sig,V)
            else
                % Stop axes from flipping
                % Find the correct orientation for D, depending on the
                % previous results
                mindf = Inf;
                for j = [1,-1]
                    for k = [1,-1]
                        [U,Sig,V]=svds(obj.D,3);
                        p = zeros(3);
                        if max(V(:,1))>=max(-V(:,1))
                            p(1,1) = 1;
                        else
                            p(1,1) = -1;
                        end
                        p(2,2) = j;
                        p(3,3) = k;
                        U = U(:,1:3)*p;
                        Sig = Sig(1:3,1:3);
                        V = V(:,1:3)*p;
                        % Calculate the boundaries
                        obj.calcBounds(U,Sig,V)
                        df = norm([obj.I1;obj.I2]-[obj.Res.I1;obj.Res.I2]);
                        if df<mindf
                            mindf = df;
                            minj = j;
                            mink = k;
                        end
                    end
                end
                [U,Sig,V]=svds(obj.D,3);
                p = zeros(3);
                if max(V(:,1))>=max(-V(:,1))
                    p(1,1) = 1;
                else
                    p(1,1) = -1;
                end
                p(2,2) = minj;
                p(3,3) = mink;
                U = U(:,1:3)*p;
                Sig = Sig(1:3,1:3);
                V = V(:,1:3)*p;
                % Calculate the boundaries
                obj.calcBounds(U,Sig,V)
            end
        end
        %% calcBounds: Calculate the inner and the outer boundary
        function calcBounds(obj,U,Sig,V)
            % Calculate the boundaries
            sizeD = [size(U,1) size(V,1)];
            US = U*Sig;
            % Primal INNPOL
            O1 = zeros(sizeD(1));
            for i=1:sizeD(1)
                O1(i,i) = 1/US(i,1);
            end
            obj.Res.I1 = O1*obj.D*V(:,2:3);
            obj.Res.kI1 = convhull(obj.Res.I1);
            % Dual FIRPOL
            [obj.Res.F2,obj.Res.kF2] = obj.calcDual(obj.Res.I1,obj.Res.kI1);
            % Dual INNPOL
            O2 = zeros(sizeD(2));
            for i=1:sizeD(2)
                O2(i,i) = 1/V(i,1);
            end
            obj.Res.I2 = O2*V(:,2:3);
            obj.Res.kI2 = convhull(obj.Res.I2);
            % Primal FIRPOL
            [obj.Res.F1,obj.Res.kF1] = obj.calcDual(obj.Res.I2,obj.Res.kI2);
        end
    end
    %% Optimization functions
    methods(Static)
        function c = recV(v1,dpV,bnd)
            % Defining matrices according to the theory
            v1 = v1';
            n1 = norm(v1)-1;
            n2 = norm(v1.*dpV(:,1))-1;
            n3 = norm(v1.*dpV(:,2))-1;
            o12 = v1'*(v1.*dpV(:,1));
            o13 = v1'*(v1.*dpV(:,2));
            o23 = (v1.*dpV(:,1))'*(v1.*dpV(:,2));
            nn = abs(min(0,v1-bnd));
            c = norm([n1;n2;n3;o12;o13;o23;nn]);
        end
        function c = recUS(us1IN,dpUS,bnd)
            us1IN = us1IN';
            Sig = us1IN(1:3);
            u1 = us1IN(4:end);
            % Defining matrices according to the theory
            n1 = norm(u1)-1;
            n2 = Sig(1)^2*(u1.*dpUS(:,1))'*(u1.*dpUS(:,1))-Sig(2)^2;
            n3 = Sig(1)^2*(u1.*dpUS(:,2))'*(u1.*dpUS(:,2))-Sig(3)^2;
            o12 = u1'*(u1.*dpUS(:,1));
            o13 = u1'*(u1.*dpUS(:,2));
            o23 = (u1.*dpUS(:,1))'*(u1.*dpUS(:,2));
            nn = abs(min(0,u1-bnd));
            s12 = min(0,Sig(1)^2-Sig(2)^2);
            s23 = min(0,Sig(2)^2-Sig(3)^2);
            c = norm([n1;n2;n3;o12;o13;o23;nn;s12;s23]);
        end
    end
end