get_thresholds[source]

get_thresholds(data)

Function that attempts to get the high and low thresholds. Not working very well

from theonerig.synchro.io import *
from theonerig.core import *
from theonerig.utils import *
import matplotlib.pyplot as plt
photodiode_data = load_adc_raw("./files/basic_synchro/photodiode_data", sampling_rate=30000)

Loading the data... 100%      

Supposidly, get_thresholds should provide low and high threshold for the data, but the low_treshold is a sensitive value that should be checked manually in a record

get_first_high[source]

get_first_high(data, threshold)

get_first_high finds the idx of the first frame higher than the threshold

detect_frames[source]

detect_frames(data, low_threshold, high_threshold, increment, do_reverse=True, precision=0.95)

Frame detection (or ON signal detection). Capable of finding frame times produced in a regular fashion:

- data: raw data
- low_threshold: threshold used to detect begginning of each frame.
- high_threshold: threshold used to assign label to the frames, and used to detect the beggining of the reading frame.
- increment: Number of timepoints separating frames. Eg, for 30KHz recording of 60hz stimulus: 30000/60
- do_reverse: boolean to indicate if the reverse detection should be done after detecting the first frame.
- precision: Value that indicates how precise are the events recorded to accelerate the detection.
DLP is very stable (.95) whereas some camera triggers have more jitter (.6). Too low value (bellow 0.5) can
result in an overdetection of frames.

reverse_detection[source]

reverse_detection(data, frame_timepoints, low_threshold, increment, precision=0.95)

Detect frames in the left direction.

extend_timepoints[source]

extend_timepoints(frame_timepoints, n=10)

Extrapolates points to the left. Not really needed now except for the signals idx that would change otherwise (and some starting index were set manually)

error_check[source]

error_check(frame_tp)

Search error by looking at the time between each frame. DLP is regular and odd time reveal misdetections.

detect_frames do frame detection. Works for camera pulses and photodiode data emitted by a DLP. It does it by:

• Finding the first frame higher than a threshold
• Detecting the frames before if flag do_reverse is set to True
• Detect frames
• Assign a binary value of if it's higher than the high threshold
• Do a quick check on the frames to spot weird results

frame_timepoints, frame_signals = detect_frames(photodiode_data, 200, 1000, increment=500)
plt.figure()
plt.plot(photodiode_data)
plt.scatter(frame_timepoints, frame_signals*800+600, c="r")

<matplotlib.collections.PathCollection at 0x1feb1c707f0>

cluster_frame_signals[source]

cluster_frame_signals(data, frame_timepoints, n_cluster=5)

Cluster the frame_timepoints in n_cluster categories depending on the area under the curve.

• data: raw data used to compute the AUC
• frame_timepoints: timepoints delimitating each frame
• n_cluster: Number of cluster for the frame signals

cluster_by_epochs[source]

cluster_by_epochs(data, frame_timepoints, frame_signals, epochs)

Does the same thing as cluster_frame_signals, but working on epochs around which the number of cluster can differ. Useful when a record contains stimuli with different signals sizes.

cluster_by_list[source]

cluster_by_list(data, frame_timepoints, frame_signals, stim_list)

Assign the stimulus identity values from stim_list to the frames in data. stim_list contains only the sequence of stimuli. Those need to be expanded. Opposed to cluster_frame_signals and cluster_by_epochs no AUC operation is performed. Input:

- data: raw data used to compute the stimulus times
- frame_timepoints: timepoints delimiting each frame
- frame_signals: binary 1-D numpy array indicating if high_threshold was passed in 'detect_frames'
- stim_list: 1-D numpy array containing the sequence of the stimuli presented

Output:

- frame_signals: [same size as frame_timepoints] stim_signals list containing the correct value from
    stim_list at every entry

Frame signals are then refined using cluster_frame_signals of the signals to attribute them a value in a defined range

frame_signals = cluster_frame_signals(photodiode_data, frame_timepoints, n_cluster=5)

plt.figure()
plt.plot(photodiode_data[120000:131800])
plt.scatter(frame_timepoints[frame_timepoints>120000]-120000, frame_signals[frame_timepoints>120000]*200+200, c='r')

<matplotlib.collections.PathCollection at 0x1feb222a198>

With the frame detected, we can create our record master, often named reM

ref_timepoints, ref_signals = extend_sync_timepoints(frame_timepoints, frame_signals, up_bound=len(photodiode_data))
reM = RecordMaster([(ref_timepoints, ref_signals)])
print(len(reM[0]))

265

Though the reM we just created is from a tiny portion of real data. From now one we will use a premade reM generated from the same dataset, in full.

reM = import_record("./files/basic_synchro/reM_basic_synchro.h5")

Importing the record master

parse_time[source]

parse_time(time_str, pattern='%y%m%d_%H%M%S')

Default parser of rhd timestamps. (serve as a template too)

get_position_estimate[source]

get_position_estimate(stim_time, record_time, sampling_rate)

Estimate where in the record should a stimulus start, in sample points

record_time = parse_time("200331_170849") #Starting time of that example record found on the filename of the record
print(record_time)

2020-03-31 17:08:49

match_starting_position[source]

match_starting_position(frame_timepoints, frame_signals, stim_signals, estimate_start, search_size=1000)

Search the best matching index between frame_signals and stim_signals. params:

- frame_timepoints: Indexes of the frames in the record
- frame_signals: Signals of the detected frames
- stim_signals: Expected stimulus signals
- estimate_start: Estimated start index
- search_size: Stimulus is searched in frame_signals[idx_estimate-search_size: idx_estimate+search_size]

return:

- best match for the starting position of the stimulus

match_starting_position seaks in the record the first frame of a stimulus. We can use functions from theonerig.synchro.extracting to find out the stimuli used in that record, and get their values

from theonerig.synchro.extracting import get_QDSpy_logs, unpack_stim_npy
log = get_QDSpy_logs("./files/basic_synchro")[0]

flickering_bars_pr WARNING dt of frame #15864 was 50.315 m
flickering_bars_pr WARNING dt of frame #19477 was 137.235 m

print(log.stimuli[2])

checkerboard             eed21bda540934a428e93897908d049e at 2020-03-31 17:09:26

unpacked_checkerboard = unpack_stim_npy("./files/basic_synchro/stimulus_data", "eed21bda540934a428e93897908d049e")
print(unpacked_checkerboard[0].shape, unpacked_checkerboard[1].shape, unpacked_checkerboard[2])

(54039, 1, 18, 32) (54039,) None

get_position_estimate can approximately tell us where the stimulus should be to reduce the search time

estimate_start = get_position_estimate(log.stimuli[2].start_time, record_time, sampling_rate=30000)
print("Estimate position in sample points", estimate_start)

Estimate position in sample points 1110000

stim_start_frame = match_starting_position(reM["main_tp"][0], reM["signals"][0], stim_signals=unpacked_checkerboard[1], estimate_start=estimate_start)
print(stim_start_frame)

2235

display_match[source]

display_match(match_position, reference=None, recorded=None, corrected=None, len_line=50)

Let's see the match we obtain

display_match(stim_start_frame, reference=unpacked_checkerboard[1], recorded=reM["signals"][0])

REF [0]  0 0 0 0 0 0 4 0 4 4 4 0 4 0 4 4 4 0 4 0 4 0 0 0 0 0 0 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1
REC [0]  0 0 0 0 0 0 4 0 4 4 4 0 4 0 4 4 4 0 4 0 4 0 0 0 0 0 0 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1

REF [27019]  1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0
REC [27019]  1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0

REF [53989]  1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 4 4 4 4 4 0 0 0 0 0 0
REC [53989]  1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 4 4 4 4 4 0 0 0 0 0 0

We have a match!! But be sure to check it everytime, as mismatches occurs. Set then stim_start_frame manually

frame_error_correction[source]

frame_error_correction(signals, unpacked, algo='nw', **kwargs)

Correcting the display stimulus frame values. Shifts are first detected with one of shift_detection_conv or shift_detection_NW and applied to the stimulus template. Then single frame mismatch are detected and corrected.

- signals: true signal values recorded
- unpacked: stimulus tuple (inten,marker,shader)
- algo: algorithm for shift detection among [nw, conv]
- **kwargs: extra parameter for shift detection functions

returns: stim_tuple_corrected, shift_log, (error_frames_idx, replacement_idx)

error_frame_matches[source]

error_frame_matches(signals, marker, range_)

Find the frames mismatching and finds in the record the closest frame with an identical signal value

apply_shifts[source]

apply_shifts(unpacked, op_log)

Applies the shifts found by either shift_detection functions

shift_detection_conv[source]

shift_detection_conv(signals, marker, range_)

Detect shifts with a convolution method. First look at how far the next closest frame are, and average it over the record. When the average cross the -1 or 1 threshold, shift the reference accordingly.

shift_detection_NW[source]

shift_detection_NW(signals, marker, simmat_basis=[1, -1, -3, -3, -1], insdel=-10, rowside=20)

Memory optimized Needleman-Wunsch algorithm. Instead of an NN matrix, it uses a N(side*2+1) matrix. Indexing goes slightly differently but result is the same, with far less memory consumption and exection speed scaling better with size of the sequences to align.

We correct the stimulus values with frame_error_correction and it gives us back the changes it made to keep track of the errors made.

signals = reM["signals"][0][stim_start_frame:stim_start_frame+len(unpacked_checkerboard[0])]
corrected_checkerboard, shift_log, error_frames = frame_error_correction(signals, unpacked_checkerboard, algo="nw")
print(shift_log, len(error_frames))

[(53516, 'ins'), (53537, 'del')] 45

chop_stim_edges[source]

chop_stim_edges(first_frame, last_frame, stim_tuple, shift_log, frame_replacement)

Cut out the stimulus parts not containing actual stimulus, and change the idx values of shift_log and frame_replacement to match the new indexing.

detect_calcium_frames[source]

detect_calcium_frames(scanning_data, epoch_threshold=-8)

Detect the timing of the 2P frames, epoch by epoch over a record.