function [A,icasig] = separateICA(data,numICs,numRuns)
%% ICA is performed on each dataset
% Inputs:
%
% 1. data - 2D matrix components by volume
% 3. numIC - number of ICs
% 4. numRuns - number of runs

% Outputs:
%
% 1. A - Estimated mixing matrix (components by components)
% 2. icasig - Sources (components by volume)

% PCA
V = icatb_calculate_pca(data',numICs, 'whiten', false);
pca_data = V'*data;

% ICA
tempW = zeros(numICs,numICs,numRuns);
for nRun = 1:numRuns
    [weights,sphere] = icatb_runica(pca_data);
    tempW(:,:,nRun) = weights*sphere;
end
bestRun = crossISIidx(tempW); 
W = tempW(:,:,bestRun);
icasig_tmp = W*pca_data;
A_tmp = pinv(W*V');
[A, ~, icasig] = changeSignOfComponents(A_tmp,icasig_tmp);