diff --git a/icatb_estimate_dimension.m b/icatb_estimate_dimension.m new file mode 100644 index 0000000..fb982e3 --- /dev/null +++ b/icatb_estimate_dimension.m @@ -0,0 +1,1005 @@ +function [comp_est_AIC, comp_est_KIC, comp_est_MDL, mdl, aic, kic] = icatb_estimate_dimension(files, maskvec, precisionType, dim_est_opts) +% function [comp_est, mdl, aic, kic] = icatb_estimate_dimension(files, maskvec) +% Select the order of the multivariate data using Information theoretic +% criteria with the option of sample dependence correction; +% Inputs: +% +% 1. files - fullfile path of data files as a character array or +% numeric array of size xdim by ydim by zdim by tdim. +% 2. maskvec - Mask file in analyze or Nifti format. Mask must have the +% same dimensions as the data or mask vector must be in logical vector format of +% of length xdim*ydim*zdim. +% 3. precisionType - Load data as double or single +% 4. dim_est_opts - Dimensionality options. Options to select method and if +% you specified method = 2 set fwhm variable as well. +% +% Output: +% 1. comp_est_AIC: estimated order using AIC +% 2. comp_est_KIC: estimated order using KIC +% 3. comp_est_MDL: estimated order using MDL +% 4. mdl - MDL vector +% 5. aic - AIC vector +% 6. kic - KIC vector + +% Please cite the following paper if you use this code for publication +% Y.-O. Li, T. Adali and V. D. Calhoun, "Estimating the number of independent +% components for fMRI data," Human Brain Mapping, vol. 28, pp. 1251--66, 2007 + + +%COPYRIGHT NOTICE +%This function is a part of GIFT software library +%Copyright (C) 2003-2009 +% +%This program is free software; you can redistribute it and/or +%modify it under the terms of the GNU General Public License +%as published by the Free Software Foundation; either version 2 +%of the License, or any later version. +% +%This program is distributed in the hope that it will be useful, +%but WITHOUT ANY WARRANTY; without even the implied warranty of +%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +%GNU General Public License for more details. +% +%You should have received a copy of the GNU General Public License +%along with this program; if not, write to the Free Software +%Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +icatb_defaults; +global DIM_ESTIMATION_OPTS; + +if (~exist('files', 'var')) + files = icatb_selectEntry('typeSelection', 'multiple', 'filter', '*.img', 'title', 'Select functional data files'); +end + +drawnow; + +if (~exist('maskvec', 'var')) + maskvec = []; +end + +if (~exist('precisionType', 'var')) + precisionType = 'double'; +end + +dim_est_method = 1; +fwhm = []; + +if (~exist('dim_est_opts', 'var') || isempty(dim_est_opts)) + dim_est_opts = DIM_ESTIMATION_OPTS; +end + +try + dim_est_method = dim_est_opts.method; +catch +end + + +% try +% dim_est_method = DIM_ESTIMATION_OPTS.method; +% catch +% if (isfield(DIM_ESTIMATION_OPTS, 'iid_sampling')) +% if (DIM_ESTIMATION_OPTS.iid_sampling == 1) +% % iid sampling on +% dim_est_method = 1; +% else +% % iid sampling off +% dim_est_method = 2; +% end +% end +% end + +try + fwhm = dim_est_opts.fwhm; +catch +end + +if (dim_est_method == 2) + if (isempty(fwhm)) + error('Please set fwhm smoothness factor using DIM_ESTIMATION_OPTS.fwhm when DIM_ESTIMATION_OPTS.method = 2'); + end +end + +if (ischar(maskvec) && strcmp(maskvec(end-3:end), '.mat')) + mask = load(maskvec); + maskvec = mask.mask; +else + % If mask is supplied as a file, read data and convert data to + % logical + maskvec = icatb_rename_4d_file(maskvec); + % Load only the first file + maskvec = deblank(maskvec(1,:)); + maskvec = icatb_loadData(maskvec); + + maskvec = icatb_loadData(deblank(maskvec(1, :))); +end + +%% Get dimensions of data by loading first time point +if (ischar(files) && ~strcmp(files(end-3:end), '.mat')) + files = icatb_rename_4d_file(files); + first_file = deblank(files(1, :)); + data = icatb_loadData(first_file); + xdim = size(data, 1); ydim = size(data, 2); zdim = size(data, 3); tdim = size(files, 1); + first_file = icatb_parseExtn(first_file); + if (strcmpi(first_file(end-3:end), '.mat')) + tdim = size(data, 4); + end +else + if (ischar(files) && strcmp(files(end-3:end), '.mat')) + data = load(files); + data = data.data; + else + data = files; + end + + if (length(size(data)) == 4) + xdim = size(data, 1); ydim = size(data, 2); zdim = size(data, 3); tdim = size(data, 4); + else + if (length(size(maskvec)) == 3) + xdim = size(maskvec, 1); ydim = size(maskvec, 2); zdim = size(maskvec, 3); tdim = size(data, 2); + if (size(data, 1) ~= xdim*ydim*zdim) + error('Mask dimensions (x, y, z) don''t match data dimensions (x, y, z)'); + end + else + error('Cannot get the x,y,z dimension information from the mask or the data. Data must be 4D array or mask must be 3D array'); + end + end + +end + +%% Load mask and check the dimensions of the data +if (isempty(maskvec)) + + if (length(size(data)) == 4) + tmp = data(:, :, :, 1); + elseif ((length(size(data)) == 3) && (numel(data) == xdim*ydim*zdim)) + tmp = data(:); + elseif ((length(size(data)) == 2) && (numel(data) == xdim*ydim*zdim)) + tmp = data(:); + else + tmp = data(:, 1); + end + + + maskvec = (tmp(:) >= abs(mean(tmp(:)))); + +else + + % Get mask vector + maskvec = (maskvec ~= 0); + maskvec = maskvec(:); + +end + +if (length(maskvec) ~= xdim*ydim*zdim) + error('Error:MaskDIM', 'Mask dimensions (%d) doesn''t match the data (%d)', length(maskvec), xdim*ydim*zdim); +end + +%% Load data by reading timepoint at a time +if (ischar(files) && ~strcmp(files(end-3:end), '.mat')) + data = icatb_read_data(files, [], find(maskvec == 1), precisionType); +else + if (strcmpi(precisionType, 'single')) + data = single(data); + end + data = reshape(data, xdim*ydim*zdim, tdim); + data = data(maskvec, :); +end + +verbose = 1; + +mdl = []; +aic = []; +kic = []; + +if ((dim_est_method == 3) || (dim_est_method == 4)) + % Use ER-FM or ER-AR + + if (verbose) + disp('Using entropy rate estimator ...'); + end + + [ER_FM, ER_AR] = orderEstimation_HdsvsEig_R(data'); + + if (dim_est_method == 3) + disp('Order estimated by entropy rate based method assumes signal has finite memory length ...'); + comp_est = ER_FM; + else + disp('Order estimated by entropy rate based method assumes AR signal ...'); + comp_est = ER_AR; + end + + return; + +end + + +if (dim_est_method == 1) + + %% Perform variance normalization + if verbose + fprintf('\n Performing variance normalization ...'); + end + + for n = 1:size(data, 2) + data(:, n) = detrend(data(:, n), 0) ./ std(data(:, n)); + end + + %% (completed) + + fprintf('\n P1:'); + [V, EigenValues] = icatb_svd(data, tdim, 'verbose', 0); + EigenValues = diag(EigenValues); + EigenValues = EigenValues(end:-1:1); + dataN = data*V(:, end:-1:1); + clear V; + + %% Select Gaussian principal components for effectively i.i.d. sample estimation + if verbose + fprintf('\n Selecting Gaussian principal components ...'); + end + + %% Using 12 gaussian components from middle, top and bottom gaussian components to determine the subsampling depth. Final subsampling depth is determined using median + kurtv1 = kurtn(dataN); + kurtv1(EigenValues > mean(EigenValues)) = 1000; + idx_gauss = find(((kurtv1(:) < 0.3) .* (EigenValues > eps)) == 1); + idx = idx_gauss(:)'; + dfs = length(find(EigenValues > eps)); % degrees of freedom + minTp = 12; + + if (length(idx) >= minTp) + middle = round(length(idx)/2) + 1; + idx = [idx(1:4), idx(middle - 1:middle + 2), idx(end-3:end)]; + elseif isempty(idx) + minTp = min([minTp, dfs]); + idx = (dfs - minTp + 1:dfs); + end + + idx = unique(idx); + + %% Estimate the subsampling depth for effectively i.i.d. samples + if verbose + fprintf('\n\n Estimate effectively i.i.d. samples: '); + end + + mask_ND = reshape(maskvec, xdim, ydim, zdim); + ms = length(idx); + s = zeros(ms,1); + for i = 1:ms + x_single = zeros(xdim*ydim*zdim, 1); + x_single(maskvec) = dataN(:, idx(i)); + x_single = reshape(x_single, xdim, ydim, zdim); + s(i) = est_indp_sp(x_single); + if (i > 6) + tmpS = s(1:i); + tmpSMedian = round(median(tmpS)); + if (length(find(tmpS == tmpSMedian)) > 6) + s = tmpS; + break; + end + end + dim_n = prod(size(size(x_single))); + clear x_single; + end + clear dataN; + fprintf('\n'); + s1 = round(median(s)); + if floor((sum(maskvec)/(1*tdim))^(1/dim_n)) < s1 + s1 = floor((sum(maskvec)/(1*tdim))^(1/dim_n)); + end + N = round(sum(maskvec)/s1^dim_n); + %% (completed) + + %% Use the subsampled dataset to calculate eigen values + if verbose + fprintf('\n Perform EVD on the effectively i.i.d. samples ...'); + end + if s1 ~= 1 + mask_s = subsampling(mask_ND, s1, [1,1,1]); + mask_s_1d = reshape(mask_s, 1, prod(size(mask_s))); + dat = zeros(length(find(mask_s_1d == 1)), tdim); + for i = 1:tdim + x_single = zeros(xdim*ydim*zdim, 1); + x_single(maskvec) = data(:, i); + x_single = reshape(x_single, xdim, ydim, zdim); + dat0 = subsampling(x_single, s1, [1, 1, 1]); + dat0 = reshape(dat0, prod(size(dat0)), 1); + dat(:, i) = dat0(mask_s_1d); + clear x_single; + end + clear data; + %% Perform variance normalization + for n = 1:size(dat, 2) + dat(:, n) = detrend(dat(:, n), 0) ./ std(dat(:, n)); + end + %% (completed) + fprintf('\n P2:'); + [V, EigenValues] = icatb_svd(dat, tdim, 'verbose', 0); + EigenValues = diag(EigenValues); + EigenValues = EigenValues(end:-1:1); + %lam = EigenValues; + end + + lam = EigenValues; + clear dat; + + if verbose + fprintf('\n Effective number of i.i.d. samples: %d \n',N); + end + + %% Make eigen spectrum adjustment + if verbose + fprintf('\n Perform eigen spectrum adjustment ...'); + end + lam = eigensp_adj(lam, N, length(lam)); + %% (completed) + + if sum(imag(lam)) + error('Invalid eigen value found for the subsampled data.'); + end + + +else + %% No iid sampling + if verbose + fprintf('\n computing eigen values ...\n'); + end + + [V, EigenValues] = icatb_svd(data, tdim, 'verbose', 0); + EigenValues = diag(EigenValues); + EigenValues = EigenValues(end:-1:1); + lam = EigenValues; + N = ceil(size(data, 1) / prod(fwhm)); + +end + +%% Correction on the ill-conditioned results (when tdim is large, some least significant eigenvalues become small negative numbers) +lam(real(lam) <= eps) = min(lam(real(lam)>= eps)); + +if verbose + fprintf('\n Estimating the dimension ...'); +end +p = tdim; +aic = zeros(1, p - 1); +kic = zeros(1, p - 1); +mdl = zeros(1, p - 1); +for k = 1:p-1 + LH = log(prod( lam(k+1:end).^(1/(p-k)) )/mean(lam(k+1:end))); + mlh = 0.5*N*(p-k)*LH; + df = 1 + 0.5*k*(2*p-k+1); + aic(k) = -2*mlh + 2*df; + kic(k) = -2*mlh + 3*df; + mdl(k) = -mlh + 0.5*df*log(N); +end + +% Find the first local minimum of each ITC +itc = zeros(3, length(mdl)); +itc(1,:) = aic; +itc(2,:) = kic; +itc(3,:) = mdl; + +%% Estimated components using MDL +dlap = squeeze(itc(3, 2:end) - itc(3, 1:end-1)); +a3 = find(dlap > 0); + +dlap2 = squeeze(itc(2, 2:end) - itc(2, 1:end-1)); +a2 = find(dlap2>0); + +dlap1 = squeeze(itc(1, 2:end) - itc(1, 1:end-1)); +a1 = find(dlap1>0); + + +if isempty(a3) + comp_est_MDL = length(squeeze(itc(3, :))); +else + comp_est_MDL = a3(1); +end + +if isempty(a1) + comp_est_AIC = length(squeeze(itc(1, :))); +else + comp_est_AIC = a1(1); +end + +if isempty(a2) + comp_est_KIC = length(squeeze(itc(2, :))); +else + comp_est_KIC = a2(1); +end + +fprintf('\n'); +disp(['Estimated components with mdl is found out to be ', num2str(comp_est_MDL)]); + +fprintf('\n'); +disp(['Estimated components with aic is found out to be ', num2str(comp_est_AIC)]); + +fprintf('\n'); +disp(['Estimated components with kic is found out to be ', num2str(comp_est_KIC)]); +%fprintf('\n'); + +function out = subsampling(x,s,x0) +% Subsampling the data evenly with space 's' + +n = size(x); + +if max(size(n)) == 2 & min(n) == 1 % 1D + out = x([x0(1):s:max(n)]); + +elseif max(size(n)) == 2 & min(n) ~= 1 % 2D + out = x([x0(1):s:n(1)],[x0(2):s:n(2)]); + +elseif max(size(n)) == 3 & min(n) ~= 1 % 3D + out = x([x0(1):s:n(1)],[x0(2):s:n(2)],[x0(3):s:n(3)]); + +else + error('Unrecognized matrix dimension!(subsampling)'); + +end + + +function [s, entrate_m] = est_indp_sp(x) +% estimate the effective number of independent samples based on the maximum +% entropy rate principle of stationary random process + + +dimv = size(x); +s0 = 0; +fprintf('\n'); +for j = 1:min(dimv)-1 + x_sb = subsampling(x,j,[1,1,1]); + if j == 1 + fprintf('\n Estimating the entropy rate of the Gaussian component with subsampling depth %d,',j); + else + fprintf(' %d,',j); + end + entrate_m = entrate_sp(x_sb,1); + + ent_ref = 1.41; + if entrate_m > ent_ref + s0 = j; break; + end +end +fprintf(' Done;'); +if s0 == 0 + error('Ill conditioned data, can not estimate independent samples.(est_indp_sp)'); +else + s = s0; +end + + +function out = entrate_sp(x, sm_window) +% Calculate the entropy rate of a stationary Gaussian random process using +% spectrum estimation with smoothing window + +n = size(x); + +% Normalize x_sb to be unit variance +x_std = std(reshape(x,prod(n),1)); +if x_std < 1e-10; x_std = 1e-10; end; +x = x/x_std; + +if( sm_window == 1) + + M = ceil(n/10); + + if(max(size(n) >= 3)) + parzen_w_3 = zeros(2*n(3)-1,1); + parzen_w_3(n(3)-M(3):n(3)+M(3)) = parzen_win(2*M(3)+1); + end + + if(max(size(n) >= 2)) + parzen_w_2 = zeros(2*n(2)-1,1); + parzen_w_2(n(2)-M(2):n(2)+M(2)) = parzen_win(2*M(2)+1); + end + + if(max(size(n) >= 1)) + parzen_w_1 = zeros(2*n(1)-1,1); + parzen_w_1(n(1)-M(1):n(1)+M(1)) = parzen_win(2*M(1)+1); + end + +end + +if max(size(n)) == 2 & min(n) == 1 % 1D + xc = xcorr(x,'unbiased'); + xc = xc.*parzen_w; + xf = fftshift(fft(xc)); + +elseif max(size(n)) == 2 & min(n) ~= 1 % 2D + xc = xcorr2(x); % default option: computes raw correlations with NO + % normalization -- Matlab help on xcorr + + % Bias correction + v1 = [1:n(1),n(1)-1:-1:1]; + v2 = [1:n(2),n(2)-1:-1:1]; + + vd = v1'*v2; + xc = xc./vd; + parzen_window_2D = parzen_w_1*parzen_w_2'; + xc = xc.*parzen_window_2D; + xf = fftshift(fft2(xc)); + +elseif max(size(n)) == 3 & min(n) ~= 1 % 3D + xc = zeros(2*n-1); + for m3 = 0:n(3)-1 + temp = zeros(2*n(1:2)-1); + for k = 1:n(3)-m3 + temp = temp + xcorr2(x(:,:,k+m3),x(:,:,k)); % default option: + % computes raw correlations with NO normalization + % -- Matlab help on xcorr + end + xc(:,:,n(3)-m3) = temp; + xc(:,:,n(3)+m3) = temp; + end + + % Bias correction + v1 = [1:n(1),n(1)-1:-1:1]; + v2 = [1:n(2),n(2)-1:-1:1]; + v3 = [n(3):-1:1]; + + vd = v1'*v2; + vcu = zeros(2*n-1); + for m3 = 0:n(3)-1 + vcu(:,:,n(3)-m3) = vd*v3(m3+1); + vcu(:,:,n(3)+m3) = vd*v3(m3+1); + end + + xc = xc./vcu; + + parzen_window_2D = parzen_w_1*parzen_w_2'; + for m3 = 0:n(3)-1 + parzen_window_3D(:,:,n(3)-m3) = parzen_window_2D*parzen_w_3(n(3)-m3); + parzen_window_3D(:,:,n(3)+m3) = parzen_window_2D*parzen_w_3(n(3)+m3); + end + xc = xc.*parzen_window_3D; + + xf = fftshift(fftn(xc)); + +else + error('Unrecognized matrix dimension.'); + +end + +xf = abs(xf); +xf(xf<1e-4) = 1e-4; +out = 0.5*log(2*pi*exp(1)) + sumN(log(abs((xf))))/2/sumN(abs(xf)); + + +function w = parzen_win(n) +% PARZENWIN Parzen window. +% PARZENWIN(N) returns the N-point Parzen (de la Valle-Poussin) window in +% a column vector. + +% Check for valid window length (i.e., n < 0) +[n,w,trivialwin] = checkOrder(n); +if trivialwin, return, end; + +% Index vectors +k = -(n-1)/2:(n-1)/2; +k1 = k(k<-(n-1)/4); +k2 = k(abs(k)<=(n-1)/4); + +% Equation 37 of [1]: window defined in three sections +w1 = 2 * (1-abs(k1)/(n/2)).^3; +w2 = 1 - 6*(abs(k2)/(n/2)).^2 + 6*(abs(k2)/(n/2)).^3; +w = [w1 w2 w1(end:-1:1)]'; + + +function [n_out, w, trivalwin] = checkOrder(n_in) +% CHECKORDER Checks the order passed to the window functions. + +w = []; +trivalwin = 0; + +% Special case of negative orders: +if n_in < 0, + error('Order cannot be less than zero.'); +end + +% Check if order is already an integer or empty +% If not, round to nearest integer. +if isempty(n_in) | n_in == floor(n_in), + n_out = n_in; +else + n_out = round(n_in); + warning('Rounding order to nearest integer.'); +end + +% Special cases: +if isempty(n_out) | n_out == 0, + w = zeros(0,1); % Empty matrix: 0-by-1 + trivalwin = 1; +elseif n_out == 1, + w = 1; + trivalwin = 1; +end + + +function lam_adj = eigensp_adj(lam,n,p) +% Eigen spectrum adjustment for EVD on finite samples + +r = p/n; + +bp = (1+sqrt(r))^2; +bm = (1-sqrt(r))^2; + +vv = [bm:(bp-bm)/(5*p-1):bp]; + +gv = 1./(2*pi*r*vv).*sqrt((vv-bm).*(bp-vv)); + +gvd = zeros(size(gv)); +for i = 1:length(gv) + gvd(i) = sum(gv(1:i)); +end +gvd = gvd/max(gvd); + +lam_emp = zeros(size(lam)); +for i = 1:p + i_norm = i/p; + [minv,minx] = min(abs(i_norm-gvd)); + lam_emp(i) = vv(minx); +end + +lam_emp = rot90(lam_emp, 2); + +lam_adj = lam./lam_emp; + + +function [sum_dat] = sumN(dat) +% sum of the all elements of the dat matrix + +sum_dat = sum(dat(:)); + + +function c = xcorr2(a,b) +% two dimensional cross correlation + +if nargin == 1 + b = a; +end + +c = conv2(a, rot90(conj(b),2)); + +function kurt = kurtn(x) +% Normalized kurtosis funtion so that for a Gaussian r.v. the kurtn(g) = 0 +% Input: x 1:N vector + +kurt = zeros(1, size(x, 2)); + +for i = 1:size(x, 2) + a = detrend(x(:, i), 0); + a = a/std(a); + kurt(i) = mean(a.^4) - 3; +end + + + +function [ER_FM, ER_AR] = orderEstimation_HdsvsEig_R( data1 ) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%ENTROPY-RATE BASED ORDER SELECTION (ER_FM & ER_AR) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Input: +% data1: mixtures X +% Outputs: +% ER_FM: Order estimated by entropy rate based method that assumes signal has finite memory length +% ER_AR: Order estimated by entropy rate based method that assumes AR signal + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% References: +% +% ER_FM & ER_AR +% G.-S. Fu, M. Anderson, and T. Adali, Likelihood estimators for dependent samples +% and their application to order detection, in Signal Process, IEEE Trans. +% on, vol. 62, no. 16, pp. 4237-4244, Aug. 2014. +% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Program by Geng-Shen Fu +% Please contact me at fugengs1@umbc.edu +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +data1 = data1 - mean(data1,2)*ones(1,size(data1,2)); % remove the mean +[N T]=size(data1); + +% estimate the downsampling depth and order detection +[U1,S1,V1] = svd(data1,'econ'); +clear V1; +data1 = U1'*data1; + +[ACS DSD dsdAll] = Entropy_rate_estimator(data1); + +ER_FM_DS=zeros(1,max(DSD,100)); +for i=1:DSD + for n=1:N + H(n)=Entropy_rate_single_estimator(ACS(n,1:i:DSD)); + end; + if sum(H==Inf)>0 + continue; + end + [H index_H] = sort(H, 'descend'); + ER_FM_DS(i)=order_estimate_entropy(N, T/i, H); +end +ER_FM_DS1=ER_FM_DS(ER_FM_DS~=0); +ER_FM=round(median(ER_FM_DS1)); + +ER_AR_DS=zeros(1,max(DSD,100)); +for i=1:DSD + for n=1:N + H(n)=Entropy_rate_single_estimator_AR(data1(n,1:i:end)); + end; + if sum(H==Inf)>0 + continue; + end + [H index_H] = sort(H, 'descend'); + ER_AR_DS(i)=order_estimate_entropy(N, T/i, H); +end +ER_AR_DS1=ER_AR_DS(ER_AR_DS~=0); +ER_AR=round(median(ER_AR_DS1)); + +maxp0 = floor(2*size(data1,2)/size(data1,1))/2; +dsd = size(data1,1)/2; % the initial downsampling depth is half the number of sensors +existing_orders = []; % save all estimated orders +while 1 + % estimate orders + E_JDS = order_estimate_complex_c(size(data1,1), size(data1,2), dsd, diag(S1)); + existing_orders = [existing_orders; E_JDS]; + if sum(existing_orders==E_JDS)>=5 % return when the order detection converges + dsd = ceil(dsd); + break; + end + + % re-estimate downsampling depth + dsd = 0; + cnt1 = 0; + cnt2 = 0; % counting the MA orders reaching the upper bound + ma_order_list=[]; + for k=1:E_JDS + ma_order1=DownSampOrderEstimationR(data1(k,:)); + ma_order_list=[ma_order_list;ma_order1]; + end + dsd=median(ma_order_list); +end + +function p0 = order_estimate_complex_c( N, T0, dsd, S1 ) +% order estimation of real-valued multivariate process +% N: number of time points +% T0: original sample size +% dsd: downsampling depth +% S1: singular values +% p0: estimated order +% Program by Geng-Shen Fu: fugengs1@umbc.edu +T = T0/dsd; % effective sample size +lambda = S1.^2/T0; % eigenvalues + +DL = zeros(N,1); +for p = 0 : N-1 + + nf = 1+p*N-0.5*p*(p-1); % number of degrees of freedom + DL(p+1) = DL(p+1) + 0.5*nf*log(T)/T; + + for r = 1 : p + DL(p+1) = DL(p+1) + 0.5*log(lambda(r)); + end + + sigma2v = sum(lambda(p+1:N))/(N-p); + DL(p+1) = DL(p+1) + 0.5*(N-p)*log(sigma2v); + +end +[dl0, p0] = min(real(DL)); +p0 = p0-1; +%figure; plot([0:N-1], real(DL),'--o') + + + +function [ACS dsd ma_order_list] = Entropy_rate_estimator( data ) +N = size(data,1); +T = size(data, 2); +ma_order_list=[]; + +for i = 1:N + dsd=DownSampOrderEstimationR(data(i,:)); + ma_order_list = [ma_order_list;dsd]; +end; +% ma_order_list=[40 35 30 25 10 10 10 10 10 10]'; + +dsd=round(mean(ma_order_list)); +dsd_max=max(ma_order_list); + +% ma_order_list(ma_order_listinf + p0=inf; + return; + end + end +end + + +function s = isWhiteR(x) +% hypothesis: white process vs colored process +% if accept white, return 1; otherwise, return 0 +% Program by Xi-Lin Li: lixilin@umbc.edu +T=length(x); + +% DL of white case +p=0; +dl0=log(mean(abs(x).^2)); + +% DL of colored case with AR order 1 +p=1; +x1 = [0,x(1:end-1)]; +a = x*x1'/(x1*x1'); +n = x-a*x1; +dl1 = log(mean(abs(n).^2)) + p*log(T)/T/2; + +s=(dl0 0] = 1 + tedana_mask_img = nib.Nifti1Image(mask_array, tedana_mask_img.affine) + + # Save tedana optimally combined data and mask into mat files + tedana_optcom_mat = os.path.join(sbj_dir, "optcom_bold.mat") + tedana_mask_mat = os.path.join(sbj_dir, "mask.mat") + print("Saving tedana optimally combined data and mask into mat files") + scipy_io.savemat(tedana_optcom_mat, {"data": tedana_optcom_img.get_fdata()}) + scipy_io.savemat(tedana_mask_mat, {"mask": tedana_mask_img.get_fdata()}) + + # Run mapca + print("Running mapca") + pca = MovingAveragePCA(normalize=True) + _ = pca.fit_transform(tedana_optcom_img, tedana_mask_img) + + # Get AIC, KIC and MDL values + aic = pca.aic_ + kic = pca.kic_ + mdl = pca.mdl_ + + # Remove tedana output directory and the anat and func directories + subprocess.run(f"rm -rf {tedana_output_dir}", shell=True, cwd=repo) + subprocess.run(f"datalad drop {sbj}/anat", shell=True, cwd=repo) + subprocess.run(f"datalad drop {sbj}/func", shell=True, cwd=repo) + + # Here run matlab script with subprocess.run + print("Running GIFT version of maPCA") + + cmd = f'matlab -nodesktop -nosplash -nojvm -logfile {sbj_dir}/giftoutput.txt -r "addpath(genpath(\'{gift}\'));[comp_est_AIC,comp_est_KIC,comp_est_MDL,mdl,aic,kic]=icatb_estimate_dimension(\'{tedana_optcom_mat}\',\'{tedana_mask_mat}\',\'double\',3);save(\'{sbj_dir}/gift.mat\',\'comp_est_AIC\',\'comp_est_KIC\',\'comp_est_MDL\');quit"' + + proc = subprocess.Popen( + cmd, shell=True, stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE + ) + output, err = proc.communicate() + print(output.decode("utf-8")) + + giftmat = scipy_io.loadmat(os.path.join(sbj_dir, "gift.mat")) + + # Append AIC, KIC and MDL values to a pandas dataframe + print("Appending AIC, KIC and MDL values to csv file") + writer.writerow([sbj, + aic["n_components"], + giftmat["comp_est_AIC"][0][0], + kic["n_components"], + giftmat["comp_est_KIC"][0][0], + mdl["n_components"], + giftmat["comp_est_MDL"][0][0]]) + csv_file.flush() + print("Subject", sbj, "done") + + csv_file.close()