function cval=cosmo_measure_clusters(sample,nbrhood_mat,cluster_stat,varargin)
% General cluster measure
%
% cval=cosmo_measure_clusters(sample,nbrhood_mat,method,['dh',dh|'threshold',t],...)
%
% Inputs:
% sample 1xP data array for k features
% nbrhood_mat MxP neighborhood matrix, if each feature has at
% most M neighbors. nbrhood_mx(:,k) should contain
% the feature ids of the neighbors of feature k
% (values of zero indicate no neighbors)
% cluster_stat 'tfce', 'max', 'maxsize', or 'maxsum'.
% 'dh',dh when method is 'tfce', the integral step size
% 'threshold',t if method is anything but 'tfce', the threshold
% level
% 'E', E (optional) when method is 'tfce': the extent
% exponent, default E=0.5
% 'H', E (optional) when method is 'tfce': the height
% exponent, default H=2
% 'feature_sizes', s (optional) 1xP element size of each feature;
% default is a vector with all elements equal to 1.
%
% Output:
% cval 1xP data array with cluster values.
% Define:
% cl_k(thr) Feature ids that are connected
% with feature k when sample is
% thresholded at thr.
% (with cl_k(thr)=[] if sample(k)<thr).
% e_k(thr) The extent of the cluster for feature
% k when thresholded at thr, i.e.
% e_k=sum[i=cl_k(thr)] s(i)
% max0[x=xs]v(x) the maximum value that v takes for
% any value in xs, or zero if xs is
% empty.
% The output for a feature k depends on cluster_stat:
% - 'tfce':
% cval(k)=sum[h=dh:dh:max(sample)] e_k(h)^E *h^H * dh
% i.e. a weighted sum of the extent and height of
% the cluster at different thresholds
% - 'max':
% cval(k)=max0[i=cl_k(thr)] samples(i))
% i.e. the maximum value in the cluster that contains k
% - 'maxsize'
% cval(k)=e_k(thr)
% i.e. the size of the cluster that contains k
% - 'maxsum'
% cval(k)=sum[i=cl_k(thr)] samples(i)
% i.e. the sum of the values in the cluster that
% contains k
%
% Examples:
% % very simple one-dimensional 'line' example
% samples= [1 2 1 1 0 2 2];
% nbrhood_mat= [1 1 2 3 4 5 6;...
% 2 2 3 4 5 6 7;...
% 0 3 4 5 6 7 0];
% % illustrate TFCE cluster stat. Note that the results vary little
% % as a function of dh, and that both features in clusters with a
% % larger extent but less extreme value and features in clusters with
% % smaller extent but more extreme values have larger values
% s=cosmo_measure_clusters(samples,nbrhood_mat,'tfce','dh',.1);
% cosmo_disp(s,'threshold',inf)
% %|| [ 0.77 2.86 0.77 0.77 0 3.49 3.49 ]
% %
% s=cosmo_measure_clusters(samples,nbrhood_mat,'tfce','dh',.05);
% cosmo_disp(s,'threshold',inf)
% %|| [ 0.618 2.88 0.618 0.618 0 3.63 3.63 ]
% %
% % illustrate other cluster stats. Note that the results varies more
% % as a function of the chosen threshold
% cosmo_measure_clusters(samples,nbrhood_mat,'max','threshold',1)
% %|| 2 2 2 2 0 2 2
% %
% cosmo_measure_clusters(samples,nbrhood_mat,'max','threshold',2)
% %|| 0 2 0 0 0 2 2
% %
% cosmo_measure_clusters(samples,nbrhood_mat,'maxsize','threshold',1)
% %|| 4 4 4 4 0 2 2
% %
% cosmo_measure_clusters(samples,nbrhood_mat,'maxsize','threshold',2)
% %|| 0 1 0 0 0 2 2
% %
% cosmo_measure_clusters(samples,nbrhood_mat,'maxsum','threshold',1)
% %|| 5 5 5 5 0 4 4
% %
% cosmo_measure_clusters(samples,nbrhood_mat,'maxsum','threshold',2)
% %|| 0 2 0 0 0 4 4
%
% Notes:
% - the 'method' argument is similar to FieldTrip's clusterstat method,
% but does support 'tfce' and does not support 'wcm'
% - unlike FieldTrip's clusterstat, a value is returned for each feature
% (rather than each cluster).
% - TFCE is advised in the general case, because it finds a compromise
% between magnitude of values and extent of clusters.
% - this function is used by cosmo_montecarlo_cluster_stat for
% non-parametric cluster-based correction for multiple comparisons.
%
% References:
% - Stephen M. Smith, Thomas E. Nichols (2009), Threshold-free
% cluster enhancement: Addressing problems of smoothing, threshold
% dependence and localisation in cluster inference, NeuroImage, Volume
% 44, 83-98.
% - Maris, E., Oostenveld, R. Nonparametric statistical testing of EEG-
% and MEG-data. Journal of Neuroscience Methods (2007).
%
% See also: cosmo_convert_neighborhood, cosmo_montecarlo_cluster_stat
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #