function distance = AMD(spikes)
%--------------------------------------------------------------------------
% AMD.m - Calculates the AMD distance between all possible pairs of a set
% of spiketrains.  Spikes is a matrix where each row is a spiketrain.
%
% Usage: distance = AMD(spikes)
%
% Input:  spikes             * MxN matrix of M spiketrains
% Output: distance           * MxM matrix of pairwise spiketrain distances
%
% Written by Marshall Crumiller
%--------------------------------------------------------------------------
[M,N] = size(spikes);

distance = zeros(M); % record our AMD distances

% For each spiketrain, we're going to duplicate it into one single long
% spiketrain, and string all the possible pairs together to form one long
% vector.  This way, we can perform the entire operation on a single
% vector.
for i = 1:M
    pair_index = 1:M; pair_index(i) = [];
    concat_len = (M-1)*N;
    num_spikes = length(spikes(i,spikes(i,:)~=0));
    current_spiketrain = repmat(spikes(i,:),1,M-1);
    compare_spiketrain = reshape(spikes(pair_index,:)',1,(M-1)*N);
    spike_diff = abs(current_spiketrain-compare_spiketrain);
    spike_diff_pos = inf(1,concat_len);
    spike_diff_neg = inf(1,concat_len);
    change = inf(1,concat_len);
    shift = 1;
       
    % perform shift operation for (i,j)
    while(~all(~(change-spike_diff)))
        change = spike_diff;
        % right shift
        spike_diff_pos(shift+1:end) = abs(current_spiketrain(shift+1:end) - compare_spiketrain(1:end-shift));
        % left shift
        spike_diff_neg(1:end-shift) = abs(current_spiketrain(1:end-shift) - compare_spiketrain(shift+1:end));
        spike_diff = min(min(spike_diff,spike_diff_pos), spike_diff_neg);
         shift=shift+1;
    end
    
    % grab sums
    spike_diff = reshape(spike_diff,M-1,N);
    sums = sum(spike_diff,2);
    sums = sums./num_spikes;
    
    distance(i,pair_index) = sums;
    % put sums into distance vector
end

distance = distance+distance';