function [phi,psi] = mnf(X)
% function [phi,psi] = mnf(X)
% Extracts the maximum noise fraction components of the data in X.
% and projects points onto these components
%
% In:
%    X - an NxP data matrix  (P is number of samples)
%
% Out:
%    phi - the basis vectors (NxP)
%    psi - the mixing coefficients (NxN)
%
%
% Error: This code will return null values for the two output parameters
%        if there is an error during execution.
%
% This follows Hundley, Kirby and Anderle (not anymore)
%

X = X';

[numSamples dim] = size(X);

%mean center the columns
%X = X - repmat(mean(X),numSamples,1);

try
	[Vx,Lx] = sorteig(X'*X);
catch
	phi = [];
	psi = [];
	fprintf('An error occured while computing the eigen decomposition of the covariance matrix.\n')
	lasterr
	return;
end

DSx = sqrt(diag(Lx));
INVDSx = diag(1./(DSx));

%estimate covariance matrix of Z = Q'*Q
SX(1,:) = X(2,:);
SX(2:size(X,1),:) = X(1:size(X,1)-1,:);
Z1 = SX'*X;
s1 = norm(X-SX,'fro')/norm(X,'fro');
Z2 = X'*SX;
s2=norm(Z2);
Z3 = X'*X;
Z4 = SX'*SX;
Z = 0.5*(Z3+Z4-Z1-Z2);
ZHAT = INVDSx*Vx'*Z*Vx*INVDSx;

try
	[PSIHAT,DPSIHAT] = sorteig(ZHAT);
catch
	phi = [];
	psi = [];
	fprintf('An error occured while computing the eigen decomposition of the ZHAT.\n')
	fprintf('%s\n\n',lasterr)
	return;
end

psi = Vx*INVDSx*PSIHAT;
phi = X*psi;

phi = phi';
psi = psi';

% this will make all signals have positive activation
%phi(find(sum(psi)<0),:) = -phi(find(sum(psi)<0),:);

