NIRS_ICA_V1.0/ica_algorithms/tentative/STICA_skew_demo/STICA_skew_demo/STICA_CORE/stica_demo.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stica_demo;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all; clear all;


global num_funevals plot_inverval num_linesearches;
global ncall hmin wmin pcs pct sval nr nc;
global U D V nt nr nc;
global nplt P Q hs ncall nic neig;
global SPATIAL_ICA TEMPORAL_ICA SPATIOTEMPORAL_ICA;
global fundat mindat temp_s temp_t mode;
global U0 V0 D0 U_orig V_orig D_orig beta;
global wm_sm;
global SKEW_PDF_s SKEW_PDF_t;

tic;

%----------------------------------------------------------------------
% SET GLOBALS
%----------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERAL GLOBALS.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RAND_SEED 	= 111; % 999 for original REAL DATA for neuroimage paper
jrand_set_seed(RAND_SEED); % 999

func_count			= 0;
num_linesearches 	= 0;
temp_s				= 0;
temp_t				= 0;
num_funevals		= 0; 
plot_inverval		= 10; % counts line searches.

init		= 1; % New everything.
init_W		= 1; % New wt matrix.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GLOBALS FOR ICA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
neig 		= 4; % 16 	% number of eigenvectors used for ICA.
mode		= 'st'; 	% s=spatial, t=temporal, st=spatiotemporal
wm_sm		= 0; 		% Use weak model (temporal smoothness).
beta 		= 0.5; 		% 0.5; weak model - see sticafg.m.
SKEW_PDF_s	= 1;        % use skew pdf for spatial model (if set, ignore hi_kurt_s & lo_kurt_s)
SKEW_PDF_t	= 0;        % use skew pdf for temporal model (if set, ignore hi_kurt_t & lo_kurt_t)

% SET KURTOSIS OF REQUIRED SPATIAL AND TEMPORAL ICS.
% CLASSIC hi_kurt_s=[1 2];lo_kurt_s=[3 4];	hi_kurt_t = [1 3];lo_kurt_t = [2 4];
% GOOD hi_kurt_s=[1 2]; lo_kurt_s=[3 4];		hi_kurt_t=[1 3 4]; lo_kurt_t=[2];
hi_kurt_s=1:neig;   
lo_kurt_s=[];		

hi_kurt_t=1:neig; 
lo_kurt_t=[];

hs	= [];
ncall 	= 0;	% counts calls to heval
hmin 	= 1e10;

nic 	= neig; 	% number of ICs to extract.
nplt	= 10;

% SET FLAGS.
SPATIAL_ICA		= 0; 
TEMPORAL_ICA		= 0;
SPATIOTEMPORAL_ICA	= 0;
if  strcmp(mode,'st')
	SPATIOTEMPORAL_ICA	= 1;
	SPATIAL_ICA			= 1; 
	TEMPORAL_ICA		= 1;
	fprintf('Using SPATIOTEMPORAL_ICA\n\n');
elseif strcmp(mode,'s')
	SPATIAL_ICA			= 1; 
	fprintf('Using SPATIAL_ICA\n\n');
elseif  strcmp(mode,'t')
	TEMPORAL_ICA		= 1;
	fprintf('Using TEMPORAL_ICA\n\n');
end;

%---------------------------------------------------------------------	

% Initialise figures.
% create_figures;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAKE DATA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

SYNTH_DATA	= 1;
LR_DATA		= 0;
BUCHEL_DATA	= 0;

if SYNTH_DATA
	nr=40; nc=nr; 
	%[m U V]=datat_skew(nr,nc,nr*nc);  % FOR SKEW % Each col of m contains an image.
    [m sm_sources tm_sources]=datat(nr,nc,nr*nc); % FOR SIN
    m = stripmean(m,'st');
	
	do_svd=0;
	if do_svd
    fprintf('\nPerforming Singular Value Decomposition (SVD), wait...\n\n');
    [U D V] = svd(m, 0);
	else
	U=sm_sources;
	V=tm_sources;
    D = eye(neig);
	end;

	U0=U;V0=V;D0=D;
	U_orig=U;V_orig=V;D_orig=D;
	U=U_orig; V=V_orig; D=D_orig;
	pcs=U; pct=V; sval=D;
elseif LR_DATA
	old=jcd('Jim_Stone_G3:'); load jim; jcd(old);
	nr=128; nc=nr; slice=4;
	pcs = reshape(Ujim,[nr nc 10 10]);
	pcs = pcs(:,:,slice,1:10);
	pcs=squeeze(pcs);
	pcs=reshape(pcs,[nr*nc 10]);
	jsize(pcs,'pcs');pr;
	pct =Vjim;
	D   =Djim*eye(10);
	sval=D;
elseif BUCHEL_DATA
	read_buchel_data; % See prepare_buchel_data.m
	neig=10;
	pcs=U;
	pct=V;
	sval=D;
end;
% REAL DATA!!!!
%load section
%U0=U;V0=V;D0=D;
%U_orig=U;V_orig=V;D_orig=D;
%U=U_orig; V=V_orig; D=D_orig;
%pcs=U; pct=V; sval=D;

% WOULD DO PCA AT THIS POINT ... BUT SKIPPED HERE BECAUSE PCA SIMPLY RETURNS
% LINEAR MIXTURE OF 4 ORIGINAL SOURCES -  SO USE THESE AS INPUT TO ICA.

% SET SPATIAL AND TEMPORAL EIGENVECTORS.
pcs 	= pcs(:, 1:neig); 
pct 	= pct(:, 1:neig); 
sval 	= sval(1:neig, 1:neig);

% PLOT RESULTS OF PCA.
% For plotting eigvecs ...
nr_plots = floor(sqrt(neig));	% num plots in row.
nc_plots = floor(sqrt(neig));	% num plots in column.

% LIMIT NUMBER OF PLOTS.
if nr_plots>2 nr_plots=2; end;
if nc_plots>2 nc_plots=2; end;
nplots=nc_plots*nr_plots;

% PLOT SPATIAL PC == each col of pcs. 
plot_pc=1;
if plot_pc
jfig(1); clf;
jsize(pcs(:,1),'pcs(:,1)');
for k = 1:nplots
	subplot(nr_plots, nc_plots, k);
	pnshow( reshape(sm_sources(:, k), nr, nc) );
	title( ['Source Image #', int2str(k)] );
	axis off; axis square
end
% fprintf('Figure 1: Original source images.');pr;

% PLOT TEMPORAL PC == EACH COL OF pct.
jfig(2); clf;
for k = 1:nplots
	subplot(nr_plots, nc_plots, k);
	plot( tm_sources(:, k) ); axis tight
	title( ['Source Sequence #', int2str(k)] );
end
fprintf('Figures 1 and 2 are spatial and temporal eigenvectors.\n'); 
end; % plot_pc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% PREPARE WEAK MODEL DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Find covariance matrix V=x'x and U=z'z where z is short term difference from mean. 
shf=4;
lhf = 1000;
jsize(pct,'pct');
if wm_sm
	[Ux Vx]=get_UxVx(pct',shf,lhf);
	% Want to maximise (w'Vx w) / (w'Ux w) so define ...
	C1=Vx;
	C2=Ux;
else
	C1=[]; C2=[];
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% ICA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if init_W
	% Initialise weight matrix.
	% Must add rand else W0 is singular and cannot be inverted.
	W0 = eye(nic,neig)+randn(nic,neig)*0.9;
	%W0 = 1:nic*neig;
   	% W0 = W0/sum(W0); W0 = reshape(W0,[nic neig]);
else
	warning('Starting from old W0 ...');
	W0=reshape(wmin,size(W0));
	jsize(W0,'W0');
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAXMIMISE ENTROPY.

	% Ensure that x=PQ; and P=kxm, Q=kxn.
	if init==1
		% Make svals sum to one. 
		d = diag(D);
		sval = D(1:neig, 1:neig)./sum(d(1:neig));
		% figure(10); plot(d(1:neig)); drawnow;
		sval = D(1:neig, 1:neig);

		P=pcs*sqrt(sval);
		Q=pct*sqrt(sval);
	end;
	jsize(P,'P');jsize(Q,'Q');
	% stica
	alpha = 0.5;

	[V1, d, S1, T1, w] = stica(P, Q, alpha, W0, hi_kurt_s,lo_kurt_s,hi_kurt_t,lo_kurt_t,C1,C2);

jtoc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

