NIRS_ICA_V1.0/ica_algorithms/tentative/STICA_skew_demo/STICA_skew_demo/STICA_CORE/basicfg.m
function [f, g, y] = basicfg(W,X,ind_hi,ind_lo)

% Basic Bell-Sejnowski gradient

%fprintf('Basicfg ');
% ind_hi  ind_lo

[n M]= size(X);
y = X*W;
[n K]=size(y);
J =  abs(det(W));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIND f.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

REORDER=1;
if REORDER
	ind = [ind_hi ind_lo];
	% Find inverse sorting indices.
	[temp ind_inv] = sort(ind);
	ya = y(:,ind_hi);
	yb = y(:,ind_lo); 
end;

p = sech(ya).^2; % High kurt
a = log(p + eps);

b = -yb.^4; % Low kurt, -yb^4 = log (exp(-yb^4)), it is exp(-yb^4) that has low k.
ab=[a b];

f = log(J) + sum(mean(ab));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIND grad.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Get deriv and take account of eps in function and deriv (p+eps).
ga = 2/n * tanh(ya) .* (p./(p+eps));
gb = 4/n * yb.^3;

gab = [ga gb];

if REORDER 
	gab = reorder_cols(gab,ind_inv); 
end;

g = inv(W') - X' * gab;

% old code .......
% f = log(abs(det(W)))+sum(mean(log(sech(Y).^2)));
% g = inv(W')-2/n*X'*tanh(Y);

