NIRS_ICA_V1.0/ica_algorithms/tentative/STICA_skew_demo/STICA_skew_demo/STICA_CORE/sticafg.m
function [f, g] = sticafg(w, data, plt)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [f, g] = sticafg(w, data)
% 
% Evaluate cost and gradient for spatio-temporal ICA
% based on vanilla B & S
%
% w     = unmixing matrix and scaling as col (k*k+k)-vector
% w contains 2 extra rows:
%	row K+1:	Scaling parameters.
%	row k+2:	Skew offset parameters or bias.
% data  = structure with 3-fields:
%         P, Q  = m x k spatial & temporal data sets
%         alpha = spatio-temporal bias (alpha = 1: all spatial)
%
% f = output entropy (tanh sigmoids)
% g = analytic gradient of f
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global SPATIAL_ICA TEMPORAL_ICA SPATIOTEMPORAL_ICA;
global num_funevals plot_inverval;
global ys yt;
global beta;
global wm_sm;
global SKEW_PDF_s SKEW_PDF_t;
global mindat;

if nargin<3 plt=1; end;

num_funevals = num_funevals+1;
override=0;
% if rem(num_funevals,50)==0
% 	override=1;
% else
% 	override=0;
% end;

plt=1;
% jsize(w,'w');jsize(data,'data');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get spatial and temporal weight matrices.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% V=spatial (kxk), W=TEMPORAL (kxk), D=scalings.
% [V W D]=get_VWD(w,data);
%[V W D B]=get_VWDB(w,data); % V=spatial (kxk), W=TEMPORAL (kxk), D=scalings, B=BIAS.

[V W D Bs Bt]=get_VWDBsBt(w,data); % V=spatial (kxk), W=TEMPORAL (kxk), D=scalings, B=BIAS.

% Keep record of params in mindat.
mindat.w=w;
mindat.V=V;
mindat.W=W;
mindat.D=D;
mindat.Bs=Bs;
mindat.Bt=Bt;

[m, k] = size(data.P); % m=num samples, k=num mixtures - one mixture per col.
% jsize(w,'w');jsize(data,'data');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spatial entropy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % NOTE: ALWYAS NEED 	TO GET ys.
%if SPATIOTEMPORAL_ICA | SPATIAL_ICA
if SKEW_PDF_s & SPATIAL_ICA
	% Augment square V wtih bias row.
	VBs = [V; Bs]; % Bs=row.
	[fs, gs, ys] = basicfg_skew(VBs,data.P,data);
	% gs is k*(k+1), the last k elements are bias grads.
	% Put gbias at end row of gs.
	%jsize(gs,'gs');
	gs=reshape(gs,k+1,k);
	% Get last row of gs for bias_grad.
	bias_grad_s = gs(k+1,:)';
	gs=gs(1:k,1:k);
	w_grad_s=gs(:);
	% jsize(V,'V');jsize(ys,'ys');pr;
else % SPATIAL
	[fs, gs, ys] = basicfg(V,data.P,data.hi_kurt_s,data.lo_kurt_s);
	% SET DERIV WRT SCALE FACTORS TO ZERO FOR SPATIAL.
	bias_grad_s =  zeros(k, 1);	
	w_grad_s=gs(:);
	% jsize(V,'V');jsize(ys,'ys');pr;
end;

% COLLATE PARTS OF SPATIAL GRAD.
scale_grad_s =  zeros(k, 1); % ALWAYS ZERO.
bias_grad_t = zeros(k, 1);
gs = [w_grad_s; scale_grad_s;  bias_grad_s; bias_grad_t];
%if plt	hplos(ys,override); end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if plt	hplos(ys); end;
% if plt	hplot(yt); end;
%end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Temporal entropy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if SKEW_PDF_t & TEMPORAL_ICA
	% Augment square W wtih bias row.
	WBt = [W; Bt]; % B=row.
	[ft, gt, yt] = basicfg_skew(WBt,data.Q,data);
	% gs is k*(k+1), the last k elements are bias grads.
	% Put gbias at end row of gs.
	gt=reshape(gt,k+1,k);
	% Get last row of gs for bias_grad.
	bias_grad_t = gt(k+1,:)';
	gt=gt(1:k,1:k);
else  % TEMPORAL 
	[ft, gt, yt] = basicfg(W,data.Q,data.hi_kurt_t,data.lo_kurt_t);
end;
gtWiD = gt'*W*inv(D);
gt1 = -W*gtWiD; 
w_grad_t = gt1(:);
% if plt	hplot(yt,override); end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GET FINAL TEMPORAL GRADIENT VECTOR WRT SCALING AND BIAS.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if SPATIOTEMPORAL_ICA
	gt2 = diag(gtWiD); % Gradient wrt scaling parameters.
	scale_grad_t = gt2(:);
elseif TEMPORAL_ICA 
	% SET DERIV WRT SCALE FACTORS TO ZERO.
	scale_grad_t =  zeros(k, 1);
elseif SPATIAL_ICA
	% SET DERIV TO ZERO.
	gt1=gt1*0;
	scale_grad_t =  zeros(k, 1);
	bias_grad_t = zeros(k, 1);
end;
% Deriv wrt bias_grad_s is always zero in gt.
bias_grad_s = zeros(k, 1);	
gt = [w_grad_t; scale_grad_t; bias_grad_s; bias_grad_t];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Combined entropy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if SPATIOTEMPORAL_ICA
	
	alpha = data.alpha; 

	f = alpha*fs+(1-alpha)*ft; 

	[gs_core, gscale, gsbias, gtbias] = get_grad_parts(gs,k);
	gs_extra = [gscale; gsbias; gtbias];
	
	[gt_core, gscale, gsbias, gtbias] = get_grad_parts(gt,k);
	gt_extra = [gscale; gsbias; gtbias];

	% Only mix core grads.
	g = alpha*gs_core+(1-alpha)*gt_core;

	ch=0;
	if ch
	% Add grad of extra bits (should OR to filled vector).
	inds=(gs_extra~=0);
	indt=(gt_extra~=0);
	ist=inds+indt;

	jprint_vec(gs,1,'gs');
	jprint_vec(gs_core,1,'gs_core');
	jprint_vec(gs_extra,1,'gs_extra');
	jprint_vec(gt_extra,1,'gt_extra');
	jprint_vec(ist,1,'ist');

	fprintf('Should be all and only ones in ist\n');
	pr;	
	end;

	g_extra=gs_extra+gt_extra;
	g=[g; g_extra];

elseif SPATIAL_ICA
 	f = fs; g = gs;
elseif TEMPORAL_ICA
 	f = ft; g = gt;
else 
	error('sticafg');
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% USE WEAK MODEL?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if wm_sm % global var.
	% W = temporal, V=spatial.
	% beta = 0.5;
	% Get w=first col of temporal wt matrix W.
	% w = W(:,1);
	v=V(:);
	[wmf wmg]=wmfg(v,data);
	f = (1-beta)*f+beta*wmf;
	% jsize(g,'g');
	% jsize(wmg,'wmg');
	gt2 = g(k*k+1:k*k+k);
	g = (1-beta)*g(1:k*k) + beta*wmg;
	g=[g; gt2];
	% jsize(g,'g'); 20 x 1
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% want to minimise entropy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = -f;
g = -g;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

