-
Notifications
You must be signed in to change notification settings - Fork 30
Expand file tree
/
Copy pathrsync.py
More file actions
186 lines (167 loc) · 10.4 KB
/
Copy pathrsync.py
File metadata and controls
186 lines (167 loc) · 10.4 KB
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
# Class for converting timestamps between recording systems using sync pulses with
# random inter-pulse intervals.
# https://pycontrol.readthedocs.io/en/latest/user-guide/synchronisation
# Dependencies: Python 3, Numpy, Matplotlib, Scikit-learn.
# (c) Thomas Akam 2018-2025. Released under the GPL-3 open source licence.
import numpy as np
import pylab as plt
from sklearn.mixture import GaussianMixture
class RsyncError(Exception):
pass
class Rsync_aligner:
def __init__(
self,
pulse_times_A,
pulse_times_B,
units_A="auto",
units_B="auto",
chunk_size=5,
plot=False,
raise_exception=True,
):
"""Class for converting timestamps between two recording systems
(e.g pyControl and an ephys) using sync pulses with random inter-pulse
intervals recorded on both systems. Typically these sync pulses are generated
by pyControl using the Rsync hardware object and sent to other systems. To use the
Rsync_aligner,instantiate it by providing the sync pulse times recorded by each
system. Timestamps from either system can then be converted into the reference frame
of the other using the A_to_B and B_to_A methods. If the hardware systems use
different units to measure time this can either be specified manually using the units
arguments when the aligner is instantiated, or estimated automatically by setting
the units arguments to 'auto'. When the aligner is instantiated it works out
which pulses in each reference frame correspond to each other by by aligning
short chunks of pulse sequence A with B by minimising the mean squared error
between inter-pulse intervals.
Arguments:
pulse_times_A: The times when sync pulses occured recorded by hardware system A.
pulse_times_B: The times when sync pulses occured recorded by hardware system B.
units_A: The time units used by system A expressed in milliseconds. E.g. if
system A uses units of seconds the *units_A* argument is 1000. If either
of the units_A or units_B arguments is set to 'auto' the units of B
relative to A are estimated automatically.
units_B: The time units used by system B expressed in milliseconds.
plot: Whether to plot information about the alignment.
raise_exception: If *True* an RsyncError exception is raised if no match is found
between the sync pulse sequences.
"""
if units_A == "auto" or units_B == "auto":
# Estimate the units of B relative to A automatically.
raw_intervals_A = np.diff(pulse_times_A)
raw_intervals_B = np.diff(pulse_times_B)
# Exclude very long intervals as likely due to missing pulses.
good_intervals_A = raw_intervals_A[raw_intervals_A < 3 * np.median(raw_intervals_A)]
good_intervals_B = raw_intervals_B[raw_intervals_B < 3 * np.median(raw_intervals_B)]
# Estimate units of B relative to A using the mean of the good intervals.
units_A = 1
units_B = np.mean(good_intervals_A) / np.mean(good_intervals_B)
# Evalute inter-pulse intervals in common units.
intervals_A = np.diff(pulse_times_A) * units_A # Inter-pulse intervals for sequence A
intervals_B = np.diff(pulse_times_B) * units_B # Inter-pulse intervals for sequence B
intervals_B2 = intervals_B**2
# Find alignments of chunks which minimise sum of squared errors.
chunk_starts_A = np.arange(0, len(pulse_times_A) - chunk_size, chunk_size) # Start indices of each chunk of A.
chunk_starts_B = np.zeros(chunk_starts_A.shape, int) # Start indicies of corresponding chunks in B.
chunk_min_mse = np.zeros(chunk_starts_A.shape) # Mean squared error for each chunks best alignment.
chunk_2nd_mse = np.zeros(chunk_starts_A.shape) # Mean sqared error for each chunks 2nd best alignment.
ones_chunk = np.ones(chunk_size)
for i, csA in enumerate(chunk_starts_A):
chunk_A = intervals_A[csA : csA + chunk_size]
mse = (
np.correlate(intervals_B2, ones_chunk, mode="valid")
+ np.sum(chunk_A**2)
- 2 * np.correlate(intervals_B, chunk_A, mode="valid")
) / chunk_size
chunk_starts_B[i] = np.argmin(mse)
sorted_chunk_min_mse = np.sort(mse)
chunk_min_mse[i] = sorted_chunk_min_mse[0]
chunk_2nd_mse[i] = sorted_chunk_min_mse[1]
# Assign chunks to matched and non-matched groups by fitting 2 component Gaussian mixture model
# to log mse distribition of best + second best alignments.
chunk_mse = np.hstack([chunk_min_mse, chunk_2nd_mse])
chunk_mse[chunk_mse == 0] = np.min(chunk_mse[chunk_mse != 0]) # Replace zeros with smallest non zero value.
log_mse = np.log(chunk_mse)
log_mse = log_mse[np.isfinite(log_mse)].reshape(-1, 1)
gmm = GaussianMixture(n_components=2, covariance_type="spherical")
gmm.fit(log_mse)
valid_matches = gmm.predict(log_mse) == np.argmin(gmm.means_) # True for chunks which are valid matches.
# Make arrays of corresponding times.
cor_times_B = np.full(pulse_times_A.shape, np.nan) # B pulse times corresponding to each A pulse.
for csA, csB, valid in zip(chunk_starts_A, chunk_starts_B, valid_matches):
if valid:
cor_times_B[csA : csA + chunk_size] = pulse_times_B[csB : csB + chunk_size]
# Store times of matched sync pulses.
self.matched_times_A = pulse_times_A[~np.isnan(cor_times_B)]
self.matched_times_B = cor_times_B[~np.isnan(cor_times_B)]
# Store empirical units_A/units_B from matched inter-pulse intervals.
self.dAdB = np.sum(np.diff(self.matched_times_A)) / np.sum(np.diff(self.matched_times_B))
# Check quality of alignment.
# Check Difference in GMM means > 3 x sum of standard deviations.
separation_OK = np.abs(gmm.means_[0] - gmm.means_[1])[0] > 3 * np.sum(np.sqrt(gmm.covariances_))
# Check corresponding times are monotonically increacing.
order_OK = np.all(np.diff(self.matched_times_A) > 0) and np.all(np.diff(self.matched_times_B) > 0)
if not (separation_OK and order_OK):
if raise_exception:
raise RsyncError("No match found between inter-pulse interval sequences.")
else:
print("Rsync warning: No match found between inter-pulse interval sequences.")
# Plotting
if plot:
plt.figure(plot if isinstance(plot, int) else 1, figsize=[7, 9]).clf()
plt.subplot2grid((3, 3), (0, 0), rowspan=1, colspan=2)
plt.hist(log_mse[valid_matches], 20, color="b", label="Match")
plt.hist(log_mse[~valid_matches], 20, color="r", label="Non-match")
plt.legend(loc="upper center")
plt.xlabel("Log mean squared error")
plt.ylabel("# chunks")
plt.subplot2grid((3, 3), (0, 2), rowspan=1, colspan=1)
timing_errors = np.diff(cor_times_B) - np.diff(pulse_times_A)
plt.hist(timing_errors[~np.isnan(timing_errors)], 100)
plt.yscale("log", nonpositive="clip")
plt.xlabel("Inter-pulse interval\ndiscrepancy (ms)")
plt.ylabel("# pulses")
plt.subplot2grid((3, 1), (1, 0), rowspan=2, colspan=1)
plt.plot(pulse_times_A, cor_times_B, ".", markersize=2)
plt.xlim(pulse_times_A[0], pulse_times_A[-1])
plt.xlabel("pulse times A")
plt.ylabel("pulse times B")
plt.tight_layout()
def A_to_B(self, times_A, extrapolate=True):
"""Convert times in A reference frame to B reference frame. If extrapolate=True, times
before the first matched sync pulse and after the last matched sync pulse will be
extrapolated, if False they will be nans.
"""
times_B = np.interp(times_A, self.matched_times_A, self.matched_times_B, left=np.nan, right=np.nan)
if extrapolate:
pf = times_A < self.matched_times_A[0] # Mask indicating times pre first matched pulse.
times_B[pf] = (times_A[pf] - self.matched_times_A[0]) / self.dAdB + self.matched_times_B[0]
pl = times_A > self.matched_times_A[-1] # Mask indicating times post last matched pulse.
times_B[pl] = (times_A[pl] - self.matched_times_A[-1]) / self.dAdB + self.matched_times_B[-1]
return times_B
def B_to_A(self, times_B, extrapolate=True):
"""Convert times in B reference frame to A reference frame. If extrapolate=True, times
before the first matched sync pulse and after the last matched sync pulse will be
extrapolated, if False they will be nans.
"""
times_A = np.interp(times_B, self.matched_times_B, self.matched_times_A, left=np.nan, right=np.nan)
if extrapolate:
pf = times_B < self.matched_times_B[0] # Mask indicating times pre first matched pulse.
times_A[pf] = (times_B[pf] - self.matched_times_B[0]) * self.dAdB + self.matched_times_A[0]
pl = times_B > self.matched_times_B[-1] # Mask indicating times post last matched pulse.
times_A[pl] = (times_B[pl] - self.matched_times_B[-1]) * self.dAdB + self.matched_times_A[-1]
return times_A
# --------------------------------------------------------------------------
def simulate_pulses(n_pulse=1000, interval=[100, 1900], units_B=2, noise_SD=2, missing_pulses=False):
"""Simulate a pair of pulse trains timestamps with drift between their timings."""
pulse_times_A = np.cumsum(np.random.randint(*interval, size=n_pulse)).astype(float)
pulse_times_B = units_B * (pulse_times_A + np.cumsum(np.random.normal(scale=noise_SD, size=n_pulse)))
if missing_pulses:
pulse_times_A = np.hstack(
[pulse_times_A[int(n_pulse * 0.05) : int(n_pulse * 0.21)], pulse_times_A[int(n_pulse * 0.33) :]]
)
pulse_times_B = np.hstack(
[
pulse_times_B[: int(n_pulse * 0.74)],
pulse_times_B[int(n_pulse * 0.85) : int(n_pulse * 0.95)],
]
)
return pulse_times_A, pulse_times_B