%% simulated datasets (a total of 22 artficial datasets)
% all using MATLAB's randn function and standard procedures

%% Data with no structure (the n-cluster data)
Ns = 100 ; % the same as 'p' in the main manuscript
Nvox = 1000 ; % the same as 'n' in the main manuscript
X = randn(Nvox, Ns) ;

%% (compact) Data with 1 cluster (the 1-cluster data)
Ns = 100 ; % the same as 'p' in the main manuscript
Nvox = 1000 ; % the same as 'n' in the main manuscript
Sigma = 1 ;
X = repmat(randn(Ns,1), [1 Nvox]) + (Sigma * randn(Ns, Nvox)) ;
X = X';

%% Data with known Copt clusters (noise = low) (total = 10 sets)
% Nc values are arbitrary ('rj' in the paper)
% under the condition sum(Nc) = Nvox. 
% Uncomment Nc value with the desired Copt 
Sigma = 1 ;

Ns = 100 ;
% Nc = [600 400] ; % Copt = 2 clusters
% Nc = [500 50 450] ; %Copt=3 clusters (only used to test the impact of m)
% Nc = [200 500 300] ; % Copt =  3 clusters
% Nc = [200 300 150 350] ; % Copt =  4 clusters
% Nc = [200 160 200 180 260] ; % Copt =  5 clusters
% Nc = [180 140 160 180 220 120] ; % Copt =  6 clusters
Nc = [140 180 150 140 120 120 150] ; % Copt =  7 clusters
% Nc = [140 120 100 120 100 170 100 150] ;% Copt =  8 clusters
% Nc = [100 110 100 120 100 140 100 130 100] ;% Copt =  9 clusters
% Nc = [100 110 90 80 120 80 80 80 160 100] ;% Copt =  10 clusters
% Nc = [100 110 90 80 100 80 80 80 90 90 100] ;% Copt =  11 clusters
X = [] ;
X_orig = randn(Ns,length(Nc)) ;
while nnz(max(moh_corr(X_orig)-eye(length(Nc))) > 0.1) > 0
    X_orig = randn(Ns,length(Nc)) ;
end
for cs=1:length(Nc)
    X = [X, repmat(X_orig(:,cs), [1 Nc(cs)])];
end
X = X + (Sigma * randn(Ns, sum(Nc))) ; 
X = X' ;

%% Repeat the same datasets with noise = high (total  = 10 sets)
Sigma = 4;
% repeat above.


%% function used to calculate correlations in Line 35 above (fast)
% % compute correlation coefficient
% % cc = moh_corr(a)
% function cc = moh_corr(a)
% 
% N = size(a,1) ;
% P = size(a,2) ;
% 
% aa = (a - repmat(sum(a)/N, N,1))'*(a - repmat(sum(a)/N, N,1)) ./(N-1) ;
% cc = aa ./ (repmat(std(a), P,1)'.*repmat(std(a), P,1)) ;


%% Core of classic FCM algorithm
% This is the same algorithm of any popular FCM analysis
% (also available in Matlab's Fuzzy Logic toolbox)
% Key steps are: V, D and U calculations

% X = fMRI data
% Nvox = size(X, 1) ; % number of voxels
% pdim = size(X, 2) ; % number of timepoints 
% c = number of expected clusters
% m = degree of fuzziness

% thi is an iterative algorithm: 
% at each itiration
Uinit = U ; % initial (step k-1 of FCM)

% centroid calculation (prototypes V)
V = (power(Uinit,m)*X) ./ (repmat(sum(power(Uinit',m))', [1,pdim])) ;

% similarity (distance D)
cc = (X' - repmat(sum(X')/pdim, pdim,1))'*(V' - repmat(sum(V')/pdim, pdim,1)) ./(pdim-1) ;
cc = cc ./ (repmat(std(X'), c,1)'.*repmat(std(V'),Nvox ,1)) ;
sqcc = sign(cc).*sqrt(abs(cc)) ;
D = (1 - sqcc) ./ (1 + sqcc) ;

% Updated membership calculation
U = ((1./power(D,2/(m-1))) ./ repmat(sum((1./power(D,2/(m-1)))'), c,1)')' ;

% update the value of FCM's objective function:
Jm = sum(sum( power(U,m) .* power(D',2))) ;

% criterion of the iterative procedure ; stop/converge when:
sum(power(U(:)-Uinit(:), 2)) < precision ; % e.g. precision = 0.0001