Home > chronux_1_50 > spikesort > ss_energy.m

ss_energy

PURPOSE ^

SS_ENERGY Interface energy based cluster similarity computation.

SYNOPSIS ^

function spikes = ss_energy(spikes)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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);

Generated on Mon 09-Oct-2006 00:54:52 by m2html © 2003