%% The structural path analysis (SPA) algorithm to extract critical individual supply chain paths.
% This is a sample code used in our JIE paper:
% Sai Liang, Sen Guo, Joshua P. Newell, Shen Qu, Yu Feng, Anthony S.F.
% Chiu, and Ming Xu. Global Drivers of Russian Timber Harvest. Journal of
% Industrial Ecology. 2016, DOI:10.1111/jiec.12417.

%% INPUT variables
% A: Direct input coefficient matrix
% R: Sectoral intensity vector
% X: Final demand vector
% ind: investigated sector (starting sector)
% coverage: percent of total coverage to be analyzed
% SelectedNO: Number of top supply chain paths to keep
% step: production layers to be searched 

%% OUTPUT variables
% SelectedList: results indicating the contribution of each top supply chain path and the sectors it passed

%% codes
function SelectedList = SPA(A,R,X,ind,coverage,SelectedNO,step)

    n=length(A); % n is the number of sector
    y=X;

    adjMax = max(max(A));
    if adjMax < 1
        total=R*inv(eye(n)-A)*y;   

        step=0;
        cov=R*y;
        while cov<=(coverage*total)
            step=step+1;
            cov=cov+R*A^step*y;
        end       
    end   
    
    SelectedList = zeros(SelectedNO,step+1);

    if step<=2 && step~=0
        
        
        for i1=n:-1:1

            for i2=n:-1:0

                if i2 == 0

                    temp = R(ind)*A(ind,i1);
                    SelectedList = Compare_2(SelectedList,temp,i1,i2);   
                    break;
                end

                temp = R(ind)*A(ind,i1)*A(i1,i2);                
                SelectedList = Compare_2(SelectedList,temp,i1,i2);  

            end
        end
    end

    if step == 3 
        SelectedList = zeros(SelectedNO,4);

        for i1=n:-1:1

            for i2=n:-1:0

                if i2 == 0

                    temp = R(ind)*A(ind,i1);
                    SelectedList = Compare_3(SelectedList,temp,i1,i2,0);   
                    break;
                end
                for i3=n:-1:0      
                    if i3 == 0
                        temp = R(ind)*A(ind,i1)*A(i1,i2);
                        SelectedList = Compare_3(SelectedList,temp,i1,i2,i3);   
                        break;
                    end

                    temp = R(ind)*A(ind,i1)*A(i1,i2)*A(i2,i3);                
                    SelectedList = Compare_3(SelectedList,temp,i1,i2,i3);  

                end
            end
        end
    end
    
    if step >=4

        for i1=n:-1:1

            for i2=n:-1:0

                if i2 == 0

                    temp = R(ind)*A(ind,i1);
                    SelectedList = Compare_3(SelectedList,temp,i1,i2,0);   
                    break;
                end
                for i3=n:-1:0      
                    if i3 == 0
                        temp = R(ind)*A(ind,i1)*A(i1,i2);
                        SelectedList = Compare_3(SelectedList,temp,i1,i2,i3);   
                        break;
                    end
                    
                    for i4=n:-1:0    

                        if i4 == 0
                            temp = R(ind)*A(ind,i1)*A(i1,i2)*A(i2,i3);
                            SelectedList = Compare_4(SelectedList,temp,i1,i2,i3,i4);   
                            break;
                            
                        end

                        temp = R(ind)*A(ind,i1)*A(i1,i2)*A(i2,i3)*A(i3,i4);                
                        SelectedList = Compare_4(SelectedList,temp,i1,i2,i3,i4);  
                    end

                end
            end
        end
        
    end

    SelectedList = sortrows(SelectedList,-1);
    
    [Srow, ~] = size(SelectedList);
    for i = 1:Srow
        temp = 1;
        pathStep = sum(SelectedList(i,:)>0) - 1;
        for j = 1+1:pathStep+1
                       
            if j == 2 && j < pathStep+1
                temp = R(ind)*A(ind,SelectedList(i,j));
            elseif j == 2 && j == pathStep+1
                temp = R(ind)*A(ind,SelectedList(i,j))*X(SelectedList(i,j));
            elseif j<pathStep+1
                temp = temp*A(SelectedList(i,j-1),SelectedList(i,j));
            elseif j == pathStep+1 && j~=2
                temp = temp*A(SelectedList(i,j-1),SelectedList(i,j))*X(SelectedList(i,j));
            end
        end 
        
        SelectedList(i,1) = temp;
    end
    SelectedList = sortrows(SelectedList,-1);
end

function SelectedList = Compare_4(SelectedList,temp,i1,i2,i3,i4)
    if temp > SelectedList(1,1)
        SelectedList(1,1) = temp;
        SelectedList(1,2) = i1;
        SelectedList(1,3) = i2;
        SelectedList(1,4) = i3;
        SelectedList(1,5) = i4;
    end

    SelectedList = sortrows(SelectedList,1);
end

function SelectedList = Compare_3(SelectedList,temp,i1,i2,i3)
    if temp > SelectedList(1,1)
        SelectedList(1,1) = temp;
        SelectedList(1,2) = i1;
        SelectedList(1,3) = i2;
        SelectedList(1,4) = i3;
    end

    SelectedList = sortrows(SelectedList,1);
end

function SelectedList = Compare_2(SelectedList,temp,i1,i2)
    if temp > SelectedList(1,1)
        SelectedList(1,1) = temp;
        SelectedList(1,2) = i1;
        SelectedList(1,3) = i2;
    end

    SelectedList = sortrows(SelectedList,1);
end
