function [FPU, FPV, essU, essV] = fpess4noise(V, US, nc, epsRand, noise, nr_ray)
% Computation of the outer polygon and the associated essential frequencies.
% Copyright (c) 2020 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.
% 
% 
% input:
% V ......... (ncxn)-matrix of right singular vectors, with
%                   nc being the number of singular vectors (min. 2),
%                   n being the number of channels (resp. spectra)
% US ......... (kxnc)-matrix of left singular vectors times singular values, with
%                   nc being the number of singular vectors (min. 2),
%                   k being the number of channels (resp. spectra)
% nc ........ number of components
% epsRand ... boundary precision (e.g. 1.E-4)
% noise ..... permitted negativity (<0)
% nr_ray .... number of rays that are used in the ray cast algorithm (only 
%             necessary for s>2
%             for s=3: min. 10, max. 1000
%             for s=4: min. 200, max. 502001 (max. 500 rays in each direction)
%             for s=5: min. 2000, max. 251503001 (max. 500 rays in each direction)
%             for s=6: min. 20000, max. 1.2600e+11 (max. 500 rays in each direction)
% 
% output:
% FPV ........ boundary points of firpol/the outer polygon in V-space
% essV ....... esential channels in V-space
% FPU ........ boundary points of firpol/the outer polygon in U-space
% essU ....... esential channels in U-space

if nc==2
    [FPV, essV] = FPcomp_2comp(V, epsRand, noise);
    [FPU, essU] = FPcomp_2comp(US, epsRand, noise);
elseif nc==3
    [FPV, essV] = FPcomp_3comp(V, epsRand, noise, nr_ray);
    [FPU, essU] = FPcomp_3comp(US, epsRand, noise, nr_ray);
elseif nc==4
    [FPV, essV] = FPcomp_4comp(V, epsRand, noise, nr_ray);
    [FPU, essU] = FPcomp_4comp(US, epsRand, noise, nr_ray);
elseif nc==5
    [FPV, essV] = FPcomp_5comp(V, epsRand, noise, nr_ray);
    [FPU, essU] = FPcomp_5comp(US, epsRand, noise, nr_ray);
elseif nc==6
    [FPV, essV] = FPcomp_6comp(V, epsRand, noise, nr_ray);
    [FPU, essU] = FPcomp_6comp(US, epsRand, noise, nr_ray);
else
    FPV = [];
    essV = [];
    FPU = [];
    essU = [];
end
end


function [FP, rel] = FPcomp_2comp(V, epsRand, noise)
% Computation of the outer polygon and the associated essential frequencies
% for 2 component systems.
% 
% input:
% V ......... (sxn)-matrix singular vectors, with
%                   s is the number of singular vectors (min. 2),
%                   n is the number of channels (resp. spectra)
% epsRand ... boundary precision
% noise ..... permitted negativity (<0)
% 
% output:
% FP ........ boundary points of firpol/the outer polygon
% rel ....... esential channels (resp. spectra)


%%
[s,n] = size(V); 

if n<s
   V=V';
end
[s,n] = size(V); 

if s<2 || n<2 
    error('Not enough singular vectors')
end

v = [-1 1]; % using 2 Rays, because it's 1D
rel = [];
FP = [];


for vv=v
    % get current ray
    v1=vv;
    outv = 0;
    
    % get a point on the ray so that it is outside firpol
    idok = 1;
    while idok
        mi = V(1,1)+V(2,1)*v1;
        ma = mi;
        for i=2:n
            mi = min(mi, V(1,i)+V(2,i)*v1);
            ma = max(ma, V(1,i)+V(2,i)*v1);
        end        
        if mi/ma > noise - epsRand
            outv = v1; % point that is still in firpol
            v1 = 2*v1;
        else
            idok=0;
        end
    end
    
    outa = v1; % point that is outside firpol
    h= v1/2;   % distance between outa and outv
    
    % bisection algorithm to find the boundary of firpol
    while h > epsRand 
        v1 = (outa + outv)/2;
        mi = V(1,1)+V(2,1)*v1;
        ma = mi;
        for i=2:n
            mi = min(mi, V(1,i)+V(2,i)*v1);
            ma = max(ma, V(1,i)+V(2,i)*v1);
        end
        if mi/ma > noise - epsRand
            outv = v1;
        else
            outa = v1;
        end
        h=h/2;
    end
    FP=[FP; outv]; % points of the boundary of firpol
    
    v1=outv;
    mi = V(1,1)+V(2,1)*v1;
    ma = mi;
    for i=1:n
        ma = max(ma, V(1,i)+V(2,i)*v1);
        mi = min(mi, V(1,i)+V(2,i)*v1);
    end
    delta = abs(1.01 * (noise * ma - mi)/ma); % control parameter
    
    % find the essential spectra
    for j=1:n
        mi =  V(1,j)+V(2,j)*v1;
        if mi/ma <= noise + delta 
            rel=[rel,j];
        end
    end    
end
rel=unique(rel);
end

function [FP, rel] = FPcomp_3comp(V, epsRand, noise, nr_ray)
% Computation of the outer polygon and the associated essential frequencies
% for 3 component systems.
% 
% input:
% V ......... (sxn)-matrix singular vectors, with
%                   s is the number of singular vectors (min. 2),
%                   n is the number of channels (resp. spectra)
% epsRand ... boundary precision
% noise ..... permitted negativity (<0)
% nr_ray .... number of rays that are used in the ray cast algorithm
%             (min. 10, max. 1000) 
% 
% output:
% FP ........ boundary points of firpol/the outer polygon
% rel ....... esential channels (resp. spectra)


%%
[s,n] = size(V);

if n<s
   V=V';
end
[s,n] = size(V); 

if s<3 || n<3 
    error('Not enough singular vectors')
end

rel = [];
FP = [];

m1 = max(10, min(2000, nr_ray));

for i1=1:m1
    % get current ray
    p1 = 2*pi*i1/m1 -pi;
    
    v1(1)=0.5*cos(p1);
    v1(2)=0.5*sin(p1);
    outv = zeros(1,2);
    
    % get a point on the ray so that it is outside firpol
    idok = 1;
    while idok
        mi = V(1,1)+V(2,1)*v1(1)+V(3,1)*v1(2);
        ma = mi;
        for i=2:n
            mi = min(mi, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2));
            ma = max(ma, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2));
        end
        
        if mi/ma > noise - epsRand
            outv = v1; % point that is still in firpol
            v1 = 2*v1;
        else
            idok=0;
        end
        
    end
    outa = v1; % point that is outside firpol
    h = norm(outa-outv);
    
    % bisection algorithm to find the boundary of firpol
    while h > epsRand
        v1 = (outa + outv)/2;
        
        mi = V(1,1)+V(2,1)*v1(1)+V(3,1)*v1(2);
        ma = mi;
        for i=2:n
            mi = min(mi, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2));
            ma = max(ma, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2));
        end
        
        if mi/ma > noise - epsRand
            outv = v1;
        else
            outa = v1;
        end
        h=h/2;
    end
    FP=[FP; outv]; % points of the boundary of firpol
    
    v1=outv;
    mi = V(1,1)+V(2,1)*v1(1)+V(3,1)*v1(2);
    ma = mi;
    
    for i=1:n
        mi = min(mi, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2));
        ma = max(ma, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2));
    end
    
    delta = abs(1.01 * (noise * ma - mi)/ma); % control parameter
    
    % find the essential spectra
    for j=1:n
        mi =  V(1,j)+V(2,j)*v1(1)+V(3,j)*v1(2);
        if mi/ma <= noise + delta 
            rel=[rel,j];
        end
    end       
end
rel=unique(rel);
end

function [FP, rel] = FPcomp_4comp(V, epsRand, noise, nr_ray)
% Computation of the outer polygon and the associated essential frequencies
% for 4 component systems.
% 
% input:
% V ......... (sxn)-matrix singular vectors, with
%                   s is the number of singular vectors (min. 2),
%                   n is the number of channels (r% control parameteresp. spectra)
% epsRand ... boundary precision
% noise ..... permitted negativity (<0)
% nr_ray .... number of rays that are used in the ray cast algorithm
%             (min. 200, max. 502001 (max. 500 rays in each direction)) 
% 
% output:
% FP ........ boundary points of firpol/the outer polygon
% rel ....... esential channels (resp. spectra)


%% 
[s,n] = size(V); 

if n<s
   V=V';
end
[s,n] = size(V); 

if s<4 || n<4
    error('Not enough singular vectors')
end

rel = [];
FP = [];

m2 = max(10, min(500, floor(sqrt(nr_ray/2))));
m1 = 2*m2;

for i1=1:m1
    for i2=1:m2
        % get current ray
        p1 = pi*i2/m2;
        p2 = 2*pi*i1/m1 -pi;
        
        v1(1)=0.5*cos(p1);
        v1(2)=0.5*sin(p1)*cos(p2);
        v1(3)=0.5*sin(p1)*sin(p2);
        
        outv = zeros(1,3);
        
        % get a point on the ray so that it is outside firpol
        idok = 1;
        while idok
            mi = V(1,1)+V(2,1)*v1(1)+V(3,1)*v1(2)+V(4,1)*v1(3);
            ma = mi;
            for i=2:n
                mi = min(mi, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3));
                ma = max(ma, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3));
            end
            if mi/ma > noise - epsRand
                outv = v1; % point that is still in firpol
                v1 = 2*v1;
            else
                idok=0;
            end
        end
        outa = v1; % point that is outside firpol
        h=norm(outa-outv);
        
        % bisection algorithm to find the boundary of firpol
        while h > epsRand
            v1 = (outa + outv)/2;
            
            mi = V(1,1)+V(2,1)*v1(1)+V(3,1)*v1(2)+V(4,1)*v1(3);
            ma = mi;
            for i=2:n
                mi = min(mi, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3));
                ma = max(ma, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3));
            end
            
            if mi/ma > noise - epsRand
                outv = v1;
            else
                outa = v1;
            end
            h=h/2;
        end
        FP=[FP; outv]; % points of the boundary of firpol
        
        v1=outv;
        mi = V(1,1)+V(2,1)*v1(1)+V(3,1)*v1(2)+V(4,1)*v1(3);
        ma = mi;
        
        for i=1:n
            ma = max(ma, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3));
            mi = min(mi, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3));      
        end
        
        delta = abs(1.01 * (noise * ma - mi)/ma); % control parameter
    
        % find the essential spectra
        for j=1:n
            mi =  V(1,j)+V(2,j)*v1(1)+V(3,j)*v1(2)+V(4,j)*v1(3);
            if mi/ma <= noise + delta 
                rel=[rel,j];
            end
        end
    end
end
rel=unique(rel);
end

function [FP, rel] = FPcomp_5comp(V, epsRand, noise, nr_ray)
% Computation of the outer polygon and the associated essential frequencies
% for 5 component systems.
% 
% input:
% V ......... (sxn)-matrix singular vectors, with
%                   s is the number of singular vectors (min. 2),
%                   n is the number of channels (resp. spectra)
% epsRand ... boundary precision
% noise ..... permitted negativity (<0)
% nr_ray .... number of rays that are used in the ray cast algorithm
%             (min. 2000, max. 251503001 (max. 500 rays in each direction))  
% 
% output:
% FP ........ boundary points of firpol/the outer polygon
% rel ....... esential channels (resp. spectra)


%% 
[s,n] = size(V); 

if n<s
   V=V';
end
[s,n] = size(V); 

if s<5 || n<5 
    error('Not enough singular vectors')
end

rel = [];
FP = [];

m2 = max(10, min(500, floor((nr_ray/2)^(1/3))));
m1 = 2*m2;
m3 = m2;

for i1=1:m1
    for i2=1:m2
        for i3=1:m3
            % get current ray
            p1 = pi*i3/m3;
            p2 = pi*i2/m2;
            p3 = 2*pi*i1/m1 -pi;
                        
            v1(1)=0.5*cos(p1);
            v1(2)=0.5*sin(p1)*cos(p2);
            v1(3)=0.5*sin(p1)*sin(p2)*cos(p3);
            v1(4)=0.5*sin(p1)*sin(p2)*sin(p3);
            
            outv = zeros(1,4);
            
            % get a point on the ray so that it is outside firpol
            idok = 1;
            while idok
                mi = V(1,1)+V(2,1)*v1(1)+V(3,1)*v1(2)+V(4,1)*v1(3)+V(5,1)*v1(4);
                ma = mi;
                for i=2:n
                    mi = min(mi, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3)+V(5,i)*v1(4));
                    ma = max(ma, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3)+V(5,i)*v1(4));
                end
                if mi/ma > noise - epsRand
                    outv = v1;  % point that is still in firpol
                    v1 = 2*v1;
                else
                    idok=0;
                end
            end
            outa = v1; % point that is outside firpol            
            h=norm(outa-outv);
            
            % bisection algorithm to find the boundary of firpol
            while h > epsRand
                v1 = (outa + outv)/2;
                
                mi = V(1,1)+V(2,1)*v1(1)+V(3,1)*v1(2)+V(4,1)*v1(3)+V(5,1)*v1(4);
                ma = mi;
                for i=2:n
                    mi = min(mi, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3)+V(5,i)*v1(4));
                    ma = max(ma, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3)+V(5,i)*v1(4));
                end
                
                if mi/ma > noise - epsRand
                    outv = v1;
                else
                    outa = v1;
                end                
                h=h/2;                
            end
            FP=[FP; outv]; % points of the boundary of firpol
            
            v1=outv;
            mi = V(1,1)+V(2,1)*v1(1)+V(3,1)*v1(2)+V(4,1)*v1(3)+V(5,1)*v1(4);
            ma = mi;
            
            for i=1:n
                ma = max(ma, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3)+V(5,i)*v1(4));
                mi = min(mi, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3)+V(5,i)*v1(4));
            end
            
            delta = abs(1.01 * (noise * ma - mi)/ma); % control parameter
            
            % find the essential spectra
            for j=1:n
                mi =  V(1,j)+V(2,j)*v1(1)+V(3,j)*v1(2)+V(4,j)*v1(3)+V(5,j)*v1(4);
                if mi/ma <= noise + delta
                    rel=[rel,j];
                end
            end            
        end
    end
end
rel=unique(rel);
end

function [FP, rel] = FPcomp_6comp(V, epsRand, noise, nr_ray)
% Computation of the outer polygon and the associated essential spectral
% for 6 components.
% 
% input:
% V ......... (sxn)-matrix singular vectors, with
%                   s is the number of singular vectors (min. 2),
%                   n is the number of channels (resp. spectra)
% epsRand ... boundary precision
% noise ..... permitted negativity (<0)
% nr_ray .... number of rays that are used in the ray cast algorithm
%             (min. 20000, max. 1.2600e+11 (max. 500 rays in each direction)) 
% 
% output:
% FP ........ boundary points of firpol/the outer polygon
% rel ....... esential channels (resp. spectra)


%%
[s,n] = size(V); 

if n<s
   V=V';
end
[s,n] = size(V); 

if s<6 || n<6 
    error('Not enough singular vectors')
end

rel = [];
FP = [];

m2 = max(10, min(500, floor((nr_ray/2)^(1/4))));
m1 = 2*m2;
m3 = m2;
m4 = m3;

for i1=1:m1
    for i2=1:m2
        for i3=1:m3
            for i4=1:m4
                % get current ray
                p1 = pi*i4/m4;
                p2 = pi*i3/m3;
                p3 = pi*i2/m2;
                p4 = 2*pi*i1/m1 -pi;
                
                v1(1)=0.5*cos(p1);
                v1(2)=0.5*sin(p1)*cos(p2);
                v1(3)=0.5*sin(p1)*sin(p2)*cos(p3);
                v1(4)=0.5*sin(p1)*sin(p2)*sin(p3)*cos(p4);
                v1(5)=0.5*sin(p1)*sin(p2)*sin(p3)*sin(p4);
                
                outv = zeros(1,5);
                
                % get a point on the ray so that it is outside firpol
                idok = 1;
                while idok
                    mi = V(1,1)+V(2,1)*v1(1)+V(3,1)*v1(2)+V(4,1)*v1(3)+V(5,1)*v1(4)+V(6,1)*v1(5);
                    ma = mi;
                    for i=2:n
                        mi = min(mi, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3)+V(5,i)*v1(4)+V(6,i)*v1(5));
                        ma = max(ma, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3)+V(5,i)*v1(4)+V(6,i)*v1(5));
                    end
                    
                    if mi/ma > noise - epsRand
                        outv = v1; % point that is still in firpol
                        v1 = 2*v1;
                    else
                        idok=0;
                    end
                    
                end
                outa = v1; % point that is outside firpol                
                h=norm(outa-outv);
                
                % bisection algorithm to find the boundary of firpol
                while h > epsRand
                    v1 = (outa + outv)/2;
                    
                    mi = V(1,1)+V(2,1)*v1(1)+V(3,1)*v1(2)+V(4,1)*v1(3)+V(5,1)*v1(4)+V(6,1)*v1(5);
                    ma = mi;
                    for i=2:n
                        mi = min(mi, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3)+V(5,i)*v1(4)+V(6,i)*v1(5));
                        ma = max(ma, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3)+V(5,i)*v1(4)+V(6,i)*v1(5));
                    end
                    
                    if mi/ma > noise - epsRand
                        outv = v1;
                    else
                        outa = v1;
                    end                    
                    h=h/2;                    
                end
                FP=[FP; outv]; % points of the boundary of firpol
                
                v1=outv;
                mi = V(1,1)+V(2,1)*v1(1)+V(3,1)*v1(2)+V(4,1)*v1(3)+V(5,1)*v1(4)+V(6,1)*v1(5);
                ma = mi;
                
                for i=1:n
                    ma = max(ma, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3)+V(5,i)*v1(4)+V(6,i)*v1(5));
                    mi = min(mi, V(1,i)+V(2,i)*v1(1)+V(3,i)*v1(2)+V(4,i)*v1(3)+V(5,i)*v1(4)+V(6,i)*v1(5));
                end
                
                delta = abs(1.01 * (noise * ma - mi)/ma); % control parameter
                
                % find the essential spectra
                for j=1:n
                    mi =  V(1,j)+V(2,j)*v1(1)+V(3,j)*v1(2)+V(4,j)*v1(3)+V(5,j)*v1(4)+V(6,j)*v1(5);
                    if mi/ma <= noise + delta
                        rel=[rel,j];
                    end
                end                
            end
        end
    end
end
rel=unique(rel);
end

