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

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

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

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

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