function theLastPairwiseDistanceComputationComparisonYoullEverNeed()
Funcs =  {@matheburgMatrixMultiplication, ...
          @DivakarVectorizedVariation1, ...
          @DivakarVectorizedVariation2, ...
          @pdist2, ...
          @statinferbsxfun, ...
          @Jonasndgrid};
    
ns = 2.^(1:11);
dims = [2,3;50,200].';
for i = 1:numel(dims) 
    d = dims(i);
    subplot(size(dims,2), size(dims,1), i);
    paramGenerators = arrayfun(@(n) @(){rand(n,d), rand(n,d)}, ns,'uni',0);
    times = compareFunctions(Funcs, paramGenerators);
    loglog(ns, times.');
    legend(cellfun(@func2str,Funcs,'uni',0),'Location','NorthWest','Interpreter','none');
    title(sprintf('DIM = %d',d));
    drawnow;
end
end

function times = compareFunctions(Funcs, paramGenerators)
times = [];
for i = 1:numel(paramGenerators)
    for j = 1:numel(Funcs)
        Params = feval(paramGenerators{i});
        times(end+1) = timeit(@()feval(Funcs{j},Params{:}));
    end
end
times = reshape(times,length(Funcs),[]);
end
%% // #################################
%% // HERE COME ALL THE FANCY FUNCTIONS
%% // #################################
function distMat = matheburgMatrixMultiplication(A,B)
numA = size(A,1);
numB = size(B,1);
d = size(A,2);
helpA = zeros(numA,3*d);
helpB = zeros(numB,3*d);
for idx = 1:d
    helpA(:,3*idx-2:3*idx) = [ones(numA,1), -2*A(:,idx), A(:,idx).^2 ];
    helpB(:,3*idx-2:3*idx) = [B(:,idx).^2 ,    B(:,idx), ones(numB,1)];
end
distMat = sqrt(helpA * helpB');
end
%%
function distMat = DivakarVectorizedVariation1(A,B)
[nA,dim] = size(A);
nB = size(B,1);

A_ext = ones(nA,dim*3);
A_ext(:,2:3:end) = -2*A;
A_ext(:,3:3:end) = A.^2;

B_ext = ones(nB,dim*3);
B_ext(:,1:3:end) = B.^2;
B_ext(:,2:3:end) = B;

distMat = sqrt(A_ext * B_ext.');
end
%%
function distMat = DivakarVectorizedVariation2(A,B)
[nA,dim] = size(A);
nB = size(B,1);

A_ext = [ones(nA*dim,1) -2*A(:) A(:).^2];
B_ext = [B(:).^2 B(:) ones(nB*dim,1)];

A_ext = reshape(permute(reshape(A_ext,nA,dim,[]),[1 3 2]),nA,[]);
B_ext = reshape(permute(reshape(B_ext,nB,dim,[]),[1 3 2]),nB,[]);

distMat = sqrt(A_ext * B_ext.');
end
%%
function distMat = statinferbsxfun(A,B)
distMat = sqrt(bsxfun(@plus,dot(A,A,2),dot(B,B,2)')-2*(A*B'));
end
%%
function distMat = Jonasndgrid(A,B)
n = size(A,1);
m = size(B,1);
[idxA,idxB] = ndgrid(1:n,1:m);
distMat = zeros(n,m);
distMat(:) = sqrt(sum((A(idxA,:) - B(idxB,:)).^2,2));
end