function arrayout = aggarr(arrayin,M,N,interval)
%
% arrayout = aggarr(arrayin,M,N,interval)
% 
% PURPOSE:
%   This function aggregates defined subsets of the m-by-n matrix
%   "arrayin".
% 
%   If the "interval" parameter is set to zero, the function aggregates
%   each M-by-N subset of "arrayin" directly,so that the resulting matrix
%   has dimensions p-by-q, where p=m/M and q=n/N.
%
%   If "interval" is set to one, it aggregates subsets composed of every
%   Mth element on the rows and every Nth element on the columns, so that
%   the resulting matrix has dimensions M-by-N.
% 
% INPUTS:
%   arrayin (m x n) - The matrix to be aggregated
%   M               - no. of rows in the subset to be aggregated
%   N               - no. of cols in the subset to be aggregated
%   interval        - 0 or 1. Switch between two aggregation modes
% 
% OUTPUTS:
%   arrayout        - The aggregated matrix
% 
% COMMENTS:
% Explanation in terms of MRIO:
% Interval = 0 is equal to aggregating all sectors in each country to get a
% table of aggregate flows between countries
% Interval = 1 is equal to aggregating all countries for each sector to get
% a table of aggregate global flows between each economic sector
% 
% Both cases require M and N to satisfy the condition
% M*p=m, N*q=n, where p,q are non-zero natural numbers.
% 
% Numerical example:
% A = [1  1  1;
%      2  2  2;
%      3  3  3;
%      4  4  4];
%
% aggarr(A,2,3,0) = [9;
%                   21]
% aggarr(A,2,3,1) = [4 4 4;
%                    6 6 6]
%
% Kjartan Steen-Olsen, NTNU, 2011

[m,n] = size(arrayin);
p=m/M;
q=n/N;

if isnan((p-floor(p))/(q-floor(q)))
    
    if interval == 0
        arrayout = zeros(p,q);
        for i=1:size(arrayout,1)
            for j=1:size(arrayout,2)
                arraysubset = arrayin(M*(i-1)+1:M*i,N*(j-1)+1:N*j);
                arrayout(i,j) = sum(sum(arraysubset));
            end
        end
        
    elseif interval==1
        arrayout = zeros(M,N);
        for i=1:size(arrayout,1)
            for j=1:size(arrayout,2)
                arraysubset = arrayin(i:M:end,j:N:end);
                arrayout(i,j) = sum(sum(arraysubset));
            end
        end
        
    else
        disp('Invalid "interval" parameter".')
    end
    
else
    disp('Invalid dimension parameter (M or N).')
end
