NIRS_ICA_V1.0/ica_algorithms/tentative/STICA_skew_demo/STICA_skew_demo/STICA_NIPS_ADDD_ONS/prepare_buchel_data.m
% jcorr_many(Desjim(:,1:4),yt)'

clear all;

% BUCHEL DATA.
% See  prepare_buchel;
old=jcd('Jims_working_disc:Jims_Matlab5_files:ICA:icademo_for_salk_buchel:');
load slices.mat; % slices     360x53x63    9616320  double array
jcd(old);[nt nr nc]=size(slices);

mixt=reshape_slices(slices,'s'); % 3339 rows  360 cols.
clear slices;

REMOVE_STATIONERY_CONDITION=1;
if REMOVE_STATIONERY_CONDITION
	ind=[1:80 91:170 181:260 271:350];
	mixt = mixt(:,ind); % jsize =  3339  320 
	[temp nt]=size(mixt);
end;

jsize(mixt); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Remove data from centre of image.
REMOVE=0;
if REMOVE
im=mixt(:,1);
n=11;
nnr=7; % height 
nnc=11; % width
r0=27;
c0=32;
im=zeros(nr,nc);
for r=r0-nnr:r0+nnr 
	for c=c0-nnc:c0+nnc 
		im(r,c)=1; 
	end; 
end; 
f(1); clf; pnshow(im); drawnow;
indblank=find(im(:));

% SET TO MEAN VALUE OF DATA.
ablate=mixt(indblank,:);
m=mean(ablate(:));
mixt(indblank,:)=m;
jsize(mixt);
im=mixt(:,1); im=reshape(im,nr,nc); f(2); clf; pnshow(im); drawnow;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SET VAR TO UNITY?

uvar = 0;

if uvar==1
tic;
if mode=='t'
	% Set each seq (row) to have zm and unitvar.
	temp=mixt'; clear mixt;
	temp=zmean_uvar_col(temp);
	mixt=temp'; clear temp;
else
	% Set each image (row) to have zm and unitvar.
	mixt=zmean_uvar_col(mixt);
end;
fprintf('... zmean_uvar_col\n');
jtoc; 
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% REMOVE MEAN FROM TIME AND/OR SPACE.
mixt = stripmean(mixt, 'st'); % s, t, st

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic;
fprintf('Doing SVD ...\n');
jsize(mixt,'mixt'); % jsize =  3339  320 
[U D V] = svd(mixt,0);
fprintf('... SVD done\n');
jtoc; % Elapsed time = 2.48  mins (G3).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f(1);clf;
a=U(:,2);
b=reshape(a,nr,nc);
pnshow(b); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SAVE SVD TO FILE.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fname = get_fname_ica(mode);
fname = 'res_UDV';

% SAVE
eval(['save ' eval('fname') ' U D V nt nr nc' ]);
fprintf('Saved to fname = %s\n',fname);
% save fname U D V nt nr nc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOT RESULTS OF PCA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pcs=U;pct=V;sval=D;

% Select number of eig vecs sufficient to account for 0.99 of data variance.

evals = diag(sval).^2; % eval is square of sval and former == variance.
% evals=evals/sum(evals);
tot_sval = sum(evals(:));
num_evals = length(evals);

f(5);clf;plot(cumsum(evals)); 
title('eigenvalues');

sub_tot = 0;
for i=1:num_evals
	sub_tot = sub_tot + evals(i);
	prop = sub_tot/tot_sval;
	if i==neig neig_prop = prop;end;
	if prop>=0.95; n95=i; end;
end;

fprintf('neig=%d eigs accounts for %.3f of the data\n',neig,neig_prop);
fprintf('Number of evecs required to account for 0.95 of data variance = %d\n',n95);
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;

% Note: pnshow displays positive as green-scale, negative as red-scale;
% gamma factor can be used to emphasise detail (help pnshow for details)

% 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(pcs(:, k), nr, nc) );
	title( ['Eigen Image #', int2str(k)] );
	axis off; axis square
end

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

jcd(old);
cd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
