%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Implementation of C-ICT
% Developed by  M.A.B.S. Akhonda, Yuri Levin-Schwartz, and Chunying Jia
% Coded by Chunying Jia (chunyin1 at umbc.edu)
% 
% Reference:
% Jia, Chunying, M.A.B.S. Akhonda, Yuri Levin-Schwartz, Qunfang Long, Vince D. Calhoun, and Tülay Adali. 2021. 
% "Consecutive Independence and Correlation Transform for Multimodal Data Fusion: Discovery of One-to-Many Associations 
% in Structural and Functional Imaging Data" Applied Sciences 11, no. 18:
% 8382. https://doi.org/10.3390/app11188382; &
%
% M. A. B. S. Akhonda, Y. Levin-Schwartz, S. Bhinge, V. D. Calhoun and T. Adali,
% "Consecutive Independence and Correlation Transform for Multimodal Fusion: Application to Eeg and Fmri Data,"
% 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Calgary, AB, Canada, 2018,
% pp. 2311-2315, doi: 10.1109/ICASSP.2018.8462031.

%% Step 1: ICA 
% Inputs:
%
% 1. data1 - 2D matrix (subjects by volume) for the 1st dataset
% 2. data1 - 2D matrix (subjects by volume) for the 1st dataset
% 3. data1 - 2D matrix (subjects by volume) for the 1st dataset
% 4. numICs_data1 - number of Sources (components by volume) for the 1st dataset
% 5. numICs_data1 - number of Sources (components by volume) for the 2nd dataset
% 6. numICs_data1 - number of Sources (components by volume) for the 3rd dataset

% Outputs:
% 1. A_data1 - estimated mixing matrix (components by components) for the 1st dataset
% 2. A_data2 - estimated mixing matrix (components by components) for the 2nd dataset
% 3. A_data3 - estimated mixing matrix (components by components) for the 3rd dataset
% 4. icasig_data1 - sources (components by volume) for the 1st dataset
% 5. icasig_data2 - sources (components by volume) for the 2nd dataset
% 6. icasig_data3 - sources (components by volume) for the 3rd dataset
[A_data1,A_data2,A_data3,icasig_data1,icasig_data2,icasig_data3] = ...
    ica_3datasets(data1,data2,data3,numICs_data1,numICs_data2,numICs_data3);

%% Step2: Artifact Elimination
% Artifact components are identified and eliminated from sources (icasig_data1/2/3) to get:
% 1. idx_ICs_data1 - indices of components retained for the 1st dataset
% 2. idx_ICs_data2 - indices of components retained for the 2nd dataset
% 3. idx_ICs_data3 - indices of components retained for the 3rd dataset

%% Step 3: IVA

A_data1_reduced = A_data1(:,idx_ICs_data1);
A_data2_reduced = A_data2(:,idx_ICs_data2);
A_data3_reduced = A_data3(:,idx_ICs_data3);

numIVs = min([size(A_data1_reduced,2),size(A_data2_reduced,2),size(A_data3_reduced,2)]);    % number of independent vectors 

pca_A_data1 = A_data1_reduced';
[V_A_data2,~] = icatb_calculate_pca(A_data2_reduced,numIVs, 'whiten', false);
pca_A_data2 = V_A_data2'*A_data2_reduced';
[V_A_data3,~] = icatb_calculate_pca(A_data3_reduced,numIVs, 'whiten', false);
pca_A_data3 = V_A_data3'*A_data3_reduced';

pca_A_data1 = zscore(pca_A_data1')';
pca_A_data2 = zscore(pca_A_data2')';
pca_A_data3 = zscore(pca_A_data3')';

X = zeros([size(pca_A_data1) 3]);
X(:,:,1) = pca_A_data1;
X(:,:,2) = pca_A_data2;
X(:,:,3) = pca_A_data3;

numRuns_IVA = 50;   % number of runs
W_IVA = zeros(numIVs,numIVs,3,numRuns_IVA);
for nRun = 1:numRuns_IVA
    W_IVA(:,:,:,nRun) = iva_second_order(X);
end
bestRun = crossISIidxIVA(W_IVA);
W = W_IVA(:,:,:,bestRun);

A_IVA_data1 = pinv(squeeze(W(:,:,1)));
A_IVA_data2 = pinv(squeeze(W(:,:,2))*V_A_data2');
A_IVA_data3 = pinv(squeeze(W(:,:,3))*V_A_data3');

%% Step 4: Tracing Back to Components
idxCICT = zeros(3,numIVs);
for nIV = 1:numIVs
    [~, idxCICT(1,nIV)] = max(abs(A_IVA_data1(:,nIV)));
    [~, idxCICT(2,nIV)] = max(abs(A_IVA_data2(:,nIV)));
    [~, idxCICT(3,nIV)] = max(abs(A_IVA_data3(:,nIV)));
end