function H = entropy2(pin,varargin);
% H = entropy2(p,q) = -sum_{i:q(i)>0} q(i)*log2(p(i))
% If q=p, use H = entropy2(p).
% p and q can be matrices; e.g.,
%
% -sum_{i,j} q(i,j)*log p(j|i)
%
% where q(i,j) may be the joint pmf.
%
% You can also compute the relative entropy or
% informational divergence D(q||p) = entropy2(p./q,q)
if nargin > 1
    qin = varargin{1};
else
    qin = pin;
end
lt = prod(size(pin));
p = reshape(pin,lt,1); % convert pin to column vector
q = reshape(qin,1,lt); % convert qin to row    vector
i = find(q>0);
H = - q(i) * log2( p(i) );