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

Generated on Sun 13-Aug-2006 11:49:44 by m2html © 2003