%NORMAL_MAP Map a dataset on normal-density classifiers or mappings
% 
%   F = NORMAL_MAP(A,W)
%
% INPUT
%   A   Dataset
%   W   Mapping
%
% OUTPUT
% 	F   Density estimation for classes in A	
%
% DESCRIPTION
% Maps the dataset A by the normal density based classifier or mapping W.
% For each object in A, F returns the densities for each of the classes or
% distributions stored in W. For classifiers, the densities are weighted
% by the class prior probabilities. This routine is automatically called for
% computing A*W if W is a normal density based classifier or a mapping.
%
% Use W = LOGDENS(W) (or W = W*LOGDENS) if absolute densities are not
% needed and a more accurate posterior probability is desired.
%
% SEE ALSO
% MAPPINGS, DATASETS, QDC, UDC, LDC, GAUSSM, LOGDENS

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: normal_map.m,v 1.15 2005/10/18 14:07:36 duin Exp $

function F = normal_map_new(A,W)

	prtrace(mfilename);

	w = +W;          % data field of W (fields: w.mean, w.cov, w.prior, w.nlab)
	                 % each of these data fields has a mean, cov, prior for separate
									 % Gaussian components. The nlab field assigns each component
									 % to a class. 
	[k,c] = size(W); % c is number of classes
	% DEG = 1 indicates a common cov. matrix and DEG = 2 - separate cov. matrices.
	deg = ndims(w.cov)-1; 		
	U = w.mean; G = w.cov; p = w.prior;
	if (abs(1-sum(p)) > 1e-6)
		error('Class or component probabilities do not sum to one.')
	end
	lablist = getlab(W);

	[m,ka] = size(A);
	if (ka ~= k), 
		error('Feature sizes of the dataset and the mapping do not match.'); 
	end

	n = size(U,1);			% Number of components.
	F = zeros(m,n); 		% Gaussian densities for each component to be computed.

	if (deg == 1)		
		H = G;
		if (rank(H) < size(H,1))
      prwarning(2,'Singular case, pseudo-inverse of the covariance matrix is used.');
			E = real(pinv(H));
		else
			E = real(inv(H));
		end
    end

    %detect circular features
    cFeaturesDomain = getfeatdom(A);

    DirectionalFeatures = [];

    for ii = 1:length(cFeaturesDomain)
        if(~isempty(cFeaturesDomain{ii}))
            DirectionalFeatures = [DirectionalFeatures ii];
        end
    end

    iCantFeatures = ka;
    
    aa_matrix = +A;
    NotDirectionalFeatures = setdiff(1:iCantFeatures, DirectionalFeatures);
    
    bNotNaN = ~any(isnan( aa_matrix ),2);
    iCantElementos = sum(bNotNaN);

	% Loop over components.
	for i=1:n

		% Shift A such that the mean lies at the origin.
        if( ~isempty(NotDirectionalFeatures) )
            X(bNotNaN,NotDirectionalFeatures) = aa_matrix(bNotNaN,NotDirectionalFeatures) - ones(iCantElementos,1)*U(i,NotDirectionalFeatures);
        end

        if( ~isempty(DirectionalFeatures) )
            for ii = DirectionalFeatures                                     
                iAux1 = aa_matrix(bNotNaN, ii );
                iAux2 = repmat(U(i,ii), iCantElementos,1);
                iAux3 = [ iAux1 - iAux2, 2*pi + iAux1 - iAux2, iAux1 - 2*pi - iAux2 ];
                [dummy, iMinIndex] = min(abs(iAux3), [], 2);
                iMinIndex = sub2ind(size(iAux3),1:size(iAux3,1),iMinIndex');
                X(bNotNaN,ii) = iAux3(iMinIndex);
            end 
        end
		
        if (deg == 2)		
			H = G(:,:,i);
			if (rank(H) < size(H,1))
				prwarning(1,'Singular case, pseudo-inverse of the covariance matrix is used.');
				E = real(pinv(H));
			else
				E = real(inv(H));
			end
		end
		% Gaussian distribution for the i-th component. Take log of density to preserve tails
		F(:,i) = -0.5*sum(X'.*(E*X'),1)' - (sum(log(real(eig(H)+1e-16)+realmin)) + k*log(2*pi))*0.5;
		if (getout_conv(W) ~= 2) % take log of density to preserve tails
			F(:,i) = exp(F(:,i));
		end
	end

	if isfield(w,'nlab') 
		% For mixtures of Gaussians. Relates components to classes
		nlab = w.nlab;
		cc = max(w.nlab);
	else
		nlab = [1:c];
		cc = c;
	end

	if (getout_conv(W) == 2)
		Cmax = max(F(:));  % scale to gain accuracy in the tails
		F = F - Cmax;
	end
	FF = zeros(m,cc);
	% Loop over true classes. Weight the probabilities by the priors.
	for j = 1:cc
		J = find(nlab == j); 
		if (getout_conv(W) == 2) % take log of density to preserve tails
			% difficult to get this right for MOGs, and moreover, probably not of
			% much help. Anyway, we give it a try.
			FF(:,j) = exp(F(:,J))*w.prior(J)';
			% like in parzen_map, use in tails just largest component
			L = find(FF(:,j) <= 1e-300); 
			N = find(FF(:,j) >  1e-300); 
			if ~isempty(L)
				[FM,R] = max(F(L,J),[],2);
				FF(L,j) = FM + w.prior(R)';
			end
			if ~isempty(N)
				FF(N,j) = log(FF(N,j));
			end
		else
			FF(:,j) = F(:,J) * w.prior(J)';
		end
	end
	
	if (getout_conv(W) == 2)
		FF = exp(FF);
	else
		FF = FF + realmin; % avoid devision by 0 in computing posterios later
	end

	if isdataset(A)
		F = setdata(A,FF,lablist);
	else
		F = FF;
	end

return;