IsolationDistance

IsolationDistance - Determine the isolation quality for one or more clusters.

 Compute isolation distance, as defined by K. Harris (see Schmitzer-Torbert
 et al. 2005).

  USAGE

    [d,units] = IsolationDistance(features,<options>)

    features       features obtained using e.g. <a href="matlab:help GetSpikeFeatures">GetSpikeFeatures</a> (or a
                   subset thereof)
    <options>      optional list of property-value pairs (see table below)

    =========================================================================
     Properties    Values
    -------------------------------------------------------------------------
     'units'       list of units to process (default = all clusters except 0
                   and 1)
     'show'        plot Mahalanobis distance PDF and CDF (default = 'off')
    =========================================================================

  OUTPUT
 
    d              For each cluster, a measure of the separation between the
                   cluster and all other spikes. Namely, the Mahalanobis
                   distance of the Nth closest spike outside the cluster
                   (where N is the number of spikes inside the cluster).
    units          corresponding (group,cluster) pairs

0001 function [distance,units] = IsolationDistance(features,varargin)
0002 
0003 %IsolationDistance - Determine the isolation quality for one or more clusters.
0004 %
0005 % Compute isolation distance, as defined by K. Harris (see Schmitzer-Torbert
0006 % et al. 2005).
0007 %
0008 %  USAGE
0009 %
0010 %    [d,units] = IsolationDistance(features,<options>)
0011 %
0012 %    features       features obtained using e.g. <a href="matlab:help GetSpikeFeatures">GetSpikeFeatures</a> (or a
0013 %                   subset thereof)
0014 %    <options>      optional list of property-value pairs (see table below)
0015 %
0016 %    =========================================================================
0017 %     Properties    Values
0018 %    -------------------------------------------------------------------------
0019 %     'units'       list of units to process (default = all clusters except 0
0020 %                   and 1)
0021 %     'show'        plot Mahalanobis distance PDF and CDF (default = 'off')
0022 %    =========================================================================
0023 %
0024 %  OUTPUT
0025 %
0026 %    d              For each cluster, a measure of the separation between the
0027 %                   cluster and all other spikes. Namely, the Mahalanobis
0028 %                   distance of the Nth closest spike outside the cluster
0029 %                   (where N is the number of spikes inside the cluster).
0030 %    units          corresponding (group,cluster) pairs
0031 %
0032 
0033 % Copyright (C) 2012 by Michaël Zugaro
0034 %
0035 % This program is free software; you can redistribute it and/or modify
0036 % it under the terms of the GNU General Public License as published by
0037 % the Free Software Foundation; either version 3 of the License, or
0038 % (at your option) any later version.
0039 
0040 % Defaults
0041 distance = [];
0042 nFeatures = size(features,2)-3;
0043 units = [];
0044 show = 'off';
0045 
0046 % Check number of parameters
0047 if nargin < 1,
0048   error('Incorrect number of parameters (type ''help <a href="matlab:help IsolationDistance">IsolationDistance</a>'' for details).');
0049 end
0050 
0051 % Check parameter sizes
0052 if length(size(features)) ~= 2 | size(features,2) < 4,
0053     error('Incorrect parameter ''features'' (type ''help <a href="matlab:help IsolationDistance">IsolationDistance</a>'' for details).');
0054 end
0055 
0056 % Parse options
0057 for i = 1:2:length(varargin),
0058     if ~ischar(varargin{i}),
0059         error(['Parameter ' num2str(i) ' is not a property (type ''help <a href="matlab:help IsolationDistance">IsolationDistance</a>'' for details).']);
0060     end
0061     switch(lower(varargin{i})),
0062         case 'units',
0063             units = varargin{i+1};
0064             if ~isimatrix(units,'>=0') | size(units,2) ~= 2,
0065                 error('Incorrect value for property ''units'' (type ''help <a href="matlab:help IsolationDistance">IsolationDistance</a>'' for details).');
0066             end
0067         case 'show',
0068             show = varargin{i+1};
0069             if ~isastring(show,'on','off'),
0070                 error('Incorrect value for property ''show'' (type ''help <a href="matlab:help IsolationDistance">IsolationDistance</a>'' for details).');
0071             end
0072         otherwise,
0073             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help IsolationDistance">IsolationDistance</a>'' for details).']);
0074     end
0075 end
0076 
0077 % Discard noise and artefacts
0078 artefacts = features(:,3)<=1;
0079 features(artefacts,:) = [];
0080 
0081 % By default, process all clusters
0082 if isempty(units),
0083     units = features(:,[2 3]);
0084     units = unique(units,'rows');
0085     units = units(units(:,2)~=0&units(:,2)~=1,:);
0086 end
0087 
0088 if strcmp(lower(show),'on'),
0089     figure;
0090 end
0091 for i = 1:size(units,1),
0092     unit = units(i,:);
0093     distance(i,1) = ComputeOnce(features,unit,show);
0094 end
0095 
0096 % --------------------------------------------------------------------------------------------------------
0097 
0098 function distance = ComputeOnce(features,unit,show)
0099 
0100 distance = NaN;
0101 nFeatures = size(features,2)-3;
0102 unitStr = ['(' int2str(unit(1)) ',' int2str(unit(2)) ')'];
0103 
0104 % Determine spikes inside vs outside cluster
0105 spikesInCluster = features(:,2)==unit(1) & features(:,3)==unit(2);
0106 nSpikesInCluster = sum(spikesInCluster);
0107 spikesInElectrodeGroup = features(:,2)==unit(1);
0108 nSpikesInElectrodeGroup = sum(spikesInElectrodeGroup);
0109 
0110 % We need at least as many spikes as features, but not more than half the total number of spikes
0111 if nSpikesInCluster < nFeatures,
0112     warning(['Not enough spikes for unit ' unitStr '.']);
0113     return
0114 end
0115 if nSpikesInCluster > nSpikesInElectrodeGroup/2,
0116     warning(['Too many spikes for unit ' unitStr '.']);
0117     return
0118 end
0119 
0120 % Compute Mahalanobis distance for spikes inside (mIn) and outside (mOut) the cluster
0121 m = mahal(features(:,4:end),features(spikesInCluster,4:end));
0122 mIn = m(spikesInCluster);
0123 mOut = m(~spikesInCluster);
0124 % Determine the Mahalanobis of the Nth closest spike outside the cluster (where N is the number of spikes inside the cluster)
0125 sOut = sort(mOut);
0126 distance = sOut(nSpikesInCluster);
0127 
0128 if strcmp(lower(show),'on'),
0129     % Compute and plot pdfs
0130     figure;
0131     bins = 0:1:300;
0132     subplot(3,1,1);
0133     pdfIn = hist(mIn,bins);
0134     pdfOut = hist(mOut,bins);
0135     plot(bins,[pdfIn',pdfOut']);
0136     xlim([0 100]);
0137     xlabel('Mahalanobis distance (PDF)');
0138     ylabel('Count');
0139     legend(['Cluster ' unitStr],'Others','location','southeast');
0140 
0141     % Compute and plot cdfs
0142     subplot(3,1,2);
0143     cdfIn = cumsum(pdfIn);
0144     cdfOut = cumsum(pdfOut);
0145     plot(bins,[cdfIn',cdfOut']);
0146     xlim([0 100]);
0147     xlabel('Mahalanobis distance (CDF)');
0148     ylabel('Count');
0149     legend(['Cluster ' unitStr],'Others','location','southeast');
0150 
0151     % Compute and plot ordered distances
0152     subplot(3,1,3);
0153     sIn = sort(mIn);
0154     n = length(sIn);
0155     plot(1:n,[sIn sOut(1:n)]);
0156     hold on;
0157     plot([1 n],[distance distance],'k:');
0158     xlim([1 n*1.1]);
0159     xlabel('Spike #');
0160     ylabel('Mahalanobis distance');
0161     t = text(n,distance,num2str(round(100*distance)/100));
0162     set(t,'fontsize',12);
0163     drawnow;
0164 end