Skip to content

Commit

Permalink
GUI scaling
Browse files Browse the repository at this point in the history
  • Loading branch information
@Julie-Fabre
Julie-Fabre committed Sep 6, 2023
1 parent 66bfa27 commit b9d35fb
Show file tree
Hide file tree
Showing 4 changed files with 64 additions and 64 deletions.
116 changes: 58 additions & 58 deletions qualityMetrics/bc_fractionRPviolations.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
% ------
% theseSpikeTimes: nSpikesforThisUnit × 1 double vector of time in seconds
% of each of the unit's spikes.
% theseAmplitudes: nSpikesforThisUnit × 1 double vector of the amplitude scaling factor
% theseAmplitudes: nSpikesforThisUnit × 1 double vector of the amplitude scaling factor
% that was applied to the template when extracting that spike
% , only needed if plotThis is set to true
% tauR: refractory period
Expand All @@ -19,69 +19,69 @@
% ------
% fractionRPVs: fraction of contamination
% nRPVs: number refractory period violations
% overestimateBool: boolean, true if the number of refractory period violations
% overestimateBool: boolean, true if the number of refractory period violations
% is too high. we then overestimate the fraction of
% contamination.
% contamination.
% ------
% Reference
% Reference
% ------
% Hill, Daniel N., Samar B. Mehta, and David Kleinfeld.
% Hill, Daniel N., Samar B. Mehta, and David Kleinfeld.
% "Quality metrics to accompany spike sorting of extracellular signals."
% Journal of Neuroscience 31.24 (2011): 8699-8705:
% r = 2*(tauR - tauC) * N^2 * (1-Fp) * Fp / T , solve for Fp , fraction
% refractory period violatons. 2 factor because rogue spikes can occur before or
% after true spike

% initialize variables
fractionRPVs = nan(length(timeChunks)-1, length(tauR));
overestimateBool = nan(length(timeChunks)-1, length(tauR));
nRPVs = nan(length(timeChunks)-1, length(tauR));
fractionRPVs = nan(length(timeChunks)-1, length(tauR));
overestimateBool = nan(length(timeChunks)-1, length(tauR));
nRPVs = nan(length(timeChunks)-1, length(tauR));

if plotThis
figure('Color','none');
figure('Color', 'none');
end

for iTimeChunk = 1:length(timeChunks) - 1 %loop through each time chunk
% number of spikes in chunk
% number of spikes in chunk
N_chunk = length(theseSpikeTimes(theseSpikeTimes >= timeChunks(iTimeChunk) & theseSpikeTimes < timeChunks(iTimeChunk+1)));
% total times at which refractory period violations can occur
for iTauR_value = 1:length(tauR)
a = 2 * (tauR(iTauR_value) - tauC) * N_chunk^2 / abs(diff(timeChunks(iTimeChunk:iTimeChunk+1)));
% observed number of refractory period violations
nRPVs = sum(diff(theseSpikeTimes(theseSpikeTimes >= timeChunks(iTimeChunk) & theseSpikeTimes < timeChunks(iTimeChunk+1))) <= tauR(iTauR_value));
a = 2 * (tauR(iTauR_value) - tauC) * N_chunk^2 / abs(diff(timeChunks(iTimeChunk:iTimeChunk+1)));
% observed number of refractory period violations
nRPVs = sum(diff(theseSpikeTimes(theseSpikeTimes >= timeChunks(iTimeChunk) & theseSpikeTimes < timeChunks(iTimeChunk+1))) <= tauR(iTauR_value));

if nRPVs == 0 % no observed refractory period violations - this can
% also be because there are no spikes in this interval - use presence ratio to weed this out
fractionRPVs(iTimeChunk,iTauR_value) = 0;
overestimateBool(iTimeChunk,iTauR_value) = 0;
else % otherwise solve the equation above
rts = roots([-1, 1, -nRPVs / a]);
fractionRPVs(iTimeChunk,iTauR_value) = min(rts);
overestimateBool(iTimeChunk,iTauR_value) = 0;
if ~isreal(fractionRPVs(iTimeChunk,iTauR_value)) % function returns imaginary number if r is too high: overestimate number.
overestimateBool(iTimeChunk,iTauR_value) = 1;
if nRPVs < N_chunk %to not get a negative wierd number or a 0 denominator
fractionRPVs(iTimeChunk,iTauR_value) = nRPVs / (2 * (tauR(iTauR_value) - tauC) * (N_chunk - nRPVs));
else
fractionRPVs(iTimeChunk,iTauR_value) = 1;
if nRPVs == 0 % no observed refractory period violations - this can
% also be because there are no spikes in this interval - use presence ratio to weed this out
fractionRPVs(iTimeChunk, iTauR_value) = 0;
overestimateBool(iTimeChunk, iTauR_value) = 0;
else % otherwise solve the equation above
rts = roots([-1, 1, -nRPVs / a]);
fractionRPVs(iTimeChunk, iTauR_value) = min(rts);
overestimateBool(iTimeChunk, iTauR_value) = 0;
if ~isreal(fractionRPVs(iTimeChunk, iTauR_value)) % function returns imaginary number if r is too high: overestimate number.
overestimateBool(iTimeChunk, iTauR_value) = 1;
if nRPVs < N_chunk %to not get a negative wierd number or a 0 denominator
fractionRPVs(iTimeChunk, iTauR_value) = nRPVs / (2 * (tauR(iTauR_value) - tauC) * (N_chunk - nRPVs));
else
fractionRPVs(iTimeChunk, iTauR_value) = 1;
end
end
if fractionRPVs(iTimeChunk, iTauR_value) > 1 %it is nonsense to have a rate >1, the assumptions are failing here
fractionRPVs(iTimeChunk, iTauR_value) = 1;
end
end
if fractionRPVs(iTimeChunk,iTauR_value) > 1 %it is nonsense to have a rate >1, the assumptions are failing here
fractionRPVs(iTimeChunk,iTauR_value) = 1;
end
end
end




if plotThis
subplot(2, length(timeChunks) - 1, (length(timeChunks) - 1)+iTimeChunk)

subplot(2, length(timeChunks)-1, (length(timeChunks) - 1)+iTimeChunk)
theseISI = diff(theseSpikeTimes(theseSpikeTimes >= timeChunks(iTimeChunk) & theseSpikeTimes < timeChunks(iTimeChunk+1)));
theseisiclean = theseISI(theseISI >= tauC); % removed duplicate spikes
[isiProba, edgesISI] = histcounts(theseisiclean*1000, [0:0.5:50]);
bar(edgesISI(1:end-1)+mean(diff(edgesISI)), isiProba, 'FaceColor', [0, 0.35, 0.71], ...
'EdgeColor', [0, 0.35, 0.71]); %Check FR
if iTimeChunk ==1
'EdgeColor', [0, 0.35, 0.71]); %Check FR
if iTimeChunk == 1
xlabel('Interspike interval (ms)')
ylabel('# of spikes')
else
Expand All @@ -90,49 +90,49 @@
end
ylims = ylim;
[fr, ~] = histcounts(theseisiclean*1000, [0:0.5:5000]);
line([0, 10], [nanmean(fr(800:1000)), nanmean(fr(800:1000))], 'Color',[0.86, 0.2, 0.13], 'LineStyle', '--');
line([0, 10], [nanmean(fr(800:1000)), nanmean(fr(800:1000))], 'Color', [0.86, 0.2, 0.13], 'LineStyle', '--');
dummyh = line(nan, nan, 'Linestyle', 'none', 'Marker', 'none', 'Color', 'none');
if isnan(RPV_tauR_estimate)
for iTauR_value = [1,length(tauR)]
line([tauR(iTauR_value)*1000, tauR(iTauR_value)*1000], [ylims(1), ylims(2)], 'Color', [0.86, 0.2, 0.13]);
for iTauR_value = [1, length(tauR)]
line([tauR(iTauR_value) * 1000, tauR(iTauR_value) * 1000], [ylims(1), ylims(2)], 'Color', [0.86, 0.2, 0.13]);
end
if length(tauR) == 1
legend([dummyh], {[num2str(round(fractionRPVs(iTimeChunk,1)*100,1)), '% rpv']},...
'Location', 'NorthEast','TextColor', [0.7, 0.7, 0.7], 'Color', 'none');
legend([dummyh], {[num2str(round(fractionRPVs(iTimeChunk, 1)*100, 1)), '% rpv']}, ...
'Location', 'NorthEast', 'TextColor', [0.7, 0.7, 0.7], 'Color', 'none');
else
legend([dummyh, dummyh], {[num2str(round(fractionRPVs(iTimeChunk,1)*100,1)), '% rpv'],...
[num2str(round(fractionRPVs(iTimeChunk,length(tauR))*100,1)), '% rpv']},'Location', 'NorthEast','TextColor', [0.7, 0.7, 0.7], 'Color', 'none');
legend([dummyh, dummyh], {[num2str(round(fractionRPVs(iTimeChunk, 1)*100, 1)), '% rpv'], ...
[num2str(round(fractionRPVs(iTimeChunk, length(tauR))*100, 1)), '% rpv']}, 'Location', 'NorthEast', 'TextColor', [0.7, 0.7, 0.7], 'Color', 'none');
end

else
line([tauR(RPV_tauR_estimate)*1000, tauR(RPV_tauR_estimate)*1000], [ylims(1), ylims(2)], 'Color', [0.86, 0.2, 0.13]);
legend([dummyh], {[num2str(fractionRPVs(RPV_tauR_estimate)*100,1), '% rpv']},'Location', 'NorthEast','TextColor', [0.7, 0.7, 0.7], 'Color', 'none');
line([tauR(RPV_tauR_estimate) * 1000, tauR(RPV_tauR_estimate) * 1000], [ylims(1), ylims(2)], 'Color', [0.86, 0.2, 0.13]);
legend([dummyh], {[num2str(fractionRPVs(RPV_tauR_estimate)*100, 1), '% rpv']}, 'Location', 'NorthEast', 'TextColor', [0.7, 0.7, 0.7], 'Color', 'none');

end

set(gca, 'XScale', 'log')
end

end
if plotThis
subplot(2,numel(timeChunks)-1, 1:numel(timeChunks)-1)
scatter(theseSpikeTimes, theseAmplitudes, 4,[0, 0.35, 0.71],'filled'); hold on;
% chunk lines
subplot(2, numel(timeChunks)-1, 1:numel(timeChunks)-1)
scatter(theseSpikeTimes, theseAmplitudes, 4, [0, 0.35, 0.71], 'filled');
hold on;
% chunk lines
ylims = ylim;
for iTimeChunk = 1:length(timeChunks)
line([timeChunks(iTimeChunk),timeChunks(iTimeChunk)],...
[ylims(1),ylims(2)], 'Color', [0.7, 0.7, 0.7])
line([timeChunks(iTimeChunk), timeChunks(iTimeChunk)], ...
[ylims(1), ylims(2)], 'Color', [0.7, 0.7, 0.7])
end
xlabel('time (s)')
ylabel(['amplitude scaling' newline 'factor'])
ylabel(['amplitude scaling', newline, 'factor'])
if exist('prettify_plot', 'file')
prettify_plot('none','none','none')
prettify_plot('none', 'none', 'none')
else
warning('https://github.com/Julie-Fabre/prettify-matlab repo missing - download it and add it to your matlab path to make plots pretty')
makepretty('none')
end
end

end


end
2 changes: 1 addition & 1 deletion qualityMetrics/bc_plotGlobalQualityMetric.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -240,6 +240,6 @@ function bc_plotGlobalQualityMetric(qMetric, param, unitType, uniqueTemplates, t
prettify_plot('none','none','w')
else
warning('https://github.com/Julie-Fabre/prettify-matlab repo missing - download it and add it to your matlab path to make plots pretty')
makepretty('w')
% makepretty('w')
end
end
2 changes: 1 addition & 1 deletion qualityMetrics/bc_runAllQualityMetrics.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,5 +231,5 @@
end

unitType = bc_getQualityUnitType(param, qMetric);
bc_plotGlobalQualityMetric(qMetric, param, unitType, uniqueTemplates, forGUI.tempWv);
%bc_plotGlobalQualityMetric(qMetric, param, unitType, uniqueTemplates, forGUI.tempWv);
end
8 changes: 4 additions & 4 deletions visualizationTools/bc_unitQualityGUI.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -330,24 +330,24 @@ function updateUnit(unitQualityGuiHandle, memMapData, ephysData, rawWaveforms, i
vals(iChanToPlot) = max(abs(squeeze(ephysData.templates(thisUnit, :, chansToPlot(iChanToPlot)))));
if maxChan == chansToPlot(iChanToPlot)
set(guiData.maxTemplateWaveformLines, 'XData', (ephysData.waveform_t + (ephysData.channel_positions(chansToPlot(iChanToPlot), 1) - 11) / 10), ...
'YData', -squeeze(ephysData.templates(thisUnit, :, chansToPlot(iChanToPlot)))'+ ...
'YData', -squeeze(ephysData.templates(thisUnit, :, chansToPlot(iChanToPlot)))'./50+ ...
(ephysData.channel_positions(chansToPlot(iChanToPlot), 2) ./ 100));
hold on;
set(guiData.peaks, 'XData', (ephysData.waveform_t(forGUI.peakLocs{iCluster}) ...
+(ephysData.channel_positions(chansToPlot(iChanToPlot), 1) - 11) / 10), ...
'YData', -squeeze(ephysData.templates(thisUnit, forGUI.peakLocs{iCluster}, chansToPlot(iChanToPlot)))'+ ...
'YData', -squeeze(ephysData.templates(thisUnit, forGUI.peakLocs{iCluster}, chansToPlot(iChanToPlot)))'./50+ ...
(ephysData.channel_positions(chansToPlot(iChanToPlot), 2) ./ 100));

set(guiData.troughs, 'XData', (ephysData.waveform_t(forGUI.troughLocs{iCluster}) ...
+(ephysData.channel_positions(chansToPlot(iChanToPlot), 1) - 11) / 10), ...
'YData', -squeeze(ephysData.templates(thisUnit, forGUI.troughLocs{iCluster}, chansToPlot(iChanToPlot)))'+ ...
'YData', -squeeze(ephysData.templates(thisUnit, forGUI.troughLocs{iCluster}, chansToPlot(iChanToPlot)))'./50+ ...
(ephysData.channel_positions(chansToPlot(iChanToPlot), 2) ./ 100));
set(guiData.templateWaveformLines(iChanToPlot), 'XData', nan(82, 1), ...
'YData', nan(82, 1));

else
set(guiData.templateWaveformLines(iChanToPlot), 'XData', (ephysData.waveform_t + (ephysData.channel_positions(chansToPlot(iChanToPlot), 1) - 11) / 10), ...
'YData', -squeeze(ephysData.templates(thisUnit, :, chansToPlot(iChanToPlot)))'+ ...
'YData', -squeeze(ephysData.templates(thisUnit, :, chansToPlot(iChanToPlot)))'./50+ ...
(ephysData.channel_positions(chansToPlot(iChanToPlot), 2) ./ 100));
end
end
Expand Down

0 comments on commit b9d35fb

Please sign in to comment.