function [l,con] = lqr_eval(x,means,cholesky_covar,powr,idiag, temp);
% function l = lqr_eval(x,means,cholesky_covar,[powr],[idiag],[temp]);
% evaluates the log Gaussian PDF 
%
% Input data is dimensioned (DIM-by-N)
% means are dimensioned DIM-by-1 or DIM-by-N

%
% NAVAL UNDERSEA WARFARE CENTER
% class-specific@npt.nuwc.navy.mil
% http://www.npt.nuwc.navy.mil/csf
% Paul M. Baggenstoss
% 401-832-8240
% Do not disseminate.
%



%  Dr. Paul M. Baggenstoss
%  Naval Undersea Warfare Center
%  Newport, RI
%  p.m.baggenstoss@ieee.org
%  Date              Reason
%  ------------      --------------
%  Mar 28, 1997      Initial release
%  Dec 7, 1999       Did math a better way
%  Dr. Paul M Baggenstoss
%  NUWC, Newport, RI
%  p.m.baggenstoss@ieee.org
%
%  Aug 6, 2004        Added powr argument
%  Jan 9 2008         Added second dimension on means
%  Jan 10, 2012       Added idiag, temp argument
                                                                        
if(nargin < 4 || isempty(powr) ), powr=2; end;
if(nargin < 5 || isempty(idiag) ), idiag=0; end;
if(nargin < 6 || isempty(temp) ), temp=0; end;

[DIM,N]=size(x);
DIM0=size(means,1);
if(DIM0 ~= DIM),  
     error('lqr_eval: data and means have different dimension'); 
end;

if(size(means,2)==1),
  means = repmat(means(:),1,N);
else,
   if(size(means,2) ~= N), error('wrong dimensions for means'); end;
end;


logdet = 2 * sum(log(abs(diag(cholesky_covar)))) ;
tmpidx = x - means;
if(idiag),
    cc=diag(cholesky_covar);
    tmpidx = sum((tmpidx ./ repmat(cc,1,N)).^2,1);
else,
   tmpidx = sum(( cholesky_covar' \ tmpidx) .^ 2,1);
end;

if(powr==2),
    con = (DIM/2)*log(2*pi) + 0.5*logdet;
    l = -.5*tmpidx' - con;
else,
     %con = -log(powr) + (DIM/powr+1)*log(2) + (DIM/2)*log(pi) ...
     %        - gammaln(DIM/2.0) + gammaln(DIM/powr) + .5*logdet;
     con=DIM*(log(2/powr)+gammaln(1/powr)-log(2)/powr)+ .5*logdet;
         
     l = -.5*(tmpidx').^(powr/2) - con;
end;
if(temp), l=l/exp(temp); end;

