function y = diwish(W,v,S)
% density function of Inverse Wishart distribution
% based on diwish function of MCMCpack in R
% v Degrees of freedom, scale
% S Scale matrix (p \times p)
% W Positive definite matrix W (pxp).
% We assmue that S and W are square matrices and have the same dimension
% We also assume v>dim(S)
% The mean of an inverse Wishart random variable with v degrees of freedom 
% and scale matrix S is S/(v-p-1)

% p <- nrow(S)
% gammapart <- sum(lgamma((v + 1 - 1:p)/2))
% ldenom <- gammapart + 0.5 * v * p * log(2) + 0.25 * p * (p - 1) * log(pi)
% cholS <- chol(S)
% cholW <- chol(W)
% halflogdetS <- sum(log(diag(cholS)))
% halflogdetW <- sum(log(diag(cholW)))
% invW <- chol2inv(cholW)
% exptrace <- sum(S * invW)
% lnum <- v * halflogdetS - (v + p + 1) * halflogdetW - 0.5 * exptrace
% lpdf <- lnum - ldenom
% return(exp(lpdf))

p = size(S,1);
gammapart = sum(log(gamma((v+1-(1:p))/2)));
ldenom = gammapart + 0.5 * v * p * log(2) + 0.25 * p * (p - 1) * log(pi);
cholS = chol(S, 'upper');
cholW = chol(W, 'upper');
halflogdetS = sum(log(diag(cholS)));
halflogdetW = sum(log(diag(cholW)));
invW = inv(cholW'*cholW);
exptrace = sum(S .* invW, 'all');
lnum = v * halflogdetS - (v + p + 1) * halflogdetW - 0.5 * exptrace;
lpdf = lnum - ldenom;
y = exp(lpdf);



end


