%2013, 6,14, duyuhui
%Originally from srinivas. And I made some comments and revision.

function  [A, icasig] = useICASSO(data, algorithm, useICAOptions, icasso_opts)
% Inputs:
% 1. data - Data of dimensions components x voxels
% 2. algorithm - Algorithm name. Please see icatb_icaAlgorithm for more
% info.
% 
% Outputs:
% A - Mixing matrix of dimensions components x components
% icasig - ICA maps of dimensions components x voxels

data = double(data);

numOfIC = size(data, 1);

if (~exist('algorithm', 'var'))
    algorithm = 1;
end

if (~exist('icasso_opts', 'var'))
    
    %icasso_opts = struct('sel_mode', 'randinit', 'num_ica_runs', 10); 
    icasso_opts = struct('sel_mode', 'both', 'num_ica_runs', 20);%duyuhui, 20 should be times of the group level ICA
end

if (~exist('useICAOptions', 'var'))
    useICAOptions = 0;
end

ICA_Options = {};

if (useICAOptions)
    ICA_Options = icaOptions(size(data), algorithm, 'on');
end


%%%%% Calculate PCA and Whitening matrix %%%%%
% PCA
[V, Lambda] = v_pca(data, 1, numOfIC, 0, 'transpose', 'yes');

% Whiten matrix
[w, White, deWhite] = v_whiten(data, V, Lambda, 'transpose');

clear V Lambda;

sR = GIGICA_icassoEst(icasso_opts.sel_mode, data, icasso_opts.num_ica_runs, 'numOfPC', numOfIC, 'algoIndex', algorithm, ...
    'dewhiteM', deWhite, 'whiteM', White, 'whitesig', w, 'icaOptions', ICA_Options);

clear data w deWhite White;

%%%% Visualization %%%%%%

sR = GIGICA_icassoExp(sR);

%%% Visualization & returning results
%%% Allow to disable visualization
disp(['Launch Icasso visualization supposing ', num2str(numOfIC), ' estimate-clusters.']);
disp('Show demixing matrix rows.');
% icassoShow(sR, 'L', numOfIC, 'estimate', 'demixing');

disp(['Launch Icasso visualization supposing ', num2str(numOfIC), ' estimate-clusters.']);
disp('Show IC source estimates (default), reduce number of lines');
disp('Collect results.');
% iq = icassoShow(sR, 'L', numOfIC, 'colorlimit', [.8 .9]);


minClusterSize = ceil(0.8*icasso_opts.num_ica_runs);

maxClusterSize = icasso_opts.num_ica_runs;

[metric_Q, A, W, icasig] = getStableRunEstimates(sR, minClusterSize, maxClusterSize);



function [metric_Q, A, W, icasig] = getStableRunEstimates(sR, minClusterSize, maxClusterSize)
%% Get stable run based on code by Sai Ma. Stable run estimates will be used instead of centrotype
%

% number of runs and ICs
numOfRun = length(sR.W);
numOfIC = size(sR.W{1},1);

% Get the centrotype for each cluster and Iq
index2centrotypes = icassoIdx2Centrotype(sR,'partition', sR.cluster.partition(numOfIC,:));
Iq = icassoStability(sR, numOfIC, 'none');

% Find IC index  within each cluster
partition = sR.cluster.partition(numOfIC, :);
clusterindex = cell(1, numOfIC);
for i = 1:numOfIC
    temp = (partition == i);
    clusterindex{i} = sR.index(temp, :);
    clear temp;
end

% Compute stability metric for each run within each cluster
eachRun = zeros(numOfRun, numOfIC);
qc = 0; % num of qualified clusters
for i = 1:numOfIC
    thisCluster = (clusterindex{i}(:,1))';
    clusterSize = length(clusterindex{i});
    if ((clusterSize >= minClusterSize) && (clusterSize <= maxClusterSize) && (Iq(i)>=0.7))
        qc = qc + 1;
        for k = 1:numOfRun
            thisRun = find(thisCluster == k);
            ICindex = (clusterindex{i}(thisRun,1)-1)*numOfIC + clusterindex{i}(thisRun,2);
            if ~isempty(thisRun)
                eachRun(k,i) = max(sR.cluster.similarity(index2centrotypes(i),ICindex'));
            end
            clear thisRun ICindex;
        end
    end
    clear thisCluster clusterSize;
end

%% Find stable run
metric_Q = sum(eachRun,2)/qc;
[dd, stableRun] = max(metric_Q);

%% Get stable run estimates
W = sR.W{stableRun};
clusters_stablerun = partition((stableRun - 1)*numOfIC + 1 : stableRun*numOfIC);
[dd, inds] = sort(clusters_stablerun);
W = W(inds, :);
A = pinv(W);
icasig = W*sR.signal;

