function val = jsdiv(p, q)
% jsdiv returns the js divergence of two discrete probability distributions

% Normalize
p = vec(p ./ sum(p));
q = vec(q ./ sum(q));

% Compute middle distribution
m = (p + q) ./ 2;

% Sum to get final answer
val = (kldiv2(p, m) + kldiv2(q, m)) ./ 2;

end

function val = kldiv2(p, q)
% kldiv2 returns the kl divergence of two discrete probability
% distributions using log2

% Normalize
p = vec(p ./ sum(p));
q = vec(q ./ sum(q));

% Compute value inside sum
rat = log2(p ./ q) .* p;

% Assume 0*log(0) = 0
rat(p == 0) = 0;

% Sum to get final answer
val = sum(rat);

end
