0001 function aggtree(spikes)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
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
0026 numorigclusts = length(unique(tree(:,[1,2])));
0027 maxorigclusts = max(unique(tree(:,[1,2])));
0028 aggsteps = size(tree,1);
0029 numnodes = maxorigclusts + aggsteps;
0030
0031
0032 node = repmat(struct('parent', 0, 'lchild', 0, 'rchild', 0), [numnodes,1]);
0033 xpos = zeros(numnodes,1);
0034 ypos = zeros(numnodes,1);
0035
0036
0037 atree = zeros(aggsteps,4);
0038 atree(:,[1,2]) = tree(:,[1,2]);
0039 atree(:,4) = 1.1*max(tree(:,3)) - tree(:,3);
0040
0041
0042
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
0053
0054 for step = 1:aggsteps
0055 lchild = atree(step,1);
0056 rchild = atree(step,2);
0057 parent = atree(step,3);
0058
0059
0060
0061 ypos(parent) = max(ypos([lchild,rchild])) + atree(step,4);
0062
0063
0064 node(parent).lchild = lchild;
0065 node(parent).rchild = rchild;
0066 node(lchild).parent = parent;
0067 node(rchild).parent = parent;
0068 end
0069
0070
0071 treeroots = find(cat(2,node.parent) == 0);
0072
0073
0074
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)
0081 stack = [node(current).lchild, node(current).rchild, stack];
0082 else
0083 xpos(current) = nextleaf;
0084 nextleaf = nextleaf + 1;
0085 end
0086 end
0087
0088
0089 for this = 1:numnodes
0090 if (xpos(this) == 0)
0091 xpos(this) = mean([xpos(node(this).lchild), xpos(node(this).rchild)]);
0092 end
0093 end
0094
0095
0096
0097
0098
0099 leaflabels = sortrows([[1:numnodes]' xpos ypos], [3,2]);
0100 leaflabels(find(ypos ~= 0),:) = [];
0101 leaflabels = num2str(leaflabels(:,1));
0102
0103
0104 cla
0105 plot(xpos,ypos,'.','MarkerFaceColor',[0 0 1]);
0106 hold on;
0107
0108
0109 cmap = winter(256);
0110 cind = floor(linspace(1,256,aggsteps));
0111 for step = 1:aggsteps
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
0121
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;
0132
0133
0134 oldname = tree((find(atree(:,3) == treeroots(root))),1);
0135 if (isempty(oldname))
0136 oldname = treeroots(root);
0137 end
0138
0139
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
0155 hold off;
0156 colormap(cmap);
0157
0158 set(gca, 'XTickLabel', leaflabels);
0159 set(gca, 'XTick', [1:size(leaflabels,1)], 'XTickLabelMode', 'manual');
0160
0161
0162
0163 set(gca,'Xlim',[0 length(leaflabels)+1]);
0164 set(gca,'Ylim',[0 yscale]);