/
plot_spike_density_multiunit.m
executable file
·256 lines (216 loc) · 7.82 KB
/
plot_spike_density_multiunit.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
function handles = plot_spike_density_multiunit (parsedData, parsedParams, varargin)
%% Plots a spike density plot from parsed multiunit data
% Usage: handles = plot_spike_density_multiunit (parsedData, parsedParams, varargin)
% Explanation:
% TODO
%
% Example(s):
% spikeDensityHzTrace1 = compute_spike_density([4, 4, 4, 5, 8, 9, 9, 10]);
% spikeDensityHzTrace2 = compute_spike_density([6, 6, 6, 7, 7, 7]);
% spikeDensityHz = {spikeDensityHzTrace1, spikeDensityHzTrace2};
% parsedData = table(spikeDensityHz);
% siSeconds = [1; 1];
% minTimeSec = [0; 0.1];
% stimStartSec = [3; 3.1];
% maxTimeSec = [20; 20.1];
% phaseNumber = [1; 2];
% figPathBase = {'test_plot_spike_density_multiunit_trace1'; 'test_plot_spike_density_multiunit_trace2'};
% parsedParams = table(siSeconds, minTimeSec, stimStartSec, maxTimeSec, phaseNumber, figPathBase);
% handles = plot_spike_density_multiunit (parsedData, parsedParams);
%
% Outputs:
% handles - a structure with fields:
% im
% boundaryLine
% stimStartLine
% specified as a scalar structure
%
% Arguments:
% parsedData - parsed data
% must be a table with fields:
% spikeDensityHz
% parsedParams- parsed parameters
% must be a table with fields:
% siSeconds
% minTimeSec
% maxTimeSec
% stimStartSec
% phaseNumber
% figPathBase
% varargin - 'FigPosition':
% must be a TODO
% default == TODO
% - 'XLimits':
% must be a TODO
% default == TODO
% - 'PlotStimStart':
% must be a TODO
% default == TODO
% - 'BoundaryType':
% must be a TODO
% default == TODO
% - 'MaxNYTicks':
% must be a TODO
% default == TODO
% - Any other parameter-value pair for imagesc()
%
% Requires:
% cd/compute_index_boundaries.m
% cd/create_indices.m
% cd/create_error_for_nargin.m
% cd/create_labels_from_numbers.m
% cd/decide_on_colormap.m
% cd/extract_common_prefix.m
% cd/force_matrix.m
% cd/hold_off.m
% cd/hold_on.m
% cd/plot_vertical_line.m
% cd/plot_window_boundaries.m
% cd/struct2arglist.m
%
% Used by:
% cd/m3ha_oscillations_analyze.m
% cd/parse_multiunit.m
% File History:
% 2019-11-30 Moved from parse_multiunit.m
% TODO: Pull out code as function plot_heat_map(spikeDensityHz);
%% Hard-coded parameters
validBoundaryTypes = {'horizontalLines', 'verticalBars'};
% TODO: Make optional argument
colorMap = [];
% Default values for optional arguments
figPositionDefault = [];
xLimitsDefault = [];
plotStimStartDefault = true;
boundaryTypeDefault = 'horizontalLines';
maxNYTicksDefault = 20;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Deal with arguments
% Check number of required arguments
if nargin < 2
error(create_error_for_nargin(mfilename));
end
% Set up Input Parser Scheme
iP = inputParser;
iP.FunctionName = mfilename;
iP.KeepUnmatched = true; % allow extraneous options
% Add required inputs to an input Parser
addRequired(iP, 'parsedData', ...
@(x) validateattributes(x, {'table'}, {'2d'}));
addRequired(iP, 'parsedParams', ...
@(x) validateattributes(x, {'table'}, {'2d'}));
% Add parameter-value pairs to the Input Parser
addParameter(iP, 'FigPosition', figPositionDefault, ...
@(x) validateattributes(x, {'numeric'}, {'2d'}));
addParameter(iP, 'XLimits', xLimitsDefault, ...
@(x) validateattributes(x, {'numeric'}, {'2d'}));
addParameter(iP, 'PlotStimStart', plotStimStartDefault, ...
@(x) validateattributes(x, {'logical', 'numeric'}, {'binary'}));
addParameter(iP, 'BoundaryType', boundaryTypeDefault, ...
@(x) any(validatestring(x, validBoundaryTypes)));
addParameter(iP, 'MaxNYTicks', maxNYTicksDefault, ...
@(x) validateattributes(x, {'numeric'}, {'2d'}));
% Read from the Input Parser
parse(iP, parsedData, parsedParams, varargin{:});
figPosition = iP.Results.FigPosition;
xLimits = iP.Results.XLimits;
plotStimStart = iP.Results.PlotStimStart;
boundaryType = validatestring(iP.Results.BoundaryType, validBoundaryTypes);
maxNYTicks = iP.Results.MaxNYTicks;
% Keep unmatched arguments for the imagesc() function
otherArguments = struct2arglist(iP.Unmatched);
%% Preparation
% Retrieve data for plotting
spikeDensityHz = parsedData.spikeDensityHz;
siSeconds = parsedParams.siSeconds;
minTimeSec = parsedParams.minTimeSec;
maxTimeSec = parsedParams.maxTimeSec;
stimStartSec = parsedParams.stimStartSec;
phaseVector = parsedParams.phaseNumber;
figPathBase = parsedParams.figPathBase;
% Extract the original file base
fileBase = extract_common_prefix(figPathBase);
% Create a figure title base
titleBase = replace(fileBase, '_', '\_');
% Compute the phase boundaries
phaseBoundaries = compute_index_boundaries('Grouping', phaseVector, ...
'TreatNaNsAsGroup', false);
%% Plot as a heatmap
% Hold on
wasHold = hold_on;
% Compute stimulation start
meanStimStartSec = mean(stimStartSec);
% Count traces
nSweeps = numel(spikeDensityHz);
% Get the average sampling interval in seconds
siSeconds = mean(siSeconds);
% Set x and y end points
xEnds = [min(minTimeSec); max(maxTimeSec)];
yEnds = [1; nSweeps];
% Set x and y limits
if isempty(xLimits)
xLimits = [xEnds(1) - 0.5 * siSeconds; xEnds(2) + 0.5 * siSeconds];
end
yLimits = [yEnds(1) - 0.5; yEnds(2) + 0.5];
% Compute bar X position if it is used
barXPosition = mean([xLimits(1), meanStimStartSec]);
% Decide on y ticks and labels
yTicks = create_indices('IndexEnd', nSweeps, 'MaxNum', maxNYTicks, ...
'AlignMethod', 'left');
yTickLabels = create_labels_from_numbers(nSweeps - yTicks + 1);
% Force as a matrix and transpose it so that
% each trace is a row
spikeDensityMatrix = transpose(force_matrix(spikeDensityHz));
% Set a gray-scale color map
% colormap(flipud(gray));
% colormap(jet);
cm = decide_on_colormap(colorMap, 'ColorMapFunc', @gray, ...
'ReverseOrder', true, 'HighContrast', true);
colormap(cm);
% Generate plot
im = imagesc(xEnds, flipud(yEnds), spikeDensityMatrix, otherArguments{:});
% Set the y ticks and labels
yticks(yTicks);
yticklabels(yTickLabels);
% Plot stimulation start
if plotStimStart
stimStartLine = plot_vertical_line(meanStimStartSec, 'Color', rgb('Green'), ...
'LineStyle', '--', 'LineWidth', 0.5, ...
'YLimits', yLimits);
end
% Plot phase boundaries
if ~isempty(phaseBoundaries)
% Compute y values for phase boundaries
yBoundaries = nSweeps - phaseBoundaries + 1;
% Plot phase boundaries
boundaryLine = ...
plot_window_boundaries(yBoundaries, 'BoundaryType', boundaryType, ...
'BarValue', barXPosition);
end
% Update x limits and labels
xlim(xLimits);
ylim(yLimits);
xlabel('Time (s)');
ylabel('Trace #');
title(['Spike density (Hz) for ', titleBase]);
% Show a scale bar
colorbar;
% Set figure position if requested
if ~isempty(figPosition)
set(gcf, 'Position', figPosition);
end
% Hold off
hold_off(wasHold);
%% Save output
handles.im = im;
if exist('boundaryLine', 'var')
handles.boundaryLine = boundaryLine;
end
if exist('stimStartLine', 'var')
handles.stimStartLine = stimStartLine;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%{
OLD CODE:
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%