Code to select the aod motion computation method based on file versio…

.gitignore

@@ -4,3 +4,4 @@ getHost.m
 .*
 *~
 ~*
+startup*.*

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.

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');

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
 

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;

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

processing/sync/detectLcdPhotodiodeFlips.m

@@ -137,3 +137,30 @@
     peakSamples = yp + ndx - p - 1 + b - 1;
     peakTimes = [peakTimes; reshape(br(peakSamples,'t'), [], 1)];    %#ok<AGROW>
 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<AGROW>
+        elseif x(i)>x(idx(end))
+            idx(end)=i;
+        end
+    end
+end
+

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

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]);

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));
 

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';

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;