diff --git a/.gitignore b/.gitignore index 17583d6..c65de63 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ getHost.m .* *~ ~* +startup*.* diff --git a/aodgui/ScanInspector.fig b/aodgui/ScanInspector.fig index f796eae..06752e7 100644 Binary files a/aodgui/ScanInspector.fig and b/aodgui/ScanInspector.fig differ diff --git a/aodgui/ScanInspector.m b/aodgui/ScanInspector.m index 7b7a227..9a98a5f 100644 --- a/aodgui/ScanInspector.m +++ b/aodgui/ScanInspector.m @@ -163,6 +163,19 @@ function AodScans_Callback(hObject, eventdata, handles) set(handles.Preprocessing, 'String', 'None'); set(handles.Preprocessing, 'Value', 1); end +cla(handles.Traces) + +if count(aod.ScanMotion & scan) == 1 + motion = fetch(aod.ScanMotion & scan, '*') + plot(handles.Motion, motion.t,motion.x+3,motion.t,motion.y,motion.t,motion.z-3); + linkaxes([handles.Traces, handles.Motion], 'x'); + xlim(motion.t([1 end])); + ylim([-6 6]); + legend('x','y','z'); +else + cla(handles.Motion) +end + % --- Executes during object creation, after setting all properties. diff --git a/james_setPath.m b/james_setPath.m index 68eaaa7..1c19538 100644 --- a/james_setPath.m +++ b/james_setPath.m @@ -2,8 +2,13 @@ warning off MATLAB:dispatcher:nameConflict -% user specific DJ connection parameters (uses Alex' credentials) -host = 'at-database.neusc.bcm.tmc.edu'; +% user specific DJ connection parameters +local = false +if local + host = '127.0.0.1:8889' % port is for MAMP server +else + host = 'at-database.ad.bcm.edu'; +end user = 'jcotton'; setenv('DJ_HOST', host) setenv('DJ_USER', user) @@ -27,14 +32,16 @@ % DataJoint library is assumed to be in the same directory as the base % diretory addpath(fullfile(base, '../DataJoint/matlab')); +%addpath('z:\users\alex\projects\datajoint\') +%addpath('Z:\users\alex\projects\mym\distribution\mexw64') addpath(fullfile(base, '../moksm')); if isequal(computer, 'PCWIN64') - addpath(fullfile(base, '../DataJoint/mym/win64')); + addpath(fullfile(base, '../mym/win64')); else - addpath(fullfile(base, '../DataJoint/mym')); + addpath(fullfile(base, '../mym/distribution/mexmaci64/')); end @@ -49,8 +56,9 @@ addpath(getLocalPath('/lab/users/james/Matlab/VariationalClustering')); % spike detection - run(getLocalPath('/lab/users/james/Matlab/spikesorting/detection/setPath')); - + %run(getLocalPath('/lab/users/james/Matlab/spikesorting/detection/setPath')); + run C:\Users\jcotton\Documents\MATLAB\spikedetection\setPath.m + warning on MATLAB:dispatcher:nameConflict else warning('at-lab not mounded. File access tools not added to path'); diff --git a/processing/RawPathMap.m b/processing/RawPathMap.m index b97162d..ab151e6 100644 --- a/processing/RawPathMap.m +++ b/processing/RawPathMap.m @@ -13,6 +13,8 @@ methods function self = RawPathMap(varargin) + self.search = self.getSearchPath(); + self.temp = self.getTempPath(); end function outFile = findFile(self, file) @@ -22,7 +24,7 @@ % the search path (map.search) is used, then the temporary % location (map.temp), and finally at_scratch. assert(~isempty(regexp(file, '^(/|\\)raw', 'once')), 'Not a raw data file: %s', file) - loc = [self.search, {getLocalPath('/at_scratch')}]; + loc = [self.search, {getLocalPath('/at_scratch')}, '/Volumes/M$']; for i = 1:numel(loc) outFile = strrep([loc{i}, file(5:end)], '\', '/'); if numel(dir(sprintf(strrep(outFile, filesep, '/'), 0))) % sprintf in case of HDF5 family file @@ -40,8 +42,8 @@ end function file = toTemp(self, file) - % Map /raw to the temporary location (typically on local drive) - % file = toTemp(RawPathMap, file) + % Map /raw to the temporary location (typically on local drive) + % file = toTemp(RawPathMap, file) assert(~isempty(regexp(file, '^(/|\\)raw', 'once')), 'Not a raw data file: %s', file) file = strrep([self.temp file(5:end)], '\', '/'); end @@ -52,7 +54,7 @@ if exist('rawmap.mat', 'file') search = getfield(load('rawmap'), 'search'); %#ok else - search = {'M:', 'N:', 'O:', 'P:', getLocalPath('/at_scratch')}; + search = {'M:', 'N:', 'O:', 'P:', 'F:', getLocalPath('/at_scratch'), getLocalPath('/raw')}; end end diff --git a/processing/cleanup.m b/processing/cleanup.m index 8d0cbaf..94970ae 100644 --- a/processing/cleanup.m +++ b/processing/cleanup.m @@ -23,8 +23,8 @@ function cleanup(sessKey) end % create stop_time entries if missing - aodStopTime = fetch1(aod, 'aod_scan_stop_time'); - if isnan(aodStopTime) || ~aodStopTime + stop_time_null_count = count(aod & 'aod_scan_stop_time is null') + if stop_time_null_count == 1 br = getFile(aod,'Temporal'); duration = 1000 * length(br) / getSamplingRate(br); close(br); @@ -75,7 +75,7 @@ function cleanup(sessKey) % create stop_time entries if missing behStopTime = fetch1(beh, 'beh_stop_time'); - if isnan(behStopTime) || ~behStopTime + if isnan(behStopTime) br = getFile(beh); duration = 1000 * length(br) / getSamplingRate(br); close(br); @@ -98,8 +98,8 @@ function cleanup(sessKey) % create stop_time entries if missing stimStopTime = fetch1(stimulation, 'stim_stop_time'); - if isnan(stimStopTime) || ~stimStopTime - stim = getfield(load(stimFile), 'stim'); + if true || isnan(stimStopTime) + stim = getfield(load(stimFile), 'stim'); %#ok if isempty(stim.events) stimStopTime = key.stim_start_time; totalTrials = 0; diff --git a/processing/processSet.m b/processing/processSet.m index 322d450..1d87d3e 100755 --- a/processing/processSet.m +++ b/processing/processSet.m @@ -147,12 +147,13 @@ function processSet(key, spikesCb, spikesFile, lfpCb, muaCb, pathCb, useTempDir) delete(fullfilefs(tempDir, dataFilePattern)); delete(fullfilefs(tempDir, '*.h5')); + localSpikesDir = fullfilefs(localProcessedDir, spikesDir); % copy output files to processed folder - if exist(localProcessedDir, 'file') - rmdir(localProcessedDir, 's'); + if exist(localSpikesDir, 'file') + rmdir(localSpikesDir, 's'); end - mkdir(localProcessedDir); - copyfile(destDir, localProcessedDir); + mkdir(localSpikesDir); + copyfile(outDir, localSpikesDir); % delete temp data try diff --git a/processing/sync/detectLcdPhotodiodeFlips.m b/processing/sync/detectLcdPhotodiodeFlips.m index bb94cad..e14638a 100644 --- a/processing/sync/detectLcdPhotodiodeFlips.m +++ b/processing/sync/detectLcdPhotodiodeFlips.m @@ -137,3 +137,30 @@ peakSamples = yp + ndx - p - 1 + b - 1; peakTimes = [peakTimes; reshape(br(peakSamples,'t'), [], 1)]; %#ok end + + +% spaced_max -- find local maxima separated by a given minimum interval +% +% Syntax: +% idx = spaced_max(x, min_interval) +% +% DY: 2010-08-16 + +function idx = spaced_max(x, min_interval, thresh) +peaks = local_max( x ); +if nargin>2 + peaks = peaks(x(peaks)>thresh); +end +if isempty(peaks) + idx = []; +else + idx=peaks(1); + for i=peaks(2:end)' + if i-idx(end)>=min_interval + idx(end+1)=i; %#ok + elseif x(i)>x(idx(end)) + idx(end)=i; + end + end +end + diff --git a/processing/sync/syncAod.m b/processing/sync/syncAod.m index 170ed97..a40baa2 100644 --- a/processing/sync/syncAod.m +++ b/processing/sync/syncAod.m @@ -4,9 +4,10 @@ params.oldFile = false; params.maxPhotodiodeErr = 3.00; % 100 us err allowed -params.behDiodeOffset = [2 17]; % [min max] in ms +params.behDiodeOffset = [-40 40]; % [min max] in ms params.behDiodeSlopeErr = 1e-5; % max deviation from 1 params.diodeThreshold = 0.04; +params.minFracSwaps = 0.99; % should find photodiode for almost all mac swap times params.minNegTime = -100; % 100 ms timing error if isempty(stim.events) @@ -35,8 +36,8 @@ % arithmetic and then smooth the result to account for jitter in the % swaptimes Fs = 10; % kHz -k = 20000; % max offset (samples) in each direction -smooth = 10; % smoothing window for finding the peak (half-width); +k = 500; % max offset (samples) in each direction +smooth = 40; % smoothing window for finding the peak (half-width); c = zeros(2 * k + 1, 1); for i = -k:k c(i + k + 1) = isectq(round(macSwapTimes * Fs + i), round(diodeSwapTimes * Fs)); @@ -51,6 +52,7 @@ % throw out swaps that don't have matches within one ms originalDiodeSwapTimes = diodeSwapTimes; +originalMacSwapTimes = macSwapTimes; [macSwapTimes, diodeSwapTimes] = matchTimes(macSwapTimes, diodeSwapTimes, offset); N = numel(macSwapTimes); @@ -64,8 +66,11 @@ t = mean(macSwapTimes); shift = macPar(2) * t - t + macPar(1); +fracSwaps = numel(macSwapTimes) / numel(originalMacSwapTimes); + if(~(abs(macPar(2) - 1) < params.behDiodeSlopeErr ... - && shift > params.behDiodeOffset(1) && shift < params.behDiodeOffset(2))) + && shift > params.behDiodeOffset(1) && shift < params.behDiodeOffset(2) ... + && fracSwaps > params.minFracSwaps)) originalMacTimes = cat(1, stim.params.trials.swapTimes) + offset; @@ -80,7 +85,7 @@ testMacSwapTimes = cat(1, stimDiodeTest.params.trials.swapTimes); testDiodeSwapTimes = originalDiodeSwapTimes; - plot(testMacSwapTimes, zeros(size(testMacSwapTimes))-0.5,'.',testDiodeSwapTimes, zeros(size(testDiodeSwapTimes))+0.5,'.'); + plot(testMacSwapTimes, zeros(size(testMacSwapTimes))+0.5,'.',testDiodeSwapTimes, zeros(size(testDiodeSwapTimes))-0.5,'.'); ylim([-5 5]); title('Synced times'); @@ -91,28 +96,32 @@ linkaxes(h,'x'); - disp('Sync failed'); + fprintf('Sync failed. Note only %f of mac swap times were kept', fracSwaps); keyboard end assert(abs(macPar(2) - 1) < params.behDiodeSlopeErr ... && shift > params.behDiodeOffset(1) && shift < params.behDiodeOffset(2), ... 'Regression between behavior clock and photodiode clock outside system tolerances'); - +assert(fracSwaps > params.minFracSwaps, 'Too many swap times could not be aligned') % convert times in stim file stimDiode = convertStimTimes(stim, macPar, [0 1]); stimDiode.synchronized = 'diode'; % plot residuals figure +subplot(211) macSwapTimes = cat(1, stimDiode.params.trials.swapTimes); diodeSwapTimes = originalDiodeSwapTimes; [macSwapTimes, diodeSwapTimes] = matchTimes(macSwapTimes, diodeSwapTimes, 0); %assert(N == numel(macSwapTimes), 'Error during timestamp conversion. Number of timestamps don''t match!') res = macSwapTimes(:) - diodeSwapTimes(:); plot(diodeSwapTimes, res, '.k'); -rms = sqrt(mean(res.^2)); +subplot(212) +plot(macSwapTimes,zeros(size(macSwapTimes))+0.5,'.',diodeSwapTimes,zeros(size(diodeSwapTimes))-0.5,'.') +ylim([-5 5]) +rms = sqrt(mean(res.^2)); if(rms > params.maxPhotodiodeErr) disp('Residuals too large'); keyboard diff --git a/processing/sync/syncEphys.m b/processing/sync/syncEphys.m index 31574d2..0fb83f2 100644 --- a/processing/sync/syncEphys.m +++ b/processing/sync/syncEphys.m @@ -4,7 +4,16 @@ params.oldFile = false; params.maxPhotodiodeErr = 0.5; % 0.5 ms err allowed -params.behDiodeOffset = [2.5 4]; % [min max] in ms + +% CAUTION: this min has been modified from 2.5 to 2.0 on +% June 8, 2017 by Edgar Y. Walker encountering cases of offset ~2.4951 +% reliably when processing +% subject_id: 34 +% setup: 1 +% session_start_time: 3579798609122 +% stim_start_time: 3579798640685 +% +params.behDiodeOffset = [2.0 4]; % [min max] in ms params.behDiodeSlopeErr = 1e-6; % max deviation from 1 params.diodeThreshold = 0.04; params.minNegTime = -100; % 100 ms timing error @@ -18,14 +27,18 @@ assert(strcmp(stim.synchronized, 'network'), 'Run network sync first!') -% Get photodiode swap times +% Get photodiode swap times from Mac Psychtoolbox time +% Mac event times are (most likely) relative to program start time +% swap times are stored in milli seconds tstart = stim.params.trials(1).swapTimes(1) - 500; tend = stim.params.trials(end).swapTimes(end) + 500; + +% reading the photodiode signal from Ephys and detect photodiode peak times br = getFile(acq.Ephys(key), 'Photodiode'); [peakTimes, peakAmps] = detectPhotodiodePeaks(br, tstart, tend); close(br); -% detect swaps +% detect swaps from the photodiode peak times da = abs(diff(peakAmps)); [mu, v] = MoG1(da(:), 2, 'cycles', 50, 'mu', [0 median(da)]); sd = sqrt(v); @@ -69,7 +82,7 @@ diodeOffset = macPar(1) + macSwapTimes(1) * diodeSlopeErr; assert(abs(diodeSlopeErr) < params.behDiodeSlopeErr ... && diodeOffset > params.behDiodeOffset(1) && diodeOffset < params.behDiodeOffset(2), ... - 'Regression between behavior clock and photodiode clock outside system tolerances'); + 'Regression between behavior clock and photodiode clock outside system tolerances: diodeSlopeErr = %f and diodeOffset = %f', diodeSlopeErr, diodeOffset); % convert times in stim file stimDiode = convertStimTimes(stim, macPar, [0 1]); diff --git a/processing/sync/syncNetwork.m b/processing/sync/syncNetwork.m index 1678c5c..d642290 100644 --- a/processing/sync/syncNetwork.m +++ b/processing/sync/syncNetwork.m @@ -17,7 +17,7 @@ % to convert from mac time to pc counter use p(2) * mac + p(1) i = (e - s) < params.maxRoundtrip; -assert(sum(i) / numel(i) > 0.8, 'Too many sync packets dropped') +assert(sum(i) / numel(i) > 0.6, 'Too many sync packets dropped') p = myrobustfit(mid(i), r(i)); rms = sqrt(mean((p(2) * mid + p(1) - r).^2)); diff --git a/processing/sync/syncNetworkFewTrials.m b/processing/sync/syncNetworkFewTrials.m new file mode 100644 index 0000000..04b15f5 --- /dev/null +++ b/processing/sync/syncNetworkFewTrials.m @@ -0,0 +1,32 @@ +function [stimNet, rms] = syncNetworkFewTrials(stim, key) + +params.maxRoundtrip = 5; + +if isempty(stim.events) + stimNet = stim; + rms = -1; + return % empty file; nothing to do +end + +sy = [stim.params.trials.sync]; +s = [sy.start]; +e = [sy.end]; +r = [sy.response]; +mid = (s + e) / 2; % these are the mac times +mid = mid / 1e3; % sync times are in ms + +% to convert from mac time to pc counter use p(2) * mac + p(1) +i = (e - s) < params.maxRoundtrip; +assert(sum(i) / numel(i) > 0.8, 'Too many sync packets dropped') +p = flipud(regress(r'-mean(r),[mid'-mean(mid) ones(length(mid),1)])); +p(1) = -mean(mid)*p(2) + mean(r); +rms = sqrt(mean((p(2) * mid + p(1) - r).^2)); + +% get the "zero" time of the counter that was used for network +% sync relative to session time +t0 = getHardwareStartTime(acq.BehaviorTraces(key)); + +assert(~isempty(t0), ['Could not find the hardware time stamp for behavioral session: ' fetch1(acq.BehaviorTraces(key),'beh_path')]) +stimNet = convertStimTimes(stim, p, [0; 1]); +stimNet = convertStimTimes(stimNet, [t0; 1], [t0; 1]); +stimNet.synchronized = 'network'; diff --git a/recovery/recoverFromFileSystem.m b/recovery/recoverFromFileSystem.m index 203b180..2fd2cce 100644 --- a/recovery/recoverFromFileSystem.m +++ b/recovery/recoverFromFileSystem.m @@ -18,7 +18,10 @@ function recoverFromFileSystem(folder) % mapping from subject ids to setup and experimenter setups = [2 2 2 3 3 2 3 3]; experimenters = {'Allison', 'Allison', 'Allison', 'Mani', 'Tori', 'Mani', 'Tori', 'Tori'}; -ephysTasks = {'UtahArray', 'UtahArray', 'UtahArray', 'Charles Chronic Left Tetrodes', 'Claude Chronic Left Tetrodes', 'UtahArray', 'SiliconProbes', 'SiliconProbes'}; +ephysTasks = {'UtahArray', 'UtahArray', 'UtahArray', 'Charles Chronic Left Tetrodes', 'Claude Chronic Left Tetrodes', 'UtahArray', 'SiliconProbes', 'SiliconProbes', 'Fuckyourface'}; + +setups(16) = 4; +experimenters(16) = {'Manolis'}; % read timestamps tsFileId = fopen(tsFile, 'r'); @@ -51,7 +54,8 @@ function recoverFromFileSystem(folder) % sessionStartTime = timestamperTime{1} - 1; sessionStartTime = toLabview(sessionDatetime); subjectName = tsFile(splitNdx(end-2)+1:splitNdx(end-1)-1); -subjectId = fetch1(acq.Subjects(struct('subject_name', subjectName)), 'subject_id'); + +subjectId = fetch1(acq.Subjects(struct('subject_name', subjectName)), 'subject_id') setup = setups(subjectId); sessKey.subject_id = subjectId; diff --git a/schemas/+acq/AodVolumeIgnore.m b/schemas/+acq/AodVolumeIgnore.m index acb32fe..efdcb41 100644 --- a/schemas/+acq/AodVolumeIgnore.m +++ b/schemas/+acq/AodVolumeIgnore.m @@ -19,8 +19,9 @@ end methods (Static) - function findShortVolumes - volumes = fetch(acq.AodVolume - acq.AodVolumeIgnore); + function findShortVolumes(key) + if nargin < 1; key = struct; end + volumes = fetch(acq.AodVolume(key) - acq.AodVolumeIgnore); for i = 1:length(volumes) try vol = getFile(acq.AodVolume & volumes(i)); diff --git a/schemas/+acq/SessionsCleanup.m b/schemas/+acq/SessionsCleanup.m index 00231b2..aecfef9 100644 --- a/schemas/+acq/SessionsCleanup.m +++ b/schemas/+acq/SessionsCleanup.m @@ -20,6 +20,7 @@ methods (Access=protected) function makeTuples(self, key) if count(acq.Sessions(key) & 'recording_software = "Acquisition2.0"') + key cleanup(key); end insert(self, key); diff --git a/schemas/+acq/StimulationSync.m b/schemas/+acq/StimulationSync.m index ec19cb1..2aeb701 100644 --- a/schemas/+acq/StimulationSync.m +++ b/schemas/+acq/StimulationSync.m @@ -84,7 +84,7 @@ function makeTuples(self, key) else % nothing to sync against % network sync if numel(stim.params.trials) <= 10 - disp 'Few trials' + disp 'Few trials and no recording. Not synchronizing.' return; end diff --git a/schemas/+aod/AodJobs.m b/schemas/+aod/AodJobs.m index a207fac..53f6d03 100644 --- a/schemas/+aod/AodJobs.m +++ b/schemas/+aod/AodJobs.m @@ -25,4 +25,4 @@ function run(self) end end end -%} \ No newline at end of file +%} diff --git a/schemas/+aod/CosineTuning.m b/schemas/+aod/CosineTuning.m index 5c199bd..b8477ef 100644 --- a/schemas/+aod/CosineTuning.m +++ b/schemas/+aod/CosineTuning.m @@ -10,7 +10,7 @@ r2 : double # The rsquared %} -classdef CosineTuning < dj.Relvar & dj.Automatic +classdef CosineTuning < dj.Relvar & dj.AutoPopulate properties(Constant) table = dj.Table('aod.CosineTuning') diff --git a/schemas/+aod/DesignMatrix.m b/schemas/+aod/DesignMatrix.m index caa4361..de65dd6 100644 --- a/schemas/+aod/DesignMatrix.m +++ b/schemas/+aod/DesignMatrix.m @@ -14,7 +14,7 @@ patched_cell_num : int # The cell number of the patched cell %} -classdef DesignMatrix < dj.Relvar & dj.Automatic +classdef DesignMatrix < dj.Relvar & dj.AutoPopulate properties(Constant) table = dj.Table('aod.DesignMatrix') @@ -28,7 +28,7 @@ self.restrict(varargin) end - function traces = removeStimulus(self) + function [traces regressors] = removeStimulus(self) % Remove any component that can be predicted from teh stimulus % from the traces @@ -37,7 +37,7 @@ x = [ones(size(dat.stimuli_matrix,1),1) dat.stimuli_matrix]; traces = zeros(size(dat.traces_matrix)); for i = 1:size(dat.traces_matrix,2); - [~,~,resid] = regress(dat.traces_matrix(:,i), x); + [regressors(:,i),~,resid] = regress(dat.traces_matrix(:,i), x); traces(:,i) = resid; end end @@ -54,9 +54,6 @@ function makeTuples(self, key) % Create a matrix of the traces [traces, tuple.cell_nums] = fetchn(aod.TracePreprocess & aod.UniqueCell & key, 'trace', 'cell_num'); tuple.traces_matrix = cat(2,traces{:}); - % Normalize the energy - tuple.traces_matrix = bsxfun(@rdivide, tuple.traces_matrix, ... - std(tuple.traces_matrix,[],1)); % Get the stimulus design matrix [tuple.stimuli_matrix startTime endTime] = ... @@ -68,6 +65,10 @@ function makeTuples(self, key) tuple.traces_matrix(delBins,:) = []; tuple.stimuli_matrix(delBins,:) = []; + % Normalize the energy + tuple.traces_matrix = bsxfun(@rdivide, tuple.traces_matrix, ... + std(tuple.traces_matrix,[],1)); + if count(aod.Spikes & key) == 1 spikes = fetch1(aod.Spikes & key, 'times'); tuple.spikes_binned = histc(spikes, tuple.times); diff --git a/schemas/+aod/EventEnergy.m b/schemas/+aod/EventEnergy.m index 05bcdda..b94ba20 100644 --- a/schemas/+aod/EventEnergy.m +++ b/schemas/+aod/EventEnergy.m @@ -10,7 +10,7 @@ iqr_ratio : double # p-val versus shuffled data %} -classdef EventEnergy < dj.Relvar & dj.Automatic +classdef EventEnergy < dj.Relvar & dj.AutoPopulate properties(Constant) table = dj.Table('aod.EventEnergy') diff --git a/schemas/+aod/ForwardTuning.m b/schemas/+aod/ForwardTuning.m index 2efe762..e8e4db3 100644 --- a/schemas/+aod/ForwardTuning.m +++ b/schemas/+aod/ForwardTuning.m @@ -82,7 +82,7 @@ function makeTuples( this, key ) for i = 1:length(trials) trial_info = fetch1(stimulation.StimTrials(trials(i)), 'trial_params'); event = fetch(stimulation.StimTrialEvents(trials(i), 'event_type="showSubStimulus"'),'*'); - onsets = sort([event.event_time]); + onsets = double(sort([event.event_time])); for j = 1:length(onsets) cond = trial_info.conditions(j); onset = onsets(j); diff --git a/schemas/+aod/ForwardVonMises.m b/schemas/+aod/ForwardVonMises.m index 7469a43..0abd6ef 100644 --- a/schemas/+aod/ForwardVonMises.m +++ b/schemas/+aod/ForwardVonMises.m @@ -12,7 +12,7 @@ von_base : double # dF/F base %} -classdef ForwardVonMises < dj.Relvar & dj.Automatic +classdef ForwardVonMises < dj.Relvar & dj.AutoPopulate properties(Constant) table = dj.Table('aod.ForwardVonMises'); popRel = aod.ForwardTuning; @@ -57,4 +57,4 @@ function makeTuples( this, key ) insert(aod.ForwardVonMises, tuple); end end -end \ No newline at end of file +end diff --git a/schemas/+aod/MotionKalman.m b/schemas/+aod/MotionKalman.m index befdf43..27237c9 100644 --- a/schemas/+aod/MotionKalman.m +++ b/schemas/+aod/MotionKalman.m @@ -13,7 +13,7 @@ resid : longblob # The residual traces %} -classdef MotionKalman < dj.Relvar & dj.Automatic +classdef MotionKalman < dj.Relvar & dj.AutoPopulate properties(Constant) table = dj.Table('aod.MotionKalman') diff --git a/schemas/+aod/MotionScore.m b/schemas/+aod/MotionScore.m index 796566e..1b2ce00 100644 --- a/schemas/+aod/MotionScore.m +++ b/schemas/+aod/MotionScore.m @@ -10,7 +10,7 @@ trace_corr : double # PC1 and motion correlation %} -classdef MotionScore < dj.Relvar & dj.Automatic +classdef MotionScore < dj.Relvar & dj.AutoPopulate properties(Constant) table = dj.Table('aod.MotionScore') diff --git a/schemas/+aod/OrientationResponseSet.m b/schemas/+aod/OrientationResponseSet.m index 91ade3c..46d3cec 100644 --- a/schemas/+aod/OrientationResponseSet.m +++ b/schemas/+aod/OrientationResponseSet.m @@ -7,7 +7,7 @@ responses : longblob # The trial responses (rows) for the cells (columns) %} -classdef OrientationResponseSet < dj.Relvar & dj.Automatic +classdef OrientationResponseSet < dj.Relvar & dj.AutoPopulate properties(Constant) table = dj.Table('aod.OrientationResponseSet') @@ -28,9 +28,9 @@ function makeTuples(self, key) disp 'Sorting data'; [trial_data tuple.cell_nums] = aod.OrientationResponseSet.tracesByOri(key, key.lag, key.bin_duration); - tuple.responses = zeros(1:length(trial_data),count(aod.UniqueCell & key)); + tuple.responses = zeros(length(trial_data),count(aod.UniqueCell & key)); for i = 1:length(trial_data) - tuple.responses(i,:) = mean(trial_data(i).traces,1) + tuple.responses(i,:) = mean(trial_data(i).traces,1); end tuple.orientation = cat(1,trial_data.ori); self.insert(tuple) @@ -60,12 +60,18 @@ function makeTuples(self, key) trial_data = repmat(struct,length(trials) * oriPerTrial,1); k = 1; + if fetch1(stimulation.MultiDimInfo & key, 'block_design') == 0 + first_sub = 2; + else + first_sub = 1; + end + for i = 1:length(trials) trial_info = fetch1(stimulation.StimTrials(trials(i)), 'trial_params'); event = fetch(stimulation.StimTrialEvents(trials(i), 'event_type="showSubStimulus"'),'*'); onset = sort([event.event_time]); - for j = 1:length(event) + for j = first_sub:length(event) cond = trial_info.conditions(j); ori = conditions(cond).condition_info.orientation; condIdx = find(ori == oris); diff --git a/schemas/+aod/PatchedLikelihood.m b/schemas/+aod/PatchedLikelihood.m index fbbff8a..b7e8475 100644 --- a/schemas/+aod/PatchedLikelihood.m +++ b/schemas/+aod/PatchedLikelihood.m @@ -21,7 +21,7 @@ tau_std : double # Expected std of spike time %} -classdef PatchedLikelihood < dj.Relvar & dj.Automatic +classdef PatchedLikelihood < dj.Relvar & dj.AutoPopulate properties(Constant) table = dj.Table('aod.PatchedLikelihood') diff --git a/schemas/+aod/QualitySet.m b/schemas/+aod/QualitySet.m index 41f8713..8a1928f 100644 --- a/schemas/+aod/QualitySet.m +++ b/schemas/+aod/QualitySet.m @@ -25,8 +25,8 @@ function makeTuples( this, key ) pre_vol = pro(acq.AodScan & key, (acq.AodScan * pro(acq.AodVolume)) & 'aod_volume_start_time < aod_scan_start_time', 'MAX(aod_volume_start_time)->aod_volume_start_time'); post_vol = pro(acq.AodScan & key, (acq.AodScan * pro(acq.AodVolume)) & '(aod_volume_start_time > aod_scan_stop_time)', 'MIN(aod_volume_start_time)->aod_volume_start_time'); - pre = fetch(acq.AodVolume & pre_vol, '*'); - post = fetch(acq.AodVolume(post_vol), '*'); + pre = fetch(acq.AodVolume & pre_vol, 'x_coordinate', 'y_coordinate', 'z_coordinate'); + post = fetch(acq.AodVolume(post_vol), 'x_coordinate', 'y_coordinate', 'z_coordinate'); as = fetch(acq.AodScan(key), '*'); diff --git a/schemas/+aod/ReconstructionAccuracy.m b/schemas/+aod/ReconstructionAccuracy.m index d94f7fd..06f54b2 100644 --- a/schemas/+aod/ReconstructionAccuracy.m +++ b/schemas/+aod/ReconstructionAccuracy.m @@ -12,7 +12,7 @@ iqr : double # interquartile range (in delta F / F) %} -classdef ReconstructionAccuracy < dj.Relvar & dj.Automatic +classdef ReconstructionAccuracy < dj.Relvar & dj.AutoPopulate properties(Constant) table = dj.Table('aod.ReconstructionAccuracy') @@ -43,6 +43,7 @@ function makeTuples(self, key) end ds = round(trace.fs / 100); + ds = 1; dt = 1 / trace.fs; ds_trace = decimate(trace.trace,ds); highPass = 0.1; @@ -54,10 +55,11 @@ function makeTuples(self, key) gam = 1 - dt / 2; lam = 10; + sigma = sqrt(var(diff(filtered_trace)) / 2 / ds); fr = fast_oopsi(filtered_trace, ... struct('dt',dt,'est_lam',0,'est_gam',0,'fast_iter_max',0),... - struct('gam',gam,'lam',lam)); + struct('gam',gam,'lam',lam,'sig',sigma)); tuple.trace = fr; tuple.trace_t = trace.t0 + (0:length(fr)-1) / fs * 1000; diff --git a/schemas/+aod/ScanMotion.m b/schemas/+aod/ScanMotion.m index 3a43288..309f3b6 100644 --- a/schemas/+aod/ScanMotion.m +++ b/schemas/+aod/ScanMotion.m @@ -8,7 +8,7 @@ z : longblob # z coordinate of the motion trace %} -classdef ScanMotion < dj.Relvar & dj.Automatic +classdef ScanMotion < dj.Relvar & dj.AutoPopulate properties(Constant) table = dj.Table('aod.ScanMotion') @@ -27,16 +27,24 @@ function makeTuples(self, key) tuple = key; - br = getFile(acq.AodScan & key,'Motion'); - motionData = br(:,:,:); - [T D P] = size(motionData); - mot = reshape(motionData,[T 10 20 P]); - mot = permute(mot, [4 2 3 1]); - t = br(:,'t'); - coordinates = br.motionCoordinates; - + fn = char(fetchn(acq.AodScan & key, 'aod_scan_filename')) ; % get the filename from db + oldversion = acq.checkFileversion(fn) ; % check if motion data is in old format, i.e. 2 planes or in new format, i.e. a volume + if (oldversion) + br = getFile(acq.AodScan & key,'Motion'); + motionData = br(:,:,:); + [T D P] = size(motionData); + mot = reshape(motionData,[T 10 20 P]); + mot = permute(mot, [4 2 3 1]); + t = br(:,'t'); + coordinates = br.motionCoordinates; + - [xpos ypos zpos details] = aod.trackMotion(mot, t, coordinates); + [xpos,ypos,zpos,details] = aod.trackMotion(mot, t, coordinates); + else + % returned positions in microns and time in secs + % the motion data channel is hardcoded to 2, should be eventually guided by some info in the database + [xpos,ypos,zpos,t] = aod.trackMotion3D(fn, 2); + end tuple.x = xpos; tuple.y = ypos; tuple.z = zpos; diff --git a/schemas/+aod/Spikes.m b/schemas/+aod/Spikes.m index f205bba..f2e2f3a 100644 --- a/schemas/+aod/Spikes.m +++ b/schemas/+aod/Spikes.m @@ -41,6 +41,24 @@ function makeTuples(self, key) end end end + + methods + function plot(self) + assert(count(self) == 1, 'Only for one matching relvar'); + + key = fetch(aod.Spikes & self); + + sp = fetch(self, '*'); + trace = fetch(aod.TracePreprocess(key) & aod.TracePreprocessMethod('preprocess_method_name="pc20"'), 'trace'); + times = getTimes(aod.TracePreprocess(key) & aod.TracePreprocessMethod('preprocess_method_name="pc20"')); + + mid = mean(times); + plot((times - mid)/1000,(trace.trace - min(trace.trace)),(sp.times - mid)/1000,zeros(size(sp.times)), 'ok'); + + xlim([0 60]); + end + + end methods (Static) @@ -525,4 +543,4 @@ function makeTuples(self, key) end end -end \ No newline at end of file +end diff --git a/schemas/+aod/TracePreprocessSet.m b/schemas/+aod/TracePreprocessSet.m index 9d72382..958d805 100644 --- a/schemas/+aod/TracePreprocessSet.m +++ b/schemas/+aod/TracePreprocessSet.m @@ -5,7 +5,7 @@ --- %} -classdef TracePreprocessSet < dj.Relvar & dj.Automatic +classdef TracePreprocessSet < dj.Relvar & dj.AutoPopulate properties(Constant) table = dj.Table('aod.TracePreprocessSet'); popRel = aod.TracePreprocessSetParam - acq.AodScanIgnore; diff --git a/schemas/+aod/Traces.m b/schemas/+aod/Traces.m index c5e6920..3577fb5 100644 --- a/schemas/+aod/Traces.m +++ b/schemas/+aod/Traces.m @@ -14,7 +14,7 @@ classdef Traces < dj.Relvar properties(Constant) - table = dj.Table('aod.Traces'); +% table = dj.Table('aod.Traces'); offset = -800 * 600; end diff --git a/schemas/+aod/UniqueCells.m b/schemas/+aod/UniqueCells.m index fe438cb..f7c75b5 100644 --- a/schemas/+aod/UniqueCells.m +++ b/schemas/+aod/UniqueCells.m @@ -5,7 +5,7 @@ duplicates : longblob # list of duplicate cells %} -classdef UniqueCells < dj.Relvar & dj.Automatic +classdef UniqueCells < dj.Relvar & dj.AutoPopulate properties(Constant) table = dj.Table('aod.UniqueCells') popRel = aod.QualitySet; diff --git a/schemas/+aod/checkFileversion.m b/schemas/+aod/checkFileversion.m new file mode 100644 index 0000000..f259dc5 --- /dev/null +++ b/schemas/+aod/checkFileversion.m @@ -0,0 +1,18 @@ + +function oldversion=checkFileVersion(fn) + +oldversion = false ; +try + version = h5readatt(fn,'/','Version') ; + if (version==2) + oldversion = false ; + % motion data are in the new 3d format + else + br = aodReader(fn, 'Motion') ; + sz = size(br.motionCoordinates) ; + if rem(nthroot(sz(1),3),1)==0 % a way to handle files older than version 2 where a volume with cubic dimension was collected + oldversion = false ; + end + end +catch +end diff --git a/schemas/+aod/segmenting.m b/schemas/+aod/segmenting.m index b563ca0..ae9e59c 100644 --- a/schemas/+aod/segmenting.m +++ b/schemas/+aod/segmenting.m @@ -178,4 +178,4 @@ as = fetch(acq.AodScan & sess) -av = fetch(acq.AodVolume & sess) \ No newline at end of file +av = fetch(acq.AodVolume & sess) diff --git a/schemas/+aod/trackMotion.m b/schemas/+aod/trackMotion.m index 098e160..ec28b61 100644 --- a/schemas/+aod/trackMotion.m +++ b/schemas/+aod/trackMotion.m @@ -3,9 +3,18 @@ % % [xpos ypos zpos t details] = trackMotion(fn) % +% mot is a 4-D matrix. The first index is the channel +% for either red or green. The second and third index +% are spatial with the multiple planes being concatenated +% in the third dimension. The fourth dimension is time. +% +% t is the index of time stamps for the frames +% +% coordinates are the locations of the pixels +% +% % JC 2010-07-12 - %% load and normalize stack mot = double(mot); diff --git a/schemas/+aod/trackMotion3D.m b/schemas/+aod/trackMotion3D.m new file mode 100644 index 0000000..d523caa --- /dev/null +++ b/schemas/+aod/trackMotion3D.m @@ -0,0 +1,180 @@ + +function [x,y,z,t] = trackMotion3D(fn, analysischan) + +%% test function computeshifts +% requires z:\libraries\matlab in path +%run(getLocalPath('/lab/libraries/hdf5matlab/setPath.m')) ; + +downsampleFactor = 1 ; % skip frames if > 1 +numFramesToAverageForRef = 10 ; % frames to average to determine the reference frame against which motion is computed +useMultiCore = true ; +postRefStartFrameTime = 0 ; % N secs to skip before analyzing motion signals, in many cases, right after starting the point scan, there is a slight drift + +fn2 = findFile(RawPathMap, fn) +br = aodReader(fn2, 'Motion'); +br1 = AodMotionCorrectionReader(fn2) ; % if motion was corrected by using aod offsets, this is where they would be read + +dat = br(:,:,:) ; +numSamples = br1.sz ; +duration = numSamples(:,1)/br1.Fs ; % sec + +motion_fs = size(dat,1)/duration ; % this is based on the temporal sampling rate +delta_t = downsampleFactor/motion_fs ; + +numPlanes = br.planes ; +coordinates = br.motionCoordinates ; + +% get pmt chan data +dat2 = squeeze(dat(:,:,analysischan)) ; +clearvars dat ; +gridSize = sqrt(size(dat2,2)/numPlanes) ; + +% frame number to start from +postRefStartFrame = max(round(postRefStartFrameTime*motion_fs),1) ; + +% chop off the first N secs +dat2 = dat2(postRefStartFrame:end, :) ; + +needReference = true ; % first get a reference frame +refparams = [] ; +calibparams = [0 1; 0 1; 0 1] ; % dont apply any calibration +shifts = zeros(size(dat2,1),3); % store motion vectors here +dat2(1,:) = mean(dat2(1:numFramesToAverageForRef,:)) ; % average frames for computing reference +[shifts(1,:), refparams] = computeshifts(dat2(1,:)',coordinates,gridSize,calibparams,refparams,needReference) ; % this position of the fitted sphere is treated as the reference + +% setup for parallel processing +if (useMultiCore) + try + mypool = parpool ; % use max processors + catch + fprintf('Unable to start parallel pool, single processor will be used') ; + useMultiCore = false ; + end +end + +% take into account frames to skip, i.e. sampling rate of motion, +% downsample with averaging +adat = zeros(size(dat2)) ; +if (downsampleFactor>1) + count = 0 ; + for ii=numFramesToAverageForRef+downsampleFactor:downsampleFactor:size(dat2,1) + count = count + 1 ; + adat(count,:) = adat(count,:) + sum(dat2(ii-downsampleFactor+1:ii,:)) ; + adat(count,:) = adat(count,:) / downsampleFactor ; + end +else + adat = dat2(numFramesToAverageForRef+1:end,:) ; + count = size(adat,1); +end + +needReference=false ; % no need for a reference, start computing positions relative to the reference of spheres in subsequent motion frames +if (useMultiCore) + parfor ii=1:count + shifts(ii+numFramesToAverageForRef,:) = computeshifts(adat(ii,:)',coordinates,gridSize,calibparams,refparams,needReference) ; + end + delete(mypool) ; +else + for ii=1:count + shifts(ii+numFramesToAverageForRef,:) = computeshifts(adat(ii,:)',coordinates,gridSize,calibparams,refparams,needReference) ; + if (rem(count,100)==0) + fprintf('Processed %d frames', count) ; + end + end +end +x = squeeze(shifts(numFramesToAverageForRef+1:end,2)) ; +y = squeeze(shifts(numFramesToAverageForRef+1:end,1)) ; +z = squeeze(shifts(numFramesToAverageForRef+1:end,3)) ; +t=((0:1:length(x)-1)*delta_t)+(numFramesToAverageForRef/motion_fs) ; % the first motion position is after the last frame used for reference calculations + +end + +% compute position shifts of a 3d object present in the sampled set of +% planes. +% data consists of a single "motion frame" +% data and refdata have to be same size +% variable argument parameters: initial_spaceconstant (def=5) +function [shifts, outrefparams] = computeshifts(data, motcoordinates, gridsize, calibparams, refparams, needsref, varargin) + +%% preprocessing +% [d1,d2] = size(data) ; +% data = reshape(data,gridsize,gridsize,gridsize) ; +% for ii=1:gridsize +% td = data(:,:,ii) ; +% data(:,:,ii) = (td - mean(td(:)))/std(td(:)) ; +% end ; +% data = reshape(data,d1,d2) ; +% data = (data - min(data(:)))/(max(data(:))-min(data(:))) ; +% th = 0.25 ; % arbitrary selection by trial and error +% data = data.*(data>th) ; + +%% +params.initial_spaceconstant = 5 ; +params = parseVarArgs(params, varargin{:}, 'strict'); + +shifts = nan(1,3) ; +outrefparams = nan(1,8) ; + +% functions to fit 3-d gaussian, non-oriented - need to later include +% orientation +gref = @(x,cx,cy,cz) x(1)+(x(2)*exp(-(((cx-x(3)).^2/(2*x(4)^2))... + +((cy-x(5)).^2/(2*x(6)^2))+((cz-x(7)).^2/(2*x(8)^2))))) ; +gdat = @(x,cx,cy,cz,sigmax,sigmay,sigmaz) x(1)+(x(2)*exp(-(((cx-x(3)).^2/(2*sigmax^2))... + +((cy-x(4)).^2/(2*sigmay^2))+((cz-x(5)).^2/(2*sigmaz^2))))) ; +fref = @(x,cx,cy,cz,y) sum((y - gref(x,cx,cy,cz)).^2) ; +fdat = @(x,cx,cy,cz,sigmax,sigmay,sigmaz,y) sum((y - gdat(x,cx,cy,cz,sigmax,sigmay,sigmaz)).^2) ; +options = optimset('MaxFunEvals', 1000000, 'MaxIter', 100000) ; + +% range of motion coordinates +x_min = min(motcoordinates(:,1)) ; +y_min = min(motcoordinates(:,2)) ; +z_min = min(motcoordinates(:,3)) ; +x_max = max(motcoordinates(:,1)) ; +y_max = max(motcoordinates(:,2)) ; +z_max = max(motcoordinates(:,3)) ; + +% interpolation grid to regularize the motion sampling +xg = linspace(x_min,x_max,gridsize) ; +yg = linspace(y_min,y_max,gridsize) ; +zg = linspace(z_min,z_max,gridsize) ; +[yg, xg, zg] = meshgrid(xg,yg,zg) ; +xg = reshape(xg, [], 1) ; +yg = reshape(yg, [], 1) ; +zg = reshape(zg, [], 1) ; + +cx=xg ; +cy=yg ; +cz=zg; + +% first motion "frame" is the reference frame, get params of reference +% frame +spaceconstant = params.initial_spaceconstant ; % initial space constant 5 for cells, 10 for small pollens +if(needsref) + y = data ; + try + refparams = fminsearch(@(x) fref(x,cx,cy,cz,y),[min(y(:)) max(y(:))-min(y(:)) (max(xg)+min(xg))/2 spaceconstant (max(yg)+min(yg))/2 spaceconstant (max(zg)+min(zg))/2 spaceconstant]',options) ; + catch + return ; + end ; +end ; + +% compute shifts in current motion frame relative to reference +sigmax = refparams(4) ; % assume that the space constants of the gaussian during simulated motion remains same as the space constants for reference +sigmay = refparams(6) ; +sigmaz = refparams(8) ; +y = data ; +try + X = fminsearch(@(x) fdat(x,cx,cy,cz,sigmax,sigmay,sigmaz,y),[min(y(:)) max(y(:))-min(y(:)) refparams(3) refparams(5) refparams(7)]') ; + shifts(1) = (X(3) - refparams(3)) ; + shifts(2) = (X(4) - refparams(5)) ; + shifts(3) = (X(5) - refparams(7)) ; +catch + return ; +end ; + +% apply calibration +for ii=1:3 + shifts(ii) = calibparams(ii,1)+(calibparams(ii,2)*shifts(ii)) ; +end ; +outrefparams = refparams ; + +end \ No newline at end of file diff --git a/schemas/+ephys/SpikesAlignedTrial.m b/schemas/+ephys/SpikesAlignedTrial.m index 1b5017f..f3c836a 100644 --- a/schemas/+ephys/SpikesAlignedTrial.m +++ b/schemas/+ephys/SpikesAlignedTrial.m @@ -28,8 +28,14 @@ function makeTuples( this, key ) tuples = dj.struct.join(key,trials); for tuple = tuples' - alignTime = fetchn(stimulation.StimTrialEvents(tuple),'event_time'); - endTime = fetch1(stimulation.StimTrialEvents(setfield(tuple,'event_type','endStimulus')),'event_time'); + alignTime = double(fetchn(stimulation.StimTrialEvents(tuple),'event_time')); + try + endTime = double(fetch1(stimulation.StimTrialEvents(setfield(tuple,'event_type','endStimulus')),'event_time')); + catch + + ev_times = double(fetchn(stimulation.StimTrialEvents(setfield(tuple,'event_type','showSubStimulus')),'event_time')); + endTime = max(ev_times) + 60; + end tuple.spikes_aligned = spikes(spikes > (alignTime - tuple.pre_stim_time) & ... (spikes < (endTime + tuple.post_stim_time))) - double(alignTime); diff --git a/schemas/+stimulation/MultiDimInfo.m b/schemas/+stimulation/MultiDimInfo.m index f62edd1..2aa1ea9 100644 --- a/schemas/+stimulation/MultiDimInfo.m +++ b/schemas/+stimulation/MultiDimInfo.m @@ -13,7 +13,7 @@ classdef MultiDimInfo < dj.Relvar & dj.AutoPopulate properties(Constant) table = dj.Table('stimulation.MultiDimInfo'); - popRel = pro(stimulation.StimTrialGroup & acq.Stimulation('exp_type="MultDimExperiment" OR exp_type="MouseMultiDim"'), stimulation.StimTrials, 'count(1) -> num_trials') & 'num_trials >= 10'; + popRel = pro(stimulation.StimTrialGroup & acq.Stimulation('exp_type="MultDimExperiment" OR exp_type="MouseMultiDim"'), stimulation.StimTrials, 'count(1) -> num_trials') & 'num_trials >= 6'; end methods diff --git a/schemas/+stimulation/NotificationSent.m b/schemas/+stimulation/NotificationSent.m new file mode 100644 index 0000000..0b9e5f8 --- /dev/null +++ b/schemas/+stimulation/NotificationSent.m @@ -0,0 +1,11 @@ +%{ +# table to record results notification to slack +-> stimulation.PerSectorStats +----- +msg = '' : varchar(1024) # message for the notification +sent_ts = CURRENT_TIMESTAMP : timestamp + +%} + +classdef NotificationSent < dj.Manual +end \ No newline at end of file diff --git a/schemas/+stimulation/PerSectorStats.m b/schemas/+stimulation/PerSectorStats.m new file mode 100644 index 0000000..1cbb728 --- /dev/null +++ b/schemas/+stimulation/PerSectorStats.m @@ -0,0 +1,12 @@ +%{ +# table that holds per sector behavior analysis in saccadeflashexperiment +-> stimulation.StimTrialGroup +dot_contrast : decimal(4,1) # contrast of stim +bipolar_contrast : tinyint # set to 1 for both poklarities +--- +session_sector_trials=null : blob # entire analysis +%} + + +classdef PerSectorStats < dj.Manual +end \ No newline at end of file