%% ExpectationMetrics
% This function embeds two pitch collections into respective expectation
% tensors and determines their distance using a standard metric.

% Authors: Andrew J. Milne and William A. Sethares
% Revision: 1.00  Date: 2011/06/24

% Reference: Milne, A.J., Sethares, W.A., Laney, R., Sharp, D.B. (2011)
% Modelling the Similarity of Pitch Collections with Expectation Tensors,
% Journal of Mathematics and Music, 5(1), 1-20.

% period = size, in cents, of the octave (or other period)
% q = number of elements for each dimension of the expectation array
% metric = metric used
% p = norm used for l_p metric 
% sigma = standard deviation of the discrete Gaussian
% win_len = the length of the smoothing kernel (in standard deviations)
% x_pc = pitch class collection x
% x_w = saliences of x
% y_pc = pitch class collection y
% y_w = saliences of y
% abs_or_rel = 'abs' or 'rel' - i.e., absolute or relative pitch embedding
% r = 1, 2, 3, or 4 - i.e., monad, dyad, triad, or tetrad embedding

% Note that the procedures run at O(d(q^r)), where d is the number of
% elements in x_pc (or y_pc). For higher values of r or d, it may be
% necessary to reduce q. High values of r and q may also lead to "out of
% memory" errors, due to the huge size of the resulting arrays.

function d = ExpectationMetrics(period, q, metric, p, sigma, win_len, x_pc, x_w, y_pc, y_w, abs_or_rel, r, exp_or_sal)

global shift_by virtual

%% Generating and transforming variables
x_pc = mod(round(x_pc*q/period),q);
if x_w == 0
    x_w = ones(1, length(x_pc));
end
d_x = length(x_pc);
if length(x_pc) ~= length(x_w)
    error('x_pc and x_w must have the same number of elements')
end

y_pc = mod(round(y_pc*q/period),q);
if y_w == 0
    y_w = ones(1, length(y_pc));
end
d_y = length(y_pc);
if length(y_pc) ~= length(y_w)
    error('y_pc and y_w must have the same number of elements')
end

windowlength = ceil(win_len*sigma*q/period);
n = -(windowlength-1)/2 : (windowlength-1)/2;
gwin = exp(-(n.^2) / (2 * (sigma*q/period)^2));
gwin = gwin/norm(gwin,1);

if isequal(metric, 'l_p') || isequal(metric, 'cosine')
else
    error('unknown metric - please enter "l_p" or "cosine"')
end

q_squared = q^2;
q_cubed = q^3;
q_fourth = q^4;

%% Pitch class response matrices X and Y
X_pcr_ij = zeros(d_x,q);
for i = 1:d_x
    X_pcr_ij(i,x_pc(i)+1) = x_w(i);
    if length(gwin) > 1
        X_pcr_ij(i,:) = real(ifft(fft(X_pcr_ij(i,:), q).*conj(fft(gwin, q)), q));
    end
end
X_pcr_ij(X_pcr_ij < 0.0000001) = 0;

Y_pcr_ij = zeros(d_y,q);
for i = 1:d_y
    Y_pcr_ij(i,y_pc(i)+1) = y_w(i);
    if length(gwin) > 1
        Y_pcr_ij(i,:) = real(ifft(fft(Y_pcr_ij(i,:), q).*conj(fft(gwin, q)), q));
    end
end
Y_pcr_ij(Y_pcr_ij < 0.0000001) = 0;

if isequal(abs_or_rel,'abs') && r==1
%% AbsMonad 
% This function embeds two pitch collections into absolute expectation
% tensors and then determines their distance using a metric.

% Expectation array X
 x_e_j = circshift(sum(X_pcr_ij),[0 floor(windowlength/2)-0]);
%x_e_j = sum(X_pcr_ij);
% figure(98)
% plot(x_e_j)

% Expectation array Y
 y_e_j = circshift(sum(Y_pcr_ij),[0 floor(windowlength/2)-0]);
%y_e_j = sum(Y_pcr_ij);
% figure(99)
% plot(y_e_j)

if isequal(exp_or_sal,'sal')
    % Salience array X
    x_e_j = 1-prod(1-X_pcr_ij);
    % Salience array Y
    y_e_j = 1-prod(1-Y_pcr_ij);
end

% prob_x = x_e_j/sum(x_e_j);
% ent_x = entropy(prob_x)

if virtual==1 % !!!!!!!!!!
    % Virtual pitch array
    x_e_j = conv(x_e_j,fliplr(circshift(y_e_j,[0 -shift_by])),'full');
    % figure(100)
    % plot(x_e_j)

    y_e_j = conv((circshift(y_e_j,[0 -shift_by])),fliplr(circshift(y_e_j,[0 -shift_by])),'full');
    y_e_j = circshift(y_e_j,[0 shift_by]);
    % figure(101)
    % plot(y_e_j)
end

% Distance calculation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isequal(metric, 'l_p')
%    d = sum(abs(x_e_j./(2*norm(x_e_j,p))...
%         - y_e_j./(2*norm(y_e_j,p))).^p)^(1/p);
    d = sum(abs(x_e_j-y_e_j).^p)^(1/p);
elseif isequal(metric, 'cosine')
    d = 1 - (x_e_j * y_e_j'/...
        sqrt((x_e_j * x_e_j') * (y_e_j * y_e_j')));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Plotting
% figure(1), subplot(3,1,1)
% stairs(x_e_j)
% axis([1,q+1,0,max(x_e_j)*1.1])
% set(gca,'XTick',1:q*100/period:q+1)
% set(gca,'XTickLabel',0:period/100)
% 
% figure(1), subplot(3,1,2)
% stairs(y_e_j)
% axis([1,q+1,0,max(y_e_j)*1.1])
% set(gca,'XTick',1:q*100/period:q+1)
% set(gca,'XTickLabel',0:period/100)

clear x_e_j y_e_j

elseif isequal(abs_or_rel,'rel') && r==1
%% RelMonadExp
% This function embeds two pitch collections into relative monad
% expectation scalars and then determines their distance using a metric.

% Expectation array X
x_e = sum(x_w);

% Expectation array Y
y_e = sum(y_w);

% Distance calculation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isequal(metric, 'l_p')
    d = abs(x_e - y_e);
elseif isequal(metric, 'cosine')
    d = 0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elseif isequal(abs_or_rel,'abs') && r==2
%% AbsDyadExp
% This function embeds two pitch collections into absolute dyad expectation 
% matrices and then determines their distance using a metric.

% Expectation array X
X_pcr_j = sum(sparse(X_pcr_ij));
% nonzerosX_pcr_j = find(X_pcr_j);
%
% x1_e_jk_2 = zeros(q, q, 'single');
% for j = nonzerosX_pcr_j
%     for k_2 = nonzerosX_pcr_j
%         x1_e_jk_2(j,k_2) = X_pcr_j(j) * X_pcr_j(k_2);
%     end
% end

x1_e_jk_2 = X_pcr_j' * X_pcr_j;

% nonzerosX_pcr_ij = zeros(d_x, windowlength);
% for i = 1:d_x
%     nonzerosX_pcr_ij(i,:) = find(X_pcr_ij(i,:));
% end

% x2_e_jk_2 = zeros(q, q, 'single');
% x21_e_jk_2 = zeros(q, q, 'single');
% for i = 1:d_x
%     for j = nonzerosX_pcr_ij(i,:)
%         for k_2 = nonzerosX_pcr_ij(i,:)
%             x21_e_jk_2(j,k_2) = X_pcr_ij(i,j) * X_pcr_ij(i,k_2);
%         end
%     end
%     x2_e_jk_2 = x2_e_jk_2 + x21_e_jk_2;
%     x21_e_jk_2 = zeros(q, q, 'single');
% end

j = 1:q;
k_2 = 1:q;
x2_e_jk_2 = sparse(q, q);
for i=1:d_x
    x21_e_jk_2 = sparse(X_pcr_ij(i,j))' * sparse(X_pcr_ij(i,k_2));
    x2_e_jk_2 = x2_e_jk_2 + x21_e_jk_2;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_e_jk_2 = x1_e_jk_2 - x2_e_jk_2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Expectation array Y
Y_pcr_j = sum(sparse(Y_pcr_ij));
% nonzerosY_pcr_j = find(Y_pcr_j);
% 
% y1_e_jk_2 = zeros(q, q, 'single');
% for j = nonzerosY_pcr_j
%     for k_2 = nonzerosY_pcr_j
%         y1_e_jk_2(j,k_2) = Y_pcr_j(j) * Y_pcr_j(k_2);
%     end
% end

y1_e_jk_2 = Y_pcr_j' * Y_pcr_j;

% nonzerosY_pcr_ij = zeros(d_y, windowlength);
% for i = 1:d_y
%     nonzerosY_pcr_ij(i,:) = find(Y_pcr_ij(i,:));
% end
% 
% y2_e_jk_2 = zeros(q, q, 'single');
% y21_e_jk_2 = zeros(q, q, 'single');
% for i = 1:d_y
%     for j = nonzerosY_pcr_ij(i,:)
%         for k_2 = nonzerosY_pcr_ij(i,:)
%             y21_e_jk_2(j,k_2) = Y_pcr_ij(i,j) * Y_pcr_ij(i,k_2);
%         end
%     end
%     y2_e_jk_2 = y2_e_jk_2 + y21_e_jk_2;
%     y21_e_jk_2 = zeros(q, q, 'single');
% end

y2_e_jk_2 = sparse(q, q);
for i=1:d_y
    y21_e_jk_2 = sparse(Y_pcr_ij(i,j))' * sparse(Y_pcr_ij(i,k_2));
    y2_e_jk_2 = y2_e_jk_2 + y21_e_jk_2;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y_e_jk_2 = y1_e_jk_2 - y2_e_jk_2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Distance calculation
vectoredx_e_jk_2 = reshape(x_e_jk_2, 1, q*q);
vectoredy_e_jk_2 = reshape(y_e_jk_2, 1, q*q);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isequal(metric, 'l_p')
    d = sum(abs(vectoredx_e_jk_2./(2*norm(vectoredx_e_jk_2,p))...
        - vectoredy_e_jk_2./(2*norm(vectoredy_e_jk_2,p))).^p)^(1/p);
elseif isequal(metric, 'cosine')
    d = 1 - ((vectoredx_e_jk_2 * vectoredy_e_jk_2')/...
        sqrt((vectoredx_e_jk_2 * vectoredx_e_jk_2')...
        * (vectoredy_e_jk_2 * vectoredy_e_jk_2')));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Plotting
figure(3), subplot(1,2,1)
surf((full(x_e_jk_2)), 'FaceAlpha',1,'LineStyle','none')
axis([1 q+1 1 q+1])
colormap bone
lighting phong
grid off
set(gca,'PlotBoxAspectRatio',[1,1,1])
set(gca,'XTick',1:q*100/period:q+1)
set(gca,'YTick',1:q*100/period:q+1)
set(gca,'XTickLabel',0:period/100)
set(gca,'YTickLabel',0:period/100)

figure(3), subplot(1,2,2)
surface(full(y_e_jk_2), 'FaceAlpha',1,'LineStyle','none')
axis([1 q+1 1 q+1])
colormap bone
lighting phong
grid off
set(gca,'PlotBoxAspectRatio',[1,1,1])
set(gca,'XTick',1:q*100/period:q+1)
set(gca,'YTick',1:q*100/period:q+1)
set(gca,'XTickLabel',0:period/100)
set(gca,'YTickLabel',0:period/100)

clear x_e_jk_2 y_e_jk_2 

elseif isequal(abs_or_rel,'rel') && r==2
%% RelDyadExp
% This function embeds two pitch collections into relative dyad expectation 
% vectors and then determines their distance using a metric.
%
% Note that this routine uses circular convolution (calculated with FFTs) 
% because this is significantly faster than the literal implementation of 
% Eqs(9) and (10) from the Expectation Tensors paper. (The literal
% implentation is included below, but commented out).

% Expectation array X - quick version using circular convolution (FFTs)
X_pcr_j = sum(X_pcr_ij);
clear X_pcr_ij
x1_e_k_2 = real(ifft(fft(X_pcr_j, q).*conj(fft(X_pcr_j, q)), q));

clear X_pcr_j

% for j = 0:q-1
%     if x1_e_k_2(j+1) < 0.0000001
%         x1_e_k_2(j+1) = 0;
%     end
% end

x1_e_k_2(x1_e_k_2<0.0000001) = 0;

if length(gwin) > 1
    longgwin = [gwin zeros(1,q-length(gwin))];
else
    longgwin = [1 zeros(1,q-1)];
end
x2_e_k_2 = real(ifft(fft(longgwin, q).*conj(fft(longgwin, q)), q));

% for j = 0:q-1
%     if x2_e_k_2(j+1) < 0.0000001
%         x2_e_k_2(j+1) = 0;
%     end
% end

x2_e_k_2(x2_e_k_2<0.0000001) = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_e_k_2 = x1_e_k_2 - ((norm(x_w,2)^2) * x2_e_k_2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear x1_e_k_2 x2_e_k_2 

% Expectation array X - slower version
% X_pcr_j = sum(X_pcr_ij);
% nonzerosX_pcr_j = find(X_pcr_j);
% 
% x1_e_k_2 = zeros(1, q);
% x11_e_k_2 = zeros(1, q);
% for j = nonzerosX_pcr_j
%     for k_2 = nonzerosX_pcr_j
%         x11_e_k_2(mod(k_2 - j, q)+1) = X_pcr_j(j) * X_pcr_j(k_2);
%     end
%     x1_e_k_2 = x1_e_k_2 + x11_e_k_2;
%     x11_e_k_2 = zeros(1, q);
% end
% clear X_pcr_j
% clear x11_e_k_2
% clear nonzerosX_pcr_j
% 
% nonzerosX_pcr_ij = zeros(d_x, windowlength);
% for i = 1:d_x
%     nonzerosX_pcr_ij(i,:) = find(X_pcr_ij(i,:));
% end
% 
% x2_e_k_2 = zeros(1, q);
% x21_e_k_2 = zeros(1, q);
% x22_e_k_2 = zeros(1, q);
% for i = 1:d_x
%     for j = nonzerosX_pcr_ij(i,:)
%         for k_2 = nonzerosX_pcr_ij(i,:)
%             x21_e_k_2(mod(k_2 - j, q)+1) = X_pcr_ij(i,j) * X_pcr_ij(i,k_2);
%         end
%         x2_e_k_2 = x2_e_k_2 + x21_e_k_2;
%         x21_e_k_2  = zeros(1, q);
%     end
%     x2_e_k_2 = x2_e_k_2 + x22_e_k_2;
%     x22_e_k_2  = zeros(1, q);
% end
% clear X_pcr_ij
% clear nonzerosX_pcr_ij
% clear x21_e_k_2
% clear x22_e_k_2
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% x_e_k_2 = x1_e_k_2 - x2_e_k_2;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear x1_e_k_2 x2_e_k_2 

% Expectation array Y - quick version using circular convolution (FFTs)
Y_pcr_j = sum(Y_pcr_ij);
clear Y_pcr_ij
y1_e_k_2 = real(ifft(fft(Y_pcr_j, q).*conj(fft(Y_pcr_j, q)), q));
clear Y_pcr_j

% for j = 0:q-1
%     if y1_e_k_2(j+1) < 0.0000001
%         y1_e_k_2(j+1) = 0;
%     end
% end

y1_e_k_2(y1_e_k_2<0.0000001) = 0;

if length(gwin) > 1
    longgwin = [gwin zeros(1,q-length(gwin))];
else
    longgwin = [1 zeros(1,q-1)];
end
y2_e_k_2 = real(ifft(fft(longgwin, q).*conj(fft(longgwin, q)), q));

% for j = 0:q-1
%     if y2_e_k_2(j+1) < 0.0000001
%         y2_e_k_2(j+1) = 0;
%     end
% end

y2_e_k_2(y2_e_k_2<0.0000001) = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y_e_k_2 = y1_e_k_2 - ((norm(y_w,2)^2) * y2_e_k_2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear y1_e_k_2 y2_e_k_2 

% Expectation array Y - slower version
% Y_pcr_j = sum(Y_pcr_ij);
% nonzerosY_pcr_j = find(Y_pcr_j);
% 
% y1_e_k_2 = zeros(1, q);
% y11_e_k_2 = zeros(1, q);
% for j = nonzerosY_pcr_j
%     for k_2 = nonzerosY_pcr_j
%         y11_e_k_2(mod(k_2 - j, q)+1) = Y_pcr_j(j) * Y_pcr_j(k_2);
%     end
%     y1_e_k_2 = y1_e_k_2 + y11_e_k_2;
%     y11_e_k_2 = zeros(1, q);
% end
% clear Y_pcr_j
% clear y11_e_k_2
% clear nonzerosY_pcr_j
% 
% nonzerosY_pcr_ij = zeros(d_y, windowlength);
% for i = 1:d_y
%     nonzerosY_pcr_ij(i,:) = find(Y_pcr_ij(i,:));
% end
% 
% y2_e_k_2 = zeros(1, q);
% y21_e_k_2 = zeros(1, q);
% y22_e_k_2 = zeros(1, q);
% for i = 1:d_y
%     for j = nonzerosY_pcr_ij(i,:)
%         for k_2 = nonzerosY_pcr_ij(i,:)
%             y21_e_k_2(mod(k_2 - j, q)+1) = Y_pcr_ij(i,j) * Y_pcr_ij(i,k_2);
%         end
%         y2_e_k_2 = y2_e_k_2 + y21_e_k_2;
%         y21_e_k_2  = zeros(1, q);
%     end
%     y2_e_k_2 = y2_e_k_2 + y22_e_k_2;
%     y22_e_k_2  = zeros(1, q);
% end
% clear Y_pcr_ij
% clear nonzerosY_pcr_ij
% clear y21_e_k_2
% clear y22_e_k_2
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% y_e_k_2 = y1_e_k_2 - y2_e_k_2;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear y1_e_k_2 y2_e_k_2 

% Distance calculation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isequal(metric, 'l_p')
    d = sum(abs(x_e_k_2./(2*norm(x_e_k_2,p))...
        - y_e_k_2./(2*norm(y_e_k_2,p))).^p)^(1/p);
elseif isequal(metric, 'cosine')
    d = 1 - (x_e_k_2 * y_e_k_2'/...
        sqrt((x_e_k_2 * x_e_k_2') * (y_e_k_2 * y_e_k_2')));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Plotting
figure(4), subplot(3,1,1)
stairs(x_e_k_2)
axis([1,q+1,0,max(x_e_k_2)*1.1])
set(gca,'XTick',1:q*100/period:q+1)
set(gca,'XTickLabel',0:period/100)

figure(4), subplot(3,1,2)
stairs(y_e_k_2)
axis([1,q+1,0,max(y_e_k_2)*1.1])
set(gca,'XTick',1:q*100/period:q+1)
set(gca,'XTickLabel',0:period/100)

clear x_e_k_2 y_e_k_2

elseif isequal(abs_or_rel,'abs') && r==3
%% AbsTriadExp
% This function embeds two pitch collections into absolute triad
% expectation tensors and then determines their distance using a metric.

% Expectation array X
X_pcr_j = sum(X_pcr_ij);
nonzerosX_pcr_j = find(X_pcr_j);

% tic
% x1_e_jk_2k_3 = sparse(q_cubed, 1);
% for j = nonzerosX_pcr_j
%     for k_2 = nonzerosX_pcr_j
%         for k_3 = nonzerosX_pcr_j(nonzerosX_pcr_j >= k_2)
%             x1_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1)...
%                 = X_pcr_j(j) * X_pcr_j(k_2) * X_pcr_j(k_3);
%             x1_e_jk_2k_3(j + q*(k_3 - 1) + (q_squared)*(k_2 - 1), 1)...
%                 = x1_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1);
%         end
%     end
% end
% x_e_jk_2k_3 = x1_e_jk_2k_3;
% clear x1_e_jk_2k_3
% toc

% this is faster
tic
X_pcr_j_sp = squeeze(sptensor(X_pcr_j));
x1_e_jk_2k_3 = ttt(X_pcr_j_sp,ttt(X_pcr_j_sp,X_pcr_j_sp));
x_e_jk_2k_3 = double(x1_e_jk_2k_3);
x_e_jk_2k_3 = x_e_jk_2k_3(:);
clear x1_e_jk_2k_3
toc


nonzerosX_pcr_ij = zeros(d_x, windowlength);
for i = 1:d_x
    nonzerosX_pcr_ij(i,:) = find(X_pcr_ij(i,:));
end

x2_e_jk_2k_3 = sparse(q_cubed, 1);
x21_e_jk_2k_3 = sparse(q_cubed, 1);
for i = 1:d_x
    for j = nonzerosX_pcr_ij(i,:)
        for k_2 = nonzerosX_pcr_ij(i,:)
            for k_3 = nonzerosX_pcr_j(nonzerosX_pcr_j >= k_2)
                x21_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1)...
                    = X_pcr_ij(i,j) * X_pcr_ij(i,k_2) * X_pcr_j(k_3);
                x21_e_jk_2k_3(j + q*(k_3 - 1) + (q_squared)*(k_2 - 1), 1)...
                    = x21_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1);
            end
        end
    end
    x2_e_jk_2k_3 = x2_e_jk_2k_3 + x21_e_jk_2k_3;
    x21_e_jk_2k_3 = sparse(q_cubed, 1);
end
x_e_jk_2k_3 = x_e_jk_2k_3 - x2_e_jk_2k_3;
clear x2_e_jk_2k_3

% % the below should be faster, but not sure how to get the i index thing 
% % happening...
% X_pcr_ij_sp = sptensor(X_pcr_ij);
% x21_e_jk_2k_3 = cell(d_x,1);
% for i=1:d_x
%     x21_e_jk_2k_3{i,:}...
%         = ttt(ttt(X_pcr_ij_sp(i,:),X_pcr_ij_sp(i,:)),X_pcr_j_sp);
% end
% for i=1:d_x
% x21_e_jk_2k_3(i,:,:,:) = x21_e_jk_2k_3{i};
% end
% stop
% x21_e_jk_2k_3 = collapse(x21_e_jk_2k_3,4);
% x21_e_jk_2k_3 = x21_e_jk_2k_3(:);
% stop

% another go...
% tic
% X_pcr_ij_sp = sptensor(X_pcr_ij);
% upper_simplex = sptensor([],[],[q q q]);
% for j=1:q
%     for k_2=j:q
%         for k_3=k_2:q
%             upper_simplex(j, k_2, k_3)=1;
%         end
%     end
% end
% toc
% x2_e_jk_2k_3
% forget it - ridiculosly slow way of generating an upper simplex logical 
% (binary) tensor


x3_e_jk_2k_3 = sparse(q_cubed, 1);
x31_e_jk_2k_3 = sparse(q_cubed, 1);
for i = 1:d_x
    for j = nonzerosX_pcr_ij(i,:)
        for k_2 = nonzerosX_pcr_j
            for k_3 = nonzerosX_pcr_ij(i, nonzerosX_pcr_ij(i,:) >= k_2)
                x31_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1)...
                    = X_pcr_ij(i,j) * X_pcr_ij(i,k_3) * X_pcr_j(k_2);
                x31_e_jk_2k_3(j + q*(k_3 - 1) + (q_squared)*(k_2 - 1), 1)...
                    = x31_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1);
            end
        end
    end
    x3_e_jk_2k_3 = x3_e_jk_2k_3 + x31_e_jk_2k_3;
    x31_e_jk_2k_3 = sparse(q_cubed, 1);
end
x_e_jk_2k_3 = x_e_jk_2k_3 - x3_e_jk_2k_3;
clear x3_e_jk_2k_3

x4_e_jk_2k_3 = sparse(q_cubed, 1);
x41_e_jk_2k_3 = sparse(q_cubed, 1);
for i = 1:d_x
    for j = nonzerosX_pcr_j
        for k_2 = nonzerosX_pcr_ij(i,:)
            for k_3 = nonzerosX_pcr_ij(i, nonzerosX_pcr_ij(i,:) >= k_2)
                x41_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1)...
                    = X_pcr_ij(i,k_2) * X_pcr_ij(i,k_3) * X_pcr_j(j);
                x41_e_jk_2k_3(j + q*(k_3 - 1) + (q_squared)*(k_2 - 1), 1)...
                    = x41_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1);
            end
        end
    end
    x4_e_jk_2k_3 = x4_e_jk_2k_3 + x41_e_jk_2k_3;
    x41_e_jk_2k_3 = sparse(q_cubed, 1);
end
x_e_jk_2k_3 = x_e_jk_2k_3 - x4_e_jk_2k_3;
clear x4_e_jk_2k_3

x5_e_jk_2k_3 = sparse(q_cubed, 1);
x51_e_jk_2k_3 = sparse(q_cubed, 1);
for i = 1:d_x
    for j = nonzerosX_pcr_ij(i,:)
        for k_2 = nonzerosX_pcr_ij(i,:)
            for k_3 = nonzerosX_pcr_ij(i, nonzerosX_pcr_ij(i,:) >= k_2)
                x51_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1)...
                    = X_pcr_ij(i,j) * X_pcr_ij(i,k_2) * X_pcr_ij(i,k_3);
                x51_e_jk_2k_3(j + q*(k_3 - 1) + (q_squared)*(k_2 - 1), 1)...
                    = x51_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1);
            end
        end
    end
    x5_e_jk_2k_3 = x5_e_jk_2k_3 + x51_e_jk_2k_3;
    clear x51_e_jk_2k_3
    x51_e_jk_2k_3 = sparse(q_cubed, 1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_e_jk_2k_3 = x_e_jk_2k_3 + (2*x5_e_jk_2k_3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear x5_e_jk_2k_3

test_x = squeeze(sum(reshape(full(x_e_jk_2k_3), [q, q, q]), 3));
figure(5), subplot(1,2,1)
surf(test_x, 'FaceAlpha',1,'LineStyle','none')
axis([1 q+1 1 q+1])
colormap bone
lighting phong
grid off
set(gca,'PlotBoxAspectRatio',[1,1,1])
set(gca,'XTick',1:q*100/period:q+1)
set(gca,'YTick',1:q*100/period:q+1)
set(gca,'XTickLabel',0:period/100)
set(gca,'YTickLabel',0:period/100)
clear test_x

% Expectation array Y
Y_pcr_j = sum(Y_pcr_ij);
nonzerosY_pcr_j = find(Y_pcr_j);

y1_e_jk_2k_3 = sparse(q_cubed, 1);
for j = nonzerosY_pcr_j
    for k_2 = nonzerosY_pcr_j
        for k_3 = nonzerosY_pcr_j(nonzerosY_pcr_j >= k_2)
            y1_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1)...
                = Y_pcr_j(j) * Y_pcr_j(k_2) * Y_pcr_j(k_3);
            y1_e_jk_2k_3(j + q*(k_3 - 1) + (q_squared)*(k_2 - 1), 1)...
                = y1_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1);
        end
    end
end
y_e_jk_2k_3 = y1_e_jk_2k_3;
clear y1_e_jk_2k_3

nonzerosY_pcr_ij = zeros(d_y, windowlength);
for i = 1:d_y
    nonzerosY_pcr_ij(i,:) = find(Y_pcr_ij(i,:));
end

y2_e_jk_2k_3 = sparse(q_cubed, 1);
y21_e_jk_2k_3 = sparse(q_cubed, 1);
for i = 1:d_y
    for j = nonzerosY_pcr_ij(i,:)
        for k_2 = nonzerosY_pcr_ij(i,:)
            for k_3 = nonzerosY_pcr_j(nonzerosY_pcr_j >= k_2)
                y21_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1)...
                    = Y_pcr_ij(i,j) * Y_pcr_ij(i,k_2) * Y_pcr_j(k_3);
                y21_e_jk_2k_3(j + q*(k_3 - 1) + (q_squared)*(k_2 - 1), 1)...
                    = y21_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1);
            end
        end
    end
    y2_e_jk_2k_3 = y2_e_jk_2k_3 + y21_e_jk_2k_3;
    y21_e_jk_2k_3 = sparse(q_cubed, 1);
end
y_e_jk_2k_3 = y_e_jk_2k_3 - y2_e_jk_2k_3;
clear y2_e_jk_2k_3

y3_e_jk_2k_3 = sparse(q_cubed, 1);
y31_e_jk_2k_3 = sparse(q_cubed, 1);
for i = 1:d_y
    for j = nonzerosY_pcr_ij(i,:)
        for k_2 = nonzerosY_pcr_j
            for k_3 = nonzerosY_pcr_ij(i, nonzerosY_pcr_ij(i,:) >= k_2)
                y31_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1)...
                    = Y_pcr_ij(i,j) * Y_pcr_ij(i,k_3) * Y_pcr_j(k_2);
                y31_e_jk_2k_3(j + q*(k_3 - 1) + (q_squared)*(k_2 - 1), 1)...
                    = y31_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1);
            end
        end
    end
    y3_e_jk_2k_3 = y3_e_jk_2k_3 + y31_e_jk_2k_3;
    y31_e_jk_2k_3 = sparse(q_cubed, 1);
end
y_e_jk_2k_3 = y_e_jk_2k_3 - y3_e_jk_2k_3;
clear y3_e_jk_2k_3

y4_e_jk_2k_3 = sparse(q_cubed, 1);
y41_e_jk_2k_3 = sparse(q_cubed, 1);
for i = 1:d_y
    for j = nonzerosY_pcr_j
        for k_2 = nonzerosY_pcr_ij(i,:)
            for k_3 = nonzerosY_pcr_ij(i, nonzerosY_pcr_ij(i,:) >= k_2)
                y41_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1)...
                    = Y_pcr_ij(i,k_2) * Y_pcr_ij(i,k_3) * Y_pcr_j(j);
                y41_e_jk_2k_3(j + q*(k_3 - 1) + (q_squared)*(k_2 - 1), 1)...
                    = y41_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1);
            end
        end
    end
    y4_e_jk_2k_3 = y4_e_jk_2k_3 + y41_e_jk_2k_3;
    y41_e_jk_2k_3 = sparse(q_cubed, 1);
end
y_e_jk_2k_3 = y_e_jk_2k_3 - y4_e_jk_2k_3;
clear y4_e_jk_2k_3

y5_e_jk_2k_3 = sparse(q_cubed, 1);
y51_e_jk_2k_3 = sparse(q_cubed, 1);
for i = 1:d_y
    for j = nonzerosY_pcr_ij(i,:)
        for k_2 = nonzerosY_pcr_ij(i,:)
            for k_3 = nonzerosY_pcr_ij(i, nonzerosY_pcr_ij(i,:) >= k_2)
                y51_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1)...
                    = Y_pcr_ij(i,j) * Y_pcr_ij(i,k_2) * Y_pcr_ij(i,k_3);
                y51_e_jk_2k_3(j + q*(k_3 - 1) + (q_squared)*(k_2 - 1), 1)...
                    = y51_e_jk_2k_3(j + q*(k_2 - 1) + (q_squared)*(k_3 - 1), 1);
            end
        end
    end
    y5_e_jk_2k_3 = y5_e_jk_2k_3 + y51_e_jk_2k_3;
    clear y51_e_jk_2k_3
    y51_e_jk_2k_3 = sparse(q_cubed, 1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y_e_jk_2k_3 = y_e_jk_2k_3 + (2*y5_e_jk_2k_3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear y5_e_jk_2k_3

test_y = squeeze(sum(reshape(full(y_e_jk_2k_3), [q, q, q]), 3));
figure(5), subplot(1,2,2)
surf(test_y, 'FaceAlpha',1,'LineStyle','none')
axis([1 q+1 1 q+1])
colormap bone
lighting phong
grid off
set(gca,'PlotBoxAspectRatio',[1,1,1])
set(gca,'XTick',1:q*100/period:q+1)
set(gca,'YTick',1:q*100/period:q+1)
set(gca,'XTickLabel',0:period/100)
set(gca,'YTickLabel',0:period/100)
clear test_y

% Distance calculation
p_norm_x_e_jk_2k_3 = sum(x_e_jk_2k_3.^p)^(1/p);
p_norm_y_e_jk_2k_3 = sum(y_e_jk_2k_3.^p)^(1/p);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isequal(metric, 'l_p')
    d = sum(abs(x_e_jk_2k_3./(2*p_norm_x_e_jk_2k_3)... 
        - y_e_jk_2k_3./(2*p_norm_y_e_jk_2k_3)).^p)^(1/p);
elseif isequal(metric, 'cosine')
    d = 1 - ((x_e_jk_2k_3' * y_e_jk_2k_3)/...
        sqrt((x_e_jk_2k_3' * x_e_jk_2k_3) * (y_e_jk_2k_3' * y_e_jk_2k_3)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear x_e_jk_2k_3 y_e_jk_2k_3 

elseif isequal(abs_or_rel,'rel') && r==3
%% RelTriadExp
% This function embeds two pitch collections into relative triad
% expectation matrices and then determines their distance using a metric.

% Expectation array X
X_pcr_j = sum(X_pcr_ij);
nonzerosX_pcr_j = find(X_pcr_j);
nonzerosX_pcr_ij = zeros(d_x, windowlength);
for i = 1:d_x
    nonzerosX_pcr_ij(i,:) = find(X_pcr_ij(i,:));
end

x1_e_k_2k_3 = zeros(q, q);
x11_e_k_2k_3 = zeros(q, q);
for j = nonzerosX_pcr_j
    for k_2 = nonzerosX_pcr_j
        for k_3 = nonzerosX_pcr_j(nonzerosX_pcr_j >= k_2)
            x11_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1)...
                = X_pcr_j(j) * X_pcr_j(k_2) * X_pcr_j(k_3);
            x11_e_k_2k_3(mod(k_3-j, q)+1, mod(k_2-j, q)+1)...
                = x11_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1);
        end
    end
    x1_e_k_2k_3 = x1_e_k_2k_3 + x11_e_k_2k_3;
    x11_e_k_2k_3 = zeros(q, q);
end
x_e_k_2k_3 = x1_e_k_2k_3;
clear x1_e_k_2k_3

x2_e_k_2k_3 = zeros(q, q);
x21_e_k_2k_3 = zeros(q, q);
x22_e_k_2k_3 = zeros(q, q);
for i = 1:d_x
    for j = nonzerosX_pcr_ij(i,:)
        for k_2 = nonzerosX_pcr_ij(i,:)
            for k_3 = nonzerosX_pcr_j(nonzerosX_pcr_j >= k_2)
                x21_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1)...
                    = X_pcr_ij(i,j) * X_pcr_ij(i,k_2) * X_pcr_j(k_3);
                x21_e_k_2k_3(mod(k_3-j, q)+1, mod(k_2-j, q)+1)...
                    = x21_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1);
            end
        end
        x2_e_k_2k_3 = x2_e_k_2k_3 + x21_e_k_2k_3;
        x21_e_k_2k_3 = zeros(q, q);
    end
    x2_e_k_2k_3 = x2_e_k_2k_3 + x22_e_k_2k_3;
    x22_e_k_2k_3 = zeros(q, q);
end
x_e_k_2k_3 = x_e_k_2k_3 - x2_e_k_2k_3;
clear x2_e_k_2k_3

x3_e_k_2k_3 = zeros(q, q);
x31_e_k_2k_3 = zeros(q, q);
x32_e_k_2k_3 = zeros(q, q);
for i = 1:d_x
    for j = nonzerosX_pcr_ij(i,:)
        for k_2 = nonzerosX_pcr_j
            for k_3 = nonzerosX_pcr_ij(i, nonzerosX_pcr_ij(i,:) >= k_2)
                x31_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1)...
                    = X_pcr_ij(i,j) * X_pcr_ij(i,k_3) * X_pcr_j(k_2);
                x31_e_k_2k_3(mod(k_3-j, q)+1, mod(k_2-j, q)+1)...
                    = x31_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1);
            end
        end
        x3_e_k_2k_3 = x3_e_k_2k_3 + x31_e_k_2k_3;
        x31_e_k_2k_3 = zeros(q, q);
    end
    x3_e_k_2k_3 = x3_e_k_2k_3 + x32_e_k_2k_3;
    x32_e_k_2k_3 = zeros(q, q);
end
x_e_k_2k_3 = x_e_k_2k_3 - x3_e_k_2k_3;
clear x3_e_k_2k_3

x4_e_k_2k_3 =zeros(q, q);
x41_e_k_2k_3 =zeros(q, q);
x42_e_k_2k_3 =zeros(q, q);
for i = 1:d_x
    for j = nonzerosX_pcr_j
        for k_2 = nonzerosX_pcr_ij(i,:)
            for k_3 = nonzerosX_pcr_ij(i, nonzerosX_pcr_ij(i,:) >= k_2)
                x41_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1)...
                    = X_pcr_ij(i,k_2) * X_pcr_ij(i,k_3) * X_pcr_j(j);
                x41_e_k_2k_3(mod(k_3-j, q)+1, mod(k_2-j, q)+1)...
                    = x41_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1);
            end
        end
        x4_e_k_2k_3 = x4_e_k_2k_3 + x41_e_k_2k_3;
        x41_e_k_2k_3 = zeros(q, q);
    end
    x4_e_k_2k_3 = x4_e_k_2k_3 + x42_e_k_2k_3;
    x42_e_k_2k_3 = zeros(q, q);
end
x_e_k_2k_3 = x_e_k_2k_3 - x4_e_k_2k_3;
clear x4_e_k_2k_3

x5_e_k_2k_3 = zeros(q, q);
x51_e_k_2k_3 = zeros(q, q);
x52_e_k_2k_3 = zeros(q, q);
for i = 1:d_x
    for j = nonzerosX_pcr_ij(i,:)
        for k_2 = nonzerosX_pcr_ij(i,:)
            for k_3 = nonzerosX_pcr_ij(i, nonzerosX_pcr_ij(i,:) >= k_2)
                x51_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1)...
                    = X_pcr_ij(i,j) * X_pcr_ij(i,k_2) * X_pcr_ij(i,k_3);
                x51_e_k_2k_3(mod(k_3-j, q)+1, mod(k_2-j, q)+1)...
                    = x51_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1);
            end
        end
        x5_e_k_2k_3 = x5_e_k_2k_3 + x51_e_k_2k_3;
        x51_e_k_2k_3 = zeros(q, q);
    end
    x5_e_k_2k_3 = x5_e_k_2k_3 + x52_e_k_2k_3;
    x52_e_k_2k_3 = zeros(q, q);
end
clear X_pcr_ij
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_e_k_2k_3 = x_e_k_2k_3 + (2 * x5_e_k_2k_3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear x5_e_k_2k_3

% Expectation array Y
Y_pcr_j = sum(Y_pcr_ij);
nonzerosY_pcr_j = find(Y_pcr_j);
nonzerosY_pcr_ij = zeros(d_y, windowlength);
for i = 1:d_y
    nonzerosY_pcr_ij(i,:) = find(Y_pcr_ij(i,:));
end

y1_e_k_2k_3 = zeros(q, q);
y11_e_k_2k_3 = zeros(q, q);
for j = nonzerosY_pcr_j
    for k_2 = nonzerosY_pcr_j
        for k_3 = nonzerosY_pcr_j(nonzerosY_pcr_j >= k_2)
            y11_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1)...
                = Y_pcr_j(j) * Y_pcr_j(k_2) * Y_pcr_j(k_3);
            y11_e_k_2k_3(mod(k_3-j, q)+1, mod(k_2-j, q)+1)...
                = y11_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1);
        end
    end
    y1_e_k_2k_3 = y1_e_k_2k_3 + y11_e_k_2k_3;
    y11_e_k_2k_3 = zeros(q, q);
end
y_e_k_2k_3 = y1_e_k_2k_3;
clear y1_e_k_2k_3

y2_e_k_2k_3 = zeros(q, q);
y21_e_k_2k_3 = zeros(q, q);
y22_e_k_2k_3 = zeros(q, q);
for i = 1:d_y
    for j = nonzerosY_pcr_ij(i,:)
        for k_2 = nonzerosY_pcr_ij(i,:)
            for k_3 = nonzerosY_pcr_j(nonzerosY_pcr_j >= k_2)
                y21_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1)...
                    = Y_pcr_ij(i,j) * Y_pcr_ij(i,k_2) * Y_pcr_j(k_3);
                y21_e_k_2k_3(mod(k_3-j, q)+1, mod(k_2-j, q)+1)...
                    = y21_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1);
            end
        end
        y2_e_k_2k_3 = y2_e_k_2k_3 + y21_e_k_2k_3;
        y21_e_k_2k_3 = zeros(q, q);
    end
    y2_e_k_2k_3 = y2_e_k_2k_3 + y22_e_k_2k_3;
    y22_e_k_2k_3 = zeros(q, q);
end
y_e_k_2k_3 = y_e_k_2k_3 - y2_e_k_2k_3;
clear y2_e_k_2k_3

y3_e_k_2k_3 = zeros(q, q);
y31_e_k_2k_3 = zeros(q, q);
y32_e_k_2k_3 = zeros(q, q);
for i = 1:d_y
    for j = nonzerosY_pcr_ij(i,:)
        for k_2 = nonzerosY_pcr_j
            for k_3 = nonzerosY_pcr_ij(i, nonzerosY_pcr_ij(i,:) >= k_2)
                y31_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1)...
                    = Y_pcr_ij(i,j) * Y_pcr_ij(i,k_3) * Y_pcr_j(k_2);
                y31_e_k_2k_3(mod(k_3-j, q)+1, mod(k_2-j, q)+1)...
                    = y31_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1);
            end
        end
        y3_e_k_2k_3 = y3_e_k_2k_3 + y31_e_k_2k_3;
        y31_e_k_2k_3 = zeros(q, q);
    end
    y3_e_k_2k_3 = y3_e_k_2k_3 + y32_e_k_2k_3;
    y32_e_k_2k_3 = zeros(q, q);
end
y_e_k_2k_3 = y_e_k_2k_3 - y3_e_k_2k_3;
clear y3_e_k_2k_3

y4_e_k_2k_3 =zeros(q, q);
y41_e_k_2k_3 =zeros(q, q);
y42_e_k_2k_3 =zeros(q, q);
for i = 1:d_y
    for j = nonzerosY_pcr_j
        for k_2 = nonzerosY_pcr_ij(i,:)
            for k_3 = nonzerosY_pcr_ij(i, nonzerosY_pcr_ij(i,:) >= k_2)
                y41_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1)...
                    = Y_pcr_ij(i,k_2) * Y_pcr_ij(i,k_3) * Y_pcr_j(j);
                y41_e_k_2k_3(mod(k_3-j, q)+1, mod(k_2-j, q)+1)...
                    = y41_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1);
            end
        end
        y4_e_k_2k_3 = y4_e_k_2k_3 + y41_e_k_2k_3;
        y41_e_k_2k_3 = zeros(q, q);
    end
    y4_e_k_2k_3 = y4_e_k_2k_3 + y42_e_k_2k_3;
    y42_e_k_2k_3 = zeros(q, q);
end
y_e_k_2k_3 = y_e_k_2k_3 - y4_e_k_2k_3;
clear y4_e_k_2k_3

y5_e_k_2k_3 = zeros(q, q);
y51_e_k_2k_3 = zeros(q, q);
y52_e_k_2k_3 = zeros(q, q);
for i = 1:d_y
    for j = nonzerosY_pcr_ij(i,:)
        for k_2 = nonzerosY_pcr_ij(i,:)
            for k_3 = nonzerosY_pcr_ij(i, nonzerosY_pcr_ij(i,:) >= k_2)
                y51_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1)...
                    = Y_pcr_ij(i,j) * Y_pcr_ij(i,k_2) * Y_pcr_ij(i,k_3);
                y51_e_k_2k_3(mod(k_3-j, q)+1, mod(k_2-j, q)+1)...
                    = y51_e_k_2k_3(mod(k_2-j, q)+1, mod(k_3-j, q)+1);
            end
        end
        y5_e_k_2k_3 = y5_e_k_2k_3 + y51_e_k_2k_3;
        y51_e_k_2k_3 = zeros(q, q);
    end
    y5_e_k_2k_3 = y5_e_k_2k_3 + y52_e_k_2k_3;
    y52_e_k_2k_3 = zeros(q, q);
end
clear Y_pcr_ij
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y_e_k_2k_3 = y_e_k_2k_3 + (2 * y5_e_k_2k_3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear y5_e_k_2k_3

% Distance calculation
vectored_x_e_k_2k_3 = reshape(x_e_k_2k_3, 1, q*q);
vectored_y_e_k_2k_3 = reshape(y_e_k_2k_3, 1, q*q);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isequal(metric, 'l_p')
    d = sum(abs(vectored_x_e_k_2k_3./(2*norm(vectored_x_e_k_2k_3,p))...
        - vectored_y_e_k_2k_3./(2*norm(vectored_y_e_k_2k_3,p))).^p)^(1/p);
elseif isequal(metric, 'cosine')
    d = 1 - ((vectored_x_e_k_2k_3 * vectored_y_e_k_2k_3')/...
        sqrt((vectored_x_e_k_2k_3 * vectored_x_e_k_2k_3')...
        * (vectored_y_e_k_2k_3 * vectored_y_e_k_2k_3')));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Plotting
hold off
figure(6), subplot(1,2,1)
surf((double(x_e_k_2k_3)), 'FaceAlpha',1,'LineStyle','none')
axis([1 q+1 1 q+1])
%axis([1 q+1 1 q+1 0 max(max(x_e_k_2k_3)) 0 max(max(x_e_k_2k_3))])
colormap bone
lighting phong
grid off
set(gca,'PlotBoxAspectRatio',[1,1,1])
set(gca,'XTick',1:q*100/period:q+1)
set(gca,'YTick',1:q*100/period:q+1)
set(gca,'XTickLabel',0:period/100)
set(gca,'YTickLabel',0:period/100)

hold off
figure(6), subplot(1,2,2)
surface(double(y_e_k_2k_3), 'FaceAlpha',1,'LineStyle','none')
axis([1 q+1 1 q+1])
colormap bone
lighting phong
grid off
set(gca,'PlotBoxAspectRatio',[1,1,1])
set(gca,'XTick',1:q*100/period:q+1)
set(gca,'YTick',1:q*100/period:q+1)
set(gca,'XTickLabel',0:period/100)
set(gca,'YTickLabel',0:period/100)

clear x_e_k_2k_3 y_e_k_2k_3

elseif isequal(abs_or_rel,'abs') && r==4
%% AbsTetradExp

% This function embeds two pitch collections into relative tetrad
% expectation tensors and then determines their distance using a metric.

% This routine can probably be sped up. When RelTriad is written in the
% same way as this routine (see ALTRelTriadExp.m), it is significantly
% slower than the implementation I actually use above. It would, however,
% be very long-winded to write the RelTetrad code in the same way as the
% RelTriad code.
one_x = ones(1,d_x);
X_pcr_j = sum(X_pcr_ij);
nonzerosX_pcr_j = find(X_pcr_j);

x_e_jk_2k_3k_4 = sparse(q_fourth, 1);
for j = nonzerosX_pcr_j
    for k_2 = nonzerosX_pcr_j
        for k_3 = nonzerosX_pcr_j(nonzerosX_pcr_j >= k_2);
            for k_4 = nonzerosX_pcr_j(nonzerosX_pcr_j >= k_3);
                one_xA = X_pcr_j(j);
                one_xB = X_pcr_j(k_2);
                one_xC = X_pcr_j(k_3);
                one_xD = X_pcr_j(k_4);
                one_xAB = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_2));
                one_xAC = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_3));
                one_xAD = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_4));
                one_xBC = one_x * (X_pcr_ij(:,k_2).*X_pcr_ij(:,k_3));
                one_xBD = one_x * (X_pcr_ij(:,k_2).*X_pcr_ij(:,k_4));
                one_xCD = one_x * (X_pcr_ij(:,k_3).*X_pcr_ij(:,k_4));
                one_xABC = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_2).*X_pcr_ij(:,k_3));
                one_xABD = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_2).*X_pcr_ij(:,k_4));
                one_xACD = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_3).*X_pcr_ij(:,k_4));
                one_xBCD = one_x * (X_pcr_ij(:,k_2).*X_pcr_ij(:,k_3).*X_pcr_ij(:,k_4));
                one_xABCD = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_2).*X_pcr_ij(:,k_3).*X_pcr_ij(:,k_4));
                
                test_x = one_xA * one_xB * one_xC * one_xD...
                    -   ( (one_xAB * one_xC * one_xD)...
                        + (one_xAC * one_xB * one_xD)...
                        + (one_xAD * one_xB * one_xC)...
                        + (one_xBC * one_xA * one_xD)...
                        + (one_xBD * one_xA * one_xC)...
                        + (one_xCD * one_xA * one_xB) )...
                    + 2*( (one_xABC * one_xD)...
                        + (one_xABD * one_xC)...
                        + (one_xACD * one_xB)...
                        + (one_xBCD * one_xA) )...
                    +   ( (one_xAB * one_xCD)...
                        + (one_xAC * one_xBD)...
                        + (one_xAD * one_xBC) )...
                    - 6*(  one_xABCD );
                if test_x > 0.0001
                    % precalc mods?
                    x_e_jk_2k_3k_4(mod(j-1 + q*(k_2-1) + q_squared*(k_3-1) + q_cubed*(k_4-1), q_fourth)+1, 1) = test_x;
                    x_e_jk_2k_3k_4(mod(j-1 + q*(k_2-1) + q_squared*(k_4-1) + q_cubed*(k_3-1), q_fourth)+1, 1) = test_x;
                    x_e_jk_2k_3k_4(mod(j-1 + q*(k_3-1) + q_squared*(k_2-1) + q_cubed*(k_4-1), q_fourth)+1, 1) = test_x;
                    x_e_jk_2k_3k_4(mod(j-1 + q*(k_3-1) + q_squared*(k_4-1) + q_cubed*(k_2-1), q_fourth)+1, 1) = test_x;
                    x_e_jk_2k_3k_4(mod(j-1 + q*(k_4-1) + q_squared*(k_2-1) + q_cubed*(k_3-1), q_fourth)+1, 1) = test_x;
                    x_e_jk_2k_3k_4(mod(j-1 + q*(k_4-1) + q_squared*(k_3-1) + q_cubed*(k_2-1), q_fourth)+1, 1) = test_x;
                end
            end
        end
    end
end

% test_x = squeeze(sum(sum(reshape(full(x_e_jk_2k_3k_4), [q q q q]), 4),3));
% figure(7), subplot(1, 2, 1)
% surf(test_x, 'FaceAlpha',1,'LineStyle','none')
% axis([1 q+1 1 q+1])
% colormap bone
% lighting phong
% grid off
% set(gca,'PlotBoxAspectRatio',[1,1,1])
% set(gca,'XTick',1:q*100/period:q+1)
% set(gca,'YTick',1:q*100/period:q+1)
% set(gca,'XTickLabel',0:period/100)
% set(gca,'YTickLabel',0:period/100)
% clear test_x

% Expectation array Y
Y_pcr_j = sum(Y_pcr_ij);
nonzerosY_pcr_j = find(Y_pcr_j);
one_y = ones(1,d_y);

y_e_jk_2k_3k_4 = sparse(q_fourth, 1);
for j = nonzerosY_pcr_j
    for k_2 = nonzerosY_pcr_j
        for k_3 = nonzerosY_pcr_j(nonzerosY_pcr_j >= k_2);
            for k_4 = nonzerosY_pcr_j(nonzerosY_pcr_j >= k_3);
                one_yA = Y_pcr_j(j);
                one_yB = Y_pcr_j(k_2);
                one_yC = Y_pcr_j(k_3);
                one_yD = Y_pcr_j(k_4);
                one_yAB = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_2));
                one_yAC = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_3));
                one_yAD = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_4));
                one_yBC = one_y * (Y_pcr_ij(:,k_2).*Y_pcr_ij(:,k_3));
                one_yBD = one_y * (Y_pcr_ij(:,k_2).*Y_pcr_ij(:,k_4));
                one_yCD = one_y * (Y_pcr_ij(:,k_3).*Y_pcr_ij(:,k_4));
                one_yABC = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_2).*Y_pcr_ij(:,k_3));
                one_yABD = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_2).*Y_pcr_ij(:,k_4));
                one_yACD = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_3).*Y_pcr_ij(:,k_4));
                one_yBCD = one_y * (Y_pcr_ij(:,k_2).*Y_pcr_ij(:,k_3).*Y_pcr_ij(:,k_4));
                one_yABCD = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_2).*Y_pcr_ij(:,k_3).*Y_pcr_ij(:,k_4));
                
                test_y = one_yA * one_yB * one_yC * one_yD...
                    -   ( (one_yAB * one_yC * one_yD)...
                        + (one_yAC * one_yB * one_yD)...
                        + (one_yAD * one_yB * one_yC)...
                        + (one_yBC * one_yA * one_yD)...
                        + (one_yBD * one_yA * one_yC)...
                        + (one_yCD * one_yA * one_yB) )...
                    + 2*( (one_yABC * one_yD)...
                        + (one_yABD * one_yC)...
                        + (one_yACD * one_yB)...
                        + (one_yBCD * one_yA) )...
                    +   ( (one_yAB * one_yCD)...
                        + (one_yAC * one_yBD)...
                        + (one_yAD * one_yBC) )...
                    - 6*(  one_yABCD );
                if test_y > 0.0001
                    % precalc mods?
                    y_e_jk_2k_3k_4(mod(j-1 + q*(k_2-1) + q_squared*(k_3-1) + q_cubed*(k_4-1), q_fourth)+1, 1) = test_y;
                    y_e_jk_2k_3k_4(mod(j-1 + q*(k_2-1) + q_squared*(k_4-1) + q_cubed*(k_3-1), q_fourth)+1, 1) = test_y;
                    y_e_jk_2k_3k_4(mod(j-1 + q*(k_3-1) + q_squared*(k_2-1) + q_cubed*(k_4-1), q_fourth)+1, 1) = test_y;
                    y_e_jk_2k_3k_4(mod(j-1 + q*(k_3-1) + q_squared*(k_4-1) + q_cubed*(k_2-1), q_fourth)+1, 1) = test_y;
                    y_e_jk_2k_3k_4(mod(j-1 + q*(k_4-1) + q_squared*(k_2-1) + q_cubed*(k_3-1), q_fourth)+1, 1) = test_y;
                    y_e_jk_2k_3k_4(mod(j-1 + q*(k_4-1) + q_squared*(k_3-1) + q_cubed*(k_2-1), q_fourth)+1, 1) = test_y;
                end
            end
        end
    end
end

% test_y = squeeze(sum(sum(reshape(full(y_e_jk_2k_3k_4), [q q q q]), 4),3));
% figure(7), subplot(1, 2, 2)
% surf(test_y, 'FaceAlpha',1,'LineStyle','none')
% axis([1 q+1 1 q+1])
% colormap bone
% lighting phong
% grid off
% set(gca,'PlotBoxAspectRatio',[1,1,1])
% set(gca,'XTick',1:q*100/period:q+1)
% set(gca,'YTick',1:q*100/period:q+1)
% set(gca,'XTickLabel',0:period/100)
% set(gca,'YTickLabel',0:period/100)
% clear test_y

% Distance calculation
p_norm_x_e_jk_2k_3k_4 = sum(x_e_jk_2k_3k_4.^p)^(1/p);
p_norm_y_e_jk_2k_3k_4 = sum(y_e_jk_2k_3k_4.^p)^(1/p);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isequal(metric, 'l_p')
    d = sum(abs(x_e_jk_2k_3k_4./(2*p_norm_x_e_jk_2k_3k_4)...
        - y_e_jk_2k_3k_4./(2*p_norm_y_e_jk_2k_3k_4)).^p)^(1/p);
elseif isequal(metric, 'cosine')
    d = 1 - ((x_e_jk_2k_3k_4' * y_e_jk_2k_3k_4)/...
        sqrt((x_e_jk_2k_3k_4' * x_e_jk_2k_3k_4)...
        * (y_e_jk_2k_3k_4' * y_e_jk_2k_3k_4)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear x_e_jk_2k_3k_4 y_e_jk_2k_3k_4

elseif isequal(abs_or_rel,'rel') && r==4
%% RelTetradExp
% This function embeds two pitch collections into relative tetrad
% expectation tensors and then determines their distance using a metric.

% This routine can probably be sped up. If RelTriad is written in the same
% way as this, it is significantly slower than the implementation I've used
% above. It would, however, be very long-winded to write the RelTetrad code
% in the same way as the RelTriad code.
one_x = ones(1,d_x);
X_pcr_j = sum(X_pcr_ij);
nonzerosX_pcr_j = find(X_pcr_j);

x_e_k_2k_3k_4 = sparse(q_cubed, 1);
x1_e_k_2k_3k_4 = sparse(q_cubed, 1);
for j = nonzerosX_pcr_j
    for k_2 = nonzerosX_pcr_j
        for k_3 = nonzerosX_pcr_j(nonzerosX_pcr_j >= k_2);
            for k_4 = nonzerosX_pcr_j(nonzerosX_pcr_j >= k_3);
                one_xA = X_pcr_j(j);
                one_xB = X_pcr_j(k_2);
                one_xC = X_pcr_j(k_3);
                one_xD = X_pcr_j(k_4);
                one_xAB = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_2));
                one_xAC = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_3));
                one_xAD = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_4));
                one_xBC = one_x * (X_pcr_ij(:,k_2).*X_pcr_ij(:,k_3));
                one_xBD = one_x * (X_pcr_ij(:,k_2).*X_pcr_ij(:,k_4));
                one_xCD = one_x * (X_pcr_ij(:,k_3).*X_pcr_ij(:,k_4));
                one_xABC = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_2).*X_pcr_ij(:,k_3));
                one_xABD = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_2).*X_pcr_ij(:,k_4));
                one_xACD = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_3).*X_pcr_ij(:,k_4));
                one_xBCD = one_x * (X_pcr_ij(:,k_2).*X_pcr_ij(:,k_3).*X_pcr_ij(:,k_4));
                one_xABCD = one_x * (X_pcr_ij(:,j).*X_pcr_ij(:,k_2).*X_pcr_ij(:,k_3).*X_pcr_ij(:,k_4));
                
                test_x = one_xA * one_xB * one_xC * one_xD...
                    -   ( (one_xAB * one_xC * one_xD)...
                        + (one_xAC * one_xB * one_xD)...
                        + (one_xAD * one_xB * one_xC)...
                        + (one_xBC * one_xA * one_xD)...
                        + (one_xBD * one_xA * one_xC)...
                        + (one_xCD * one_xA * one_xB) )...
                    + 2*( (one_xABC * one_xD)...
                        + (one_xABD * one_xC)...
                        + (one_xACD * one_xB)...
                        + (one_xBCD * one_xA) )...
                    +   ( (one_xAB * one_xCD)...
                        + (one_xAC * one_xBD)...
                        + (one_xAD * one_xBC) )...
                    - 6*(  one_xABCD );
                if test_x > 0.0001
                    modk_2 = k_2 - j;
                    while modk_2 < 0
                        modk_2 = modk_2 + q;                        
                    end
                    modk_3 = k_3 - j;
                    while modk_3 < 0
                        modk_3 = modk_3 + q;
                    end
                    modk_4 = k_4 - j;
                    while modk_4 < 0
                        modk_4 = modk_4 + q;
                    end
                    x1_e_k_2k_3k_4(1 + modk_2 + q*modk_3 + q_squared*modk_4, 1) = test_x;
                    x1_e_k_2k_3k_4(1 + modk_2 + q*modk_4 + q_squared*modk_3, 1) = test_x;
                    x1_e_k_2k_3k_4(1 + modk_3 + q*modk_2 + q_squared*modk_4, 1) = test_x;
                    x1_e_k_2k_3k_4(1 + modk_3 + q*modk_4 + q_squared*modk_2, 1) = test_x;
                    x1_e_k_2k_3k_4(1 + modk_4 + q*modk_2 + q_squared*modk_3, 1) = test_x;
                    x1_e_k_2k_3k_4(1 + modk_4 + q*modk_3 + q_squared*modk_2, 1) = test_x;
                end
            end
        end
    end
    x_e_k_2k_3k_4 = x_e_k_2k_3k_4 + x1_e_k_2k_3k_4;
    x1_e_k_2k_3k_4 = sparse(q_cubed, 1);
end
clear x1_e_k_2k_3k_4

% test_x = squeeze(sum(reshape(full(x_e_k_2k_3k_4), [q, q, q]), 3));
% figure(8), subplot(1, 2, 1)
% surf(test_x, 'FaceAlpha',1,'LineStyle','none')
% axis([1 q+1 1 q+1])
% colormap bone
% lighting phong
% grid off
% set(gca,'PlotBoxAspectRatio',[1,1,1])
% set(gca,'XTick',1:q*100/period:q+1)
% set(gca,'YTick',1:q*100/period:q+1)
% set(gca,'XTickLabel',0:period/100)
% set(gca,'YTickLabel',0:period/100)
% clear test_x

% Expectation array Y
Y_pcr_j = sum(Y_pcr_ij);
nonzerosY_pcr_j = find(Y_pcr_j);
one_y = ones(1,d_y);

y_e_k_2k_3k_4 = sparse(q_cubed, 1);
y1_e_k_2k_3k_4 = sparse(q_cubed, 1);
for j = nonzerosY_pcr_j
    for k_2 = nonzerosY_pcr_j
        for k_3 = nonzerosY_pcr_j(nonzerosY_pcr_j >= k_2);
            for k_4 = nonzerosY_pcr_j(nonzerosY_pcr_j >= k_3);
                one_yA = Y_pcr_j(j);
                one_yB = Y_pcr_j(k_2);
                one_yC = Y_pcr_j(k_3);
                one_yD = Y_pcr_j(k_4);
                one_yAB = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_2));
                one_yAC = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_3));
                one_yAD = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_4));
                one_yBC = one_y * (Y_pcr_ij(:,k_2).*Y_pcr_ij(:,k_3));
                one_yBD = one_y * (Y_pcr_ij(:,k_2).*Y_pcr_ij(:,k_4));
                one_yCD = one_y * (Y_pcr_ij(:,k_3).*Y_pcr_ij(:,k_4));
                one_yABC = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_2).*Y_pcr_ij(:,k_3));
                one_yABD = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_2).*Y_pcr_ij(:,k_4));
                one_yACD = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_3).*Y_pcr_ij(:,k_4));
                one_yBCD = one_y * (Y_pcr_ij(:,k_2).*Y_pcr_ij(:,k_3).*Y_pcr_ij(:,k_4));
                one_yABCD = one_y * (Y_pcr_ij(:,j).*Y_pcr_ij(:,k_2).*Y_pcr_ij(:,k_3).*Y_pcr_ij(:,k_4));
                
                test_y = one_yA * one_yB * one_yC * one_yD...
                    -   ( (one_yAB * one_yC * one_yD)...
                        + (one_yAC * one_yB * one_yD)...
                        + (one_yAD * one_yB * one_yC)...
                        + (one_yBC * one_yA * one_yD)...
                        + (one_yBD * one_yA * one_yC)...
                        + (one_yCD * one_yA * one_yB) )...
                    + 2*( (one_yABC * one_yD)...
                        + (one_yABD * one_yC)...
                        + (one_yACD * one_yB)...
                        + (one_yBCD * one_yA) )...
                    +   ( (one_yAB * one_yCD)...
                        + (one_yAC * one_yBD)...
                        + (one_yAD * one_yBC) )...
                    - 6*(  one_yABCD );
                if test_y > 0.0001
                    modk_2 = k_2 - j;
                    while modk_2 < 0
                        modk_2 = modk_2 + q;
                    end
                    modk_3 = k_3 - j;
                    while modk_3 < 0
                        modk_3 = modk_3 + q;
                    end
                    modk_4 = k_4 - j;
                    while modk_4 < 0
                        modk_4 = modk_4 + q;
                    end
                    y1_e_k_2k_3k_4(1 + modk_2 + q*modk_3 + q_squared*modk_4, 1) = test_y;
                    y1_e_k_2k_3k_4(1 + modk_2 + q*modk_4 + q_squared*modk_3, 1) = test_y;
                    y1_e_k_2k_3k_4(1 + modk_3 + q*modk_2 + q_squared*modk_4, 1) = test_y;
                    y1_e_k_2k_3k_4(1 + modk_3 + q*modk_4 + q_squared*modk_2, 1) = test_y;
                    y1_e_k_2k_3k_4(1 + modk_4 + q*modk_2 + q_squared*modk_3, 1) = test_y;
                    y1_e_k_2k_3k_4(1 + modk_4 + q*modk_3 + q_squared*modk_2, 1) = test_y;
                end
            end
        end
    end
    y_e_k_2k_3k_4 = y_e_k_2k_3k_4 + y1_e_k_2k_3k_4;
    y1_e_k_2k_3k_4 = sparse(q_cubed, 1);
end
clear y1_e_k_2k_3k_4

% test_y = squeeze(sum(reshape(full(y_e_k_2k_3k_4), [q, q, q]), 3));
% figure(8), subplot(1, 2, 2)
% surf(test_y, 'FaceAlpha',1,'LineStyle','none')
% axis([1 q+1 1 q+1])
% colormap bone
% lighting phong
% grid off
% set(gca,'PlotBoxAspectRatio',[1,1,1])
% set(gca,'XTick',1:q*100/period:q+1)
% set(gca,'YTick',1:q*100/period:q+1)
% set(gca,'XTickLabel',0:period/100)
% set(gca,'YTickLabel',0:period/100)
% clear test_y

% Distance calculation
p_norm_x_e_k_2k_3k_4 = sum(x_e_k_2k_3k_4.^p)^(1/p);
p_norm_y_e_k_2k_3k_4 = sum(y_e_k_2k_3k_4.^p)^(1/p);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isequal(metric, 'l_p')
    d = sum(abs(x_e_k_2k_3k_4./(2*p_norm_x_e_k_2k_3k_4)...
        - y_e_k_2k_3k_4./(2*p_norm_y_e_k_2k_3k_4)).^p)^(1/p);
elseif isequal(metric, 'cosine')
    d = 1 - ((x_e_k_2k_3k_4' * y_e_k_2k_3k_4)/...
        sqrt((x_e_k_2k_3k_4' * x_e_k_2k_3k_4)...
        * (y_e_k_2k_3k_4' * y_e_k_2k_3k_4)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%clear x_e_k_2k_3k_4 y_e_k_2k_3k_4 

else
error('Function does not currently exist - use a lower value for "r".');
end
end
