Home > chronux_2_00 > spikesort > graphics > aggtree.m

aggtree

PURPOSE ^

AGGTREE Generates a graphical summary of an aggregation procedure.

SYNOPSIS ^

function aggtree(spikes)

DESCRIPTION ^

  AGGTREE  Generates a graphical summary of an aggregation procedure.
      AGGTREE(spikes)
 Takes a spike-sorting data structure that has undergone hierarchical
   aggregation and visualizes the hierarchy.
 The function draws its output over the current axis.  The labels of the
   initial clusters are shown at the bottom of the axis and bracket lines
   indicate cluster merging.  The height of a bracket indicates the degree
   of overlap between the two daughter clusters; short brackets imply high
   overlap.  The bracket color evolves over the 'winter' colormap, with blue
   representing the first aggregation step and green representing the last.
   Finally, the number of spikes and isi statistic for each final cluster
   are given at the top of the plot.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function aggtree(spikes)
0002 
0003 %  AGGTREE  Generates a graphical summary of an aggregation procedure.
0004 %      AGGTREE(spikes)
0005 % Takes a spike-sorting data structure that has undergone hierarchical
0006 %   aggregation and visualizes the hierarchy.
0007 % The function draws its output over the current axis.  The labels of the
0008 %   initial clusters are shown at the bottom of the axis and bracket lines
0009 %   indicate cluster merging.  The height of a bracket indicates the degree
0010 %   of overlap between the two daughter clusters; short brackets imply high
0011 %   overlap.  The bracket color evolves over the 'winter' colormap, with blue
0012 %   representing the first aggregation step and green representing the last.
0013 %   Finally, the number of spikes and isi statistic for each final cluster
0014 %   are given at the top of the plot.
0015 
0016 if (~isfield(spikes, 'hierarchy'))
0017     error('The input data structure to AGGTREE does not have a hierarchy defined.');
0018 end
0019 tree = spikes.hierarchy.tree;
0020 assignments = spikes.hierarchy.assigns;
0021 times = spikes.spiketimes;
0022 
0023 % Get some numbers set up . . .
0024 numorigclusts = length(unique(tree(:,[1,2])));  % note that maxorigclusts can be
0025 maxorigclusts = max(unique(tree(:,[1,2])));     %    higher if there are empty clusters
0026 aggsteps = size(tree,1);
0027 numnodes = maxorigclusts + aggsteps;
0028 
0029 % . . . and make some room to store tree info.
0030 node = repmat(struct('parent', 0, 'lchild', 0, 'rchild', 0), [numnodes,1]);
0031 xpos = zeros(numnodes,1);
0032 ypos = zeros(numnodes,1);
0033 
0034 % modify tree . . .
0035 atree = zeros(aggsteps,4);
0036 atree(:,[1,2]) = tree(:,[1,2]);   % original cluster names
0037 atree(:,4) = 1.1*max(tree(:,3)) - tree(:,3);  % small when connect_strength is big
0038 
0039 % Give every node a unique name; the 'aggregate' function output reuses label
0040 %  names, but we'll want a separate identifier for each.
0041 for step = 1:aggsteps
0042     oldname = atree(step,1);
0043     newname = maxorigclusts + step;
0044     atree(step,3) = newname;
0045     
0046     atree(find(atree(step+1:end, 1) == oldname) + step, 1) = newname;
0047     atree(find(atree(step+1:end, 2) == oldname) + step, 2) = newname;
0048 end
0049 
0050 % Build tree structure by assigning left/right child and parent indices
0051 % to each aggregation node, determining node y locations along the way.
0052 for step = 1:aggsteps
0053     lchild = atree(step,1);
0054     rchild = atree(step,2);
0055     parent = atree(step,3);
0056  
0057     % y position is the greater of the parents positions + an amount
0058     % derived from connection strength (computed above as 'atree(:,4)')
0059     ypos(parent) = max(ypos([lchild,rchild])) + atree(step,4);
0060 
0061     % store parent/child indices
0062     node(parent).lchild = lchild;
0063     node(parent).rchild = rchild;
0064     node(lchild).parent = parent;
0065     node(rchild).parent = parent;
0066 end
0067 
0068 % Final clusters are those with no parents (i.e., tree roots).
0069 treeroots = find(cat(2,node.parent) == 0);
0070 
0071 % Assign node x locations to leafs (i.e., original clusters labels)
0072 % using a (LIFO) stack to walk the tree left to right.
0073 nextleaf = 1;
0074 stack = treeroots;
0075 while (~isempty(stack))
0076     current = stack(1);
0077     stack = stack(2:end);
0078     if (node(current).lchild ~= 0)  % if not a leaf, push children (left first)
0079         stack = [node(current).lchild, node(current).rchild, stack];
0080     else   % if we're at a leaf, take the next available label
0081         xpos(current) = nextleaf;
0082         nextleaf = nextleaf + 1;
0083     end
0084 end
0085 
0086 % For interior (i.e., non-leaf nodes), x pos is just average of children's xpos
0087 for this = 1:numnodes
0088     if (xpos(this) == 0)  % only interior nodes are unassigned
0089         xpos(this) = mean([xpos(node(this).lchild), xpos(node(this).rchild)]);
0090     end
0091 end
0092 
0093 % Make labels for the leafs (easy because they're ordered by their
0094 % original assignment labels -- so we just write these down and resort
0095 % based on y position to isolate leaves and x position to take the leaves
0096 % from left to right.)
0097 leaflabels = sortrows([(1:numnodes)' xpos ypos], [3,2]); % sort y then x
0098 leaflabels(ypos ~= 0,:) = [];    % leafs have y pos == 0
0099 leaflabels = num2str(leaflabels(:,1));
0100 
0101 % Draw dots at the nodes . . .
0102 cla
0103 plot(xpos,ypos,'.','MarkerFaceColor',[0 0 1]);
0104 hold on;
0105 
0106 % Now the tree itself; lines with colors indicating aggregation order
0107 cmap = winter(256);
0108 cind = floor(linspace(1,256,aggsteps));  % make use of the entire colormap
0109 for step = 1:aggsteps  % draw the brackets for each step
0110     x1 = xpos(atree(step,1));
0111     y1 = ypos(atree(step,1));
0112     x2 = xpos(atree(step,2));
0113     y2 = ypos(atree(step,2));
0114     y3 = ypos(atree(step,3));
0115     line([x1 x1 x2 x2], [y1 y3 y3 y2], 'Color', cmap(cind(step),:));
0116 end
0117 
0118 % Finally, draw extra lines to highlight the final clusters and
0119 % display info for them.
0120 tmin = size(spikes.waveforms,2)./spikes.Fs;
0121 if (isfield(spikes, 'options') && isfield(spikes.options, 'refractory_period'))
0122     tref = spikes.options.refractory_period;
0123 else   tref = max(0.002, tmin*1.5);
0124 end
0125 yscale = max(max(ypos)*1.5, 0.1);
0126 for root = 1:length(treeroots)
0127     x1 = xpos(treeroots(root));
0128     y1 = ypos(treeroots(root));
0129     y2 = max(ypos) * 1.1;  % place final cluster node higher than the rest
0130     
0131     % figure out the original label that this root matches
0132     oldname = tree((atree(:,3) == treeroots(root)),1);
0133     if (isempty(oldname))
0134         oldname = treeroots(root);
0135     end
0136 
0137     % get cluster size and timing information for this final cluster
0138     members = find(assignments == oldname);
0139     clustsize = length(members);
0140     membertimes = times(members);
0141     [a, scores] = isiQuality(membertimes, membertimes, tmin, 0.010, tref, spikes.Fs);
0142     
0143     plot(x1,y2,'o','MarkerFaceColor',[0 0 0]);
0144     line([x1 x1],[y1 
0145         y2], 'LineWidth', 2, 'Color', [0 0 0]);
0146     stagger = (yscale/20) * rem(root-1,3);
0147     text(x1 - 0.61, y2*1.05 + stagger,    ['Clust #' num2str(oldname)], 'Color', [0 0.5 0.5], 'FontSize', 8);
0148     text(x1 - 0.59, y2*1.03 + stagger,    ['N=' num2str(clustsize)], 'Color', [1 0 0], 'FontSize', 8);
0149     text(x1 + 0.22, y2*0.90 + stagger/3,  ['isi=' num2str(scores(1),2)], 'Color', [0 0 1], 'FontSize', 8);
0150 end
0151 
0152 % Prettify the axes
0153 hold off;
0154 colormap(cmap);
0155 %if (length(leaflabels) <= 32)
0156     set(gca, 'XTickLabel', leaflabels);
0157     set(gca, 'XTick', (1:size(leaflabels,1)), 'XTickLabelMode', 'manual');
0158 %else
0159 %    set(gca, 'XTickLabel', '');
0160 %end
0161 set(gca,'Xlim',[0 length(leaflabels)+1]);
0162 set(gca,'Ylim',[0 yscale]);

Generated on Fri 15-Aug-2008 11:35:42 by m2html © 2003