NIRS_ICA_V1.0/ica_algorithms/tentative/STICA_skew_demo/STICA_skew_demo/STICA_CORE/stica.m
function [V, d, S, T, w] = stica(P, Q, alpha, V0, hi_kurt_s,lo_kurt_s, hi_kurt_t, lo_kurt_t,C1,C2)

% function [V, S, T] = stica(P, Q, alpha, a, b)
%
% spatio-temporal ICA 
% based on vanilla B & S
%
% P, Q = m x k spatio-temporal data to unmix
% P = Spatial, Q=temporal - data in columns.
% W (V0?)   = (k+1) x k spatial unmixing matrix
% d    = temporal scaling k-vector
% S, T = unmixed signals
% specify hi/lo kurtosis signals.
% V = spatial unmix matrix.
% C1 and C2 are cov matrices of pre and post data 
% or are defined by F=log(V/U) = log(wC1w)/(wC2w)function.

global fundat mindat;
global plot_inverval;
global global_data;

k = size(P, 2); % number of components

% data structure used by sticafg
fundat.P = P;  % One vector per col, P = spatial.
fundat.Q = Q; 
fundat.alpha = alpha;
fundat.k = k;
fundat.nC = 1;
fundat.C1 = C1;
fundat.C2 = C2;

%bias = zeros(1,k); % row vector 1xk.
%fundat.Pdf_bias=bias; % biases to ensure pdf of data has zero mean.
% Set skew params.
%fundat.a = 4; % 4 Large => steep curve on lhs.
%fundat.b=1; % 1	Large => steep curve on rhs.

fundat.a = 4*ones(1, k); % 4 	Large => steep curve on lhs.
fundat.b = 1*ones(1, k); % 1	Large => steep curve on rhs.

% Specify hi/lo kurtosis.
fundat.hi_kurt_s = hi_kurt_s;
fundat.lo_kurt_s = lo_kurt_s;
fundat.hi_kurt_t = hi_kurt_t;
fundat.lo_kurt_t = lo_kurt_t;
% SET GLOBAL DATA
global_data = fundat;

if nargin < 4
    V0 = eye(k);
    d0 = ones(k, 1);
else
    T0 = Q*inv(V0');
    d0 = 0.9069./std(T0);
end;

b0=zeros(k,1);

jsize(V0,'stica: V0');
jsize(d0,'stica: d0');

w0 = [V0(:); d0(:); b0(:); b0(:)]; % [weights, scaling, bias_s bias_t].
    
% Gauss-newton mnimisation
mindat = []; % default parameters for G-N
mindat.grat = 0.1;
mindat.nmax = 500;
mindat.gtol = 1e-3; %1e-3;
mindat.w0=w0;
mindat.w=w0;
plt=1; %% MOD

% GAUSSNEWT
%w = gaussnewt(w0, 'sticafg', mindat, fundat);

GRADTEST=0;
if GRADTEST
	% GRAD TEST.
	fprintf('Gradtest ...\n');
	% STICAFG - GOOD.
	gradtest('sticafg', w0, fundat, sqrt(eps));
	% basicfg_skew = GOOD.
	% ww = [V0(:);  b0(:)]+eps; gradtest('basicfg_skew', ww, fundat, sqrt(eps));
	return;
end;

% CONJ GRAD
w = conjgrad_jvs(w0,'sticaf','sticag', fundat, plot_inverval, mindat.nmax, mindat.gtol);

% NO GRAD.
%options=foptions;
%w=fmins('sticaf',w0,options,[],fundat);

% FMINU
% unconstrained minimisation over w
options=foptions;
tol = 1e-6;
options(2)=tol;
options(3)=tol;
options = foptions;
options(14) = 1e6;	% max number of calls to heval
% options(9) = 1;	% check gradient calculation
%w = jfminu('sticaf', w0, options, 'sticag', fundat);

% FMINUNC
% unconstrained minimisation over w
%w = fminunc(sticaf, w0, options, sticag, fundat);


% recover optimal unmixing from weght vector 
V = reshape(w(1:k*k), k, k);
d = w(k*k+1:k);

% unmixed data (no temporal weighting so X = S*T')
S = P*V;
T = Q*inv(V');

fprintf('... stica done.\n');
