


SS_ENERGY Interface energy based cluster similarity computation.
SPIKES = SS_ENERGY(SPIKES) adds an interface energy matrix to a
spike-sorting object in SPIKES.HIERARCHY.INTERFACE_ENERGY.
The energy similarity matrix is calculated by applying an exponential
decay to all pairwise euclidean distances between waveforms from two
clusters (or within a single cluster for intra-cluster energy) and
summing these distances.
The calculation ignores the energy due to the zero distance between
points and themselves; this removes a dependence of the density on
the absolute size of the cluster. As a result, singleton clusters
do not have a well-defined energy and will cause an error.
When each entry is normalized by the number of distinct contributing
pairs (Na*Nb for off diagonal entries and Na*(Na-1)/2 on the diagonal),
it approximates the fraction of pairs in a given cluster whose distance
is not much greater than the length constant of the exponential and thus
provides an estimate of local density. This function does not, however,
normalize SPIKES.HIERARCHY.INTERFACE_ENERGY, since the normalized form is
inconvenient during cluster aggregation. The normalization can readily
be done, however, with
normalize = ((numpts * numpts') - diag(numpts));
normalize = normalize - diag(0.5*diag(normalize));
normalized_energy = interface_energy ./ normalize;
where 'numpts' is a vector of cluster sizes.
The unnormalized energy matrix can be updated during aggregation without
the need to recompute it from scratch. The intra-cluster energy E(AB,AB)
of a cluster AB formed by aggregating clusters A and B is given by
E(AB,AB) = E(A,A) + E(B,B) + E(A,B)
and the inter-cluster energy between any cluster C and an aggregate AB is
E(AB,C) = E(A,C) + E(B,C)
References:
Fee MS et al (1996). J. Neurosci Methods (69): 175-88

0001 function spikes = ss_energy(spikes) 0002 % SS_ENERGY Interface energy based cluster similarity computation. 0003 % SPIKES = SS_ENERGY(SPIKES) adds an interface energy matrix to a 0004 % spike-sorting object in SPIKES.HIERARCHY.INTERFACE_ENERGY. 0005 % 0006 % The energy similarity matrix is calculated by applying an exponential 0007 % decay to all pairwise euclidean distances between waveforms from two 0008 % clusters (or within a single cluster for intra-cluster energy) and 0009 % summing these distances. 0010 % 0011 % The calculation ignores the energy due to the zero distance between 0012 % points and themselves; this removes a dependence of the density on 0013 % the absolute size of the cluster. As a result, singleton clusters 0014 % do not have a well-defined energy and will cause an error. 0015 % 0016 % When each entry is normalized by the number of distinct contributing 0017 % pairs (Na*Nb for off diagonal entries and Na*(Na-1)/2 on the diagonal), 0018 % it approximates the fraction of pairs in a given cluster whose distance 0019 % is not much greater than the length constant of the exponential and thus 0020 % provides an estimate of local density. This function does not, however, 0021 % normalize SPIKES.HIERARCHY.INTERFACE_ENERGY, since the normalized form is 0022 % inconvenient during cluster aggregation. The normalization can readily 0023 % be done, however, with 0024 % normalize = ((numpts * numpts') - diag(numpts)); 0025 % normalize = normalize - diag(0.5*diag(normalize)); 0026 % normalized_energy = interface_energy ./ normalize; 0027 % where 'numpts' is a vector of cluster sizes. 0028 % 0029 % The unnormalized energy matrix can be updated during aggregation without 0030 % the need to recompute it from scratch. The intra-cluster energy E(AB,AB) 0031 % of a cluster AB formed by aggregating clusters A and B is given by 0032 % E(AB,AB) = E(A,A) + E(B,B) + E(A,B) 0033 % and the inter-cluster energy between any cluster C and an aggregate AB is 0034 % E(AB,C) = E(A,C) + E(B,C) 0035 % 0036 % References: 0037 % Fee MS et al (1996). J. Neurosci Methods (69): 175-88 0038 0039 % Last Modified By: sbm on Fri Oct 7 21:35:16 2005 0040 0041 starttime = clock; 0042 0043 %%%%% ARGUMENT CHECKING 0044 if (~isfield(spikes, 'waveforms') | (size(spikes.waveforms, 1) < 1)) 0045 error('SS:waveforms_undefined', 'The SS object does not contain any waveforms!'); 0046 elseif (~isfield(spikes, 'overcluster')) 0047 error('SS:overcluster_not_computed', 'The data must be overclustered before computing energy'); 0048 end 0049 numclusts = length(unique(spikes.overcluster.assigns)); 0050 waves = spikes.waveforms; 0051 0052 %%%%% PREPARE SOME INFORMATION 0053 normsqr = sum(waves.^2,2); 0054 pts = cell(numclusts,1); % collect spike indices for each cluster 0055 for clust = 1:numclusts 0056 pts{clust} = find(spikes.overcluster.assigns == clust); 0057 end 0058 numpts = cellfun('length', pts); 0059 if (any(numpts < 2)) 0060 error('SS:energy_ill_defined', 'Clusters with fewer than 2 points do not have a defined energy.'); 0061 end 0062 0063 %%%%% HEURISTIC DISTANCE SCALE that seems to work. The calculation is not too sensitive to this parameter. 0064 scale = sqrt(sum(diag(spikes.overcluster.W))) / 10; 0065 0066 %%%%% PREPARE TO LOOP 0067 total = (numclusts^2 + numclusts) / 2; 0068 k = 1; 0069 progressBar(0, max(floor(total/100),1), 'Computing Interaction Energies . . .') 0070 interface_energy = zeros(numclusts); 0071 0072 %%%%% PAIRWISE DISTANCES LOOP 0073 assigns = spikes.overcluster.assigns; 0074 for clust1 = 1:numclusts 0075 % Deal with self-case first, because PAIRDIST works better with a 0076 % different syntax in this case. 0077 dists = pairdist(waves(pts{clust1},:), 'reuse'); 0078 interface_energy(clust1,clust1) = fast_interface_energy(dists,scale); 0079 k = k + 1; 0080 for clust2 = (clust1+1):numclusts % now for the rest ... 0081 dists = pairdist(waves(pts{clust1},:), waves(pts{clust2},:), 'reuse'); 0082 interface_energy(clust1,clust2) = fast_interface_energy(dists,scale); 0083 k = k + 1; 0084 progressBar(k/total); 0085 end 0086 end 0087 0088 %%%%% CORRECTION TERMS 0089 % The energy matrix so far includes a contribution in the intra-cluster 0090 % energies that is not found in the inter-cluster energies; namely, the 0091 % computation of sum_(all x) sum_(all y) e^(-dist/scale) for 0092 % intra-cluster energy includes cases where x == y (so dist == 0). 0093 interface_energy = interface_energy - diag(numpts); % So subtract this out. 0094 0095 % Also, we've double counted pairs in the intra-energy case, since dist(a,b) 0096 % and dist(b,a) are not treated as distinct; 0097 interface_energy = interface_energy - diag(0.5*diag(interface_energy)); 0098 0099 %%%%% FINISH UP 0100 spikes.hierarchy.interface_energy = interface_energy; 0101 spikes.tictoc.energy = etime(clock, starttime);