%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Implementation of complete model order selection using IVA (CMI-IVA)
% Developed by  M.A.B.S. Akhonda
% Coded by M.A.B.S. Akhonda (mo32 at umbc.edu)
% 
% Reference:
%  M.A.B.S. Akhonda, B.  Gabrielson, S.  Bhinge, V. Calhoun, and T. Adali,
%  “Disjoint Sub-spaces for Common and Distinct Component Analysis: Application to the Fusion of Multi-task fMRI Data.” 
%  Journal of Neuroscience Methods, vol. 358, pages 109214, (2021).


function [w_c,w_d,S_c,S_d,data_dim,w_iva,w_jICA,w_ica] = DS_ICA(Data,C,N)
%%Input
% Data => Dimensionality reduced or not. If not reduce the dimension using order N
% C => Common Subspace order
% N => Total number of components 

%% Output
% W_c,W_d => W for common and distinct subspaces
% S_c,S_d => Common and Distinct components
% Data_d => IVA transformed data or 1st level sources
% W_iva, W_j, W_ica => Saved for debugging purpose

[K] = size(Data,3); % Number of Datasets 
num_comp1 = N; 

% %% Dim reduction, un-comment it if needed
% Num_Comp = N;%max([d1,d2,d3]);
% for i = 1:K
% [U2(:,:,i), D2(:,:,i), V2(:,:,i)]=svd(Data(:,:,i),'econ');
% X_d(:,:,i)=(inv(D2(1:Num_Comp,1:Num_Comp,i)))*U2(:,1:Num_Comp,i)'*Data(:,:,i);
% end
% If Data is dimensionality reduced use X_d = Data
x_d = Data;

%% CCA for 2 datasets and MCCA for more than 2 datasets
%  [C1,C2,r,E1,E2] = canoncorr(Data(:,:,1)',Data(:,:,2)');
%  W_mcca(:,:,1) = C1;
%  W_mcca(:,:,2) = C2;
%W_mcca = mcca(X_d,N,'maxvar');

%% IVA to idenity common space for more than 2 datasets
w_iva = iva_second_order(x_d);
for i = 1:K
    data_dim(:,:,i) = w_iva(:,:,i)*x_d(:,:,i);
end

%% jICA for common subspace
temp = data_dim(1:C,:,1);
for i = 2:K
temp = cat(2,temp,data_dim(1:C,:,i));
end
x_jIca = temp;

[w_jICA] = ICA_EBM(x_jIca); % Change ICA_EBM to use other algorithms

for i = 1:K
    w_c(:,:,i) = pinv(w_iva(:,1:C,i)')*pinv(w_jICA);
    %A_c(:,:,i) = (pinv(U2(:,1:Num_Comp,i)')*(D2(1:Num_Comp,1:Num_Comp,i)))*W_c(:,:,i);
end

for i = 1:K
    S_c(:,:,i)  = w_jICA*data_dim(1:C,:,i);
end

%% Distinct ICAs for distinct subspace  
for i = 1:K
x_dis(:,:,i) = pinv(w_iva(:,C+1:num_comp1,i)')* data_dim(C+1:num_comp1,:,i);
end
% Dim reduction to remove the redundent space
Num_Comp2 = num_comp1-C;
for i = 1:K
[U1(:,:,i), D1(:,:,i), V1(:,:,i)]=svd(x_dis(:,:,i),'econ');
data_dis(:,:,i)=(inv(D1(1:Num_Comp2,1:Num_Comp2,i)))*U1(:,1:Num_Comp2,i)'*x_dis(:,:,i);
end
% ICA to estimate the distinct components
for i = 1:K
    [w_ica(:,:,i)] = ICA_EBM(data_dis(:,:,i));
end
for i = 1:K
    w_d(:,:,i) = (pinv(U1(:,1:Num_Comp2,i)')*(D1(1:Num_Comp2,1:Num_Comp2,i)))*pinv(w_ica(:,:,i));
   % A_d(:,:,i) = (pinv(U2(:,1:N,i)')*(D2(1:N,1:N,i)))*W_d(:,:,i);
end
for i = 1:K
    S_d(:,:,i) = w_ica(:,:,i)*data_dis(:,:,i);
end
end