function [lg,idx,lcon]=lqr_evp(parm,x,flag,powr, idiag, temp)
%function lg=lqr_evp(parm,x,flag,[powr=2], [idiag=0],[temp=0])
% Returns log-likelihood of a Gaussian mixture.
% Flag=1,
% if flag==1, returns log mode PDF in columns, that is the
% output of lqr_eval only, not including mixing weights
% for each column.
% Flag=0,
% Computes the total log-PDF in a single
% output column (i.e. the log of the weighted sum of the rows of
% the flag=1 output)
% idx is the most likely mixture component

%
% 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
% Jul 19, 1997 Accounted for data scaling
% Oct 20, 1998 Fixed bug occuring when nmode=1
% Dec 6, 1999 MATLAB 5 Upgrade
% Aug 6, 2004 Added powr argument
% Feb 2008 added idx output, sped up code
% June 2010 Fixed comments. LQR_EVP now obsolete use eval_gmix.m
% Dec 19, 2010 expanded powr to each mode
% Jan 10, 2012 Added idiag, temp argument
% Dec 18, 2017 Added lcon output argument

nmode = length(parm.modes);

if (nargin < 4 || isempty(powr) ),
if(isfield(parm.modes(1),'powr')),
powr=[parm.modes.powr];
else,
powr=ones(1,nmode)*2;
end;
else,
if(length(powr)==1),
powr=ones(1,nmode)*powr;
end;
end;

if (nargin < 5 || isempty(idiag) ), idiag=0; end;
if (nargin < 6 || isempty(temp) ), temp=0; end;

wts = [parm.modes.weight]';
DIM0 = length(parm.features);
[DIM,N] = size(x);

if( DIM ~= DIM0) ,
error('LQREVAL_P: data and parms have different DIM');
end;

lg = zeros(N,nmode);
lcon=zeros(1,nmode);
for k=1:nmode,
[lg(:,k),lcon(k)] = lqr_eval(x,parm.modes(k).mean,parm.modes(k).cholesky_covar, powr(k), idiag, temp);
end;

[mx,idx] = max(lg,[],2);
if(flag==1), return; end;
if(nmode > 1),
lg=lg-repmat(mx,1,nmode);

% compute the exponent
g = exp(lg);

% Compute combined Gaussian sum function at each sample.
y = g * wts;

% add back in the scaling
lg = log(y) + mx;
end;