Skip to content

Commit

Permalink
feat: second order coherence with delays
Browse files Browse the repository at this point in the history
  • Loading branch information
@Peter-Barrow
Peter-Barrow committed Jul 11, 2024
1 parent ee2a2f7 commit 2ea20f3
Show file tree
Hide file tree
Showing 4 changed files with 118 additions and 24 deletions.
2 changes: 1 addition & 1 deletion pyproject.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"

[project]
name = "tangy"
version = "0.6.0"
version = "0.6.1"
description = "Timetag analysing library"
authors = [
{name = "Peter Thomas Barrow", email = "[email protected]"}
Expand Down
21 changes: 21 additions & 0 deletions tangy_src/_tangy.pxd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -331,3 +331,24 @@ cdef extern from "./src/tangy.h":
const f64 radius,
const f64 read_time,
u64* intensities)

void tangy_second_order_coherence(tangy_buffer* t_buf,
const u64 start,
const u64 stop,
const f64 correlation_window,
const f64 resolution,
const u8 signal,
const u8 idler,
const u64 length,
u64* intensities)


void tangy_second_order_coherence_delays(tangy_buffer* t_buf,
const f64 read_time,
const f64 correlation_window,
const f64 resolution,
const u8 signal,
const u8 idler,
const f64* delays,
const u64 length,
u64* intensities)
40 changes: 40 additions & 0 deletions tangy_src/_tangy.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1058,6 +1058,9 @@ def relative_delay(self, channel_a: int, channel_b: int, read_time: float,
(delay_result): Dataclass containing histogram data and fitting results
"""

# PERF: replace time trace with singles()
# NOTE: decay (tau) is ≈ 1 / singles(channelX), so we have:
# tau_a = 1 / singles(channel_a) and tau_b = 1 / singles(channel_b)
resolution_trace: f64n = 1.0
trace: u64n[:] = self.timetrace(
[channel_a, channel_b], read_time, resolution_trace)
Expand All @@ -1066,6 +1069,9 @@ def relative_delay(self, channel_a: int, channel_b: int, read_time: float,

correlation_window: f64n = window
if window is None:
# NOTE: This should then be...
# correlation_window =
# (1 / singles(channel_a)) + (1 / singles(channel_b))
correlation_window = 2 / intensity_average ** 2

length: u64n = round(correlation_window / resolution) - 1
Expand Down Expand Up @@ -1208,6 +1214,40 @@ def joint_delay_histogram(self, signal: int, idler: int, channels: List[int],
(bin_width, bin_width), histogram,
mi, ms, axis, axis)

@cython.ccall
def second_order_coherence(self, signal: int, idler: int, read_time: float,
radius: float, resolution: float,
delays: Optional[List[float]] = None):

length: u64n = u64n(radius / resolution)
correlation_window: f64n = radius / self.resolution

resolution_hist: u64n = u64n(correlation_window / f64n(length))
correlation_window = length * resolution_hist
length *= 2

intensities: u64n[:] = zeros(length, dtype=u64n)
intensities_view: u64[:] = intensities

times = (arange(length) - (length // 2)) * resolution

if delays is None:
start: u64 = self.lower_bound(self.current_time() - read_time)
stop: u64 = self.count

_tangy.tangy_second_order_coherence(
self._ptr_buf, start, stop, correlation_window,
resolution_hist, signal, idler,
length, cython.address(intensities_view[0]))
return (times, intensities)

_delays_view: cython.double[::1] = array(delays, dtype=f64n)
_tangy.tangy_second_order_coherence_delays(
self._ptr_buf, read_time, correlation_window, resolution_hist,
signal, idler, cython.address(_delays_view[0]),
length, cython.address(intensities_view[0]))
return (times, intensities)


_ptu_header_tag_types = {
"tyEmpty8": 0xFFFF0008,
Expand Down
79 changes: 56 additions & 23 deletions tangy_src/src/analysis_base.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -813,6 +813,7 @@ JOIN(stub, joint_delay_histogram)(shared_ring_buffer* const buf,
return count;
}

// TODO: test this
static inline void
JOIN(stub, second_order_coherence)(shared_ring_buffer* const buf,
const slice* data,
Expand All @@ -822,7 +823,6 @@ JOIN(stub, second_order_coherence)(shared_ring_buffer* const buf,
const f64 resolution,
const u8 signal,
const u8 idler,
const f64* delays,
const u64 length,
u64* intensities) {

Expand Down Expand Up @@ -893,7 +893,7 @@ JOIN(stub, second_order_coherence)(shared_ring_buffer* const buf,
for (u64 i = 0; i < N; i++) {
delta =
time_of_arrival -
ringbuffer_u64_get(buffer_idler, i + buffer_idler->head);
ringbuffer_u64_get(buffer_idler, buffer_idler->head - i - 1);
if (delta < correlation_window) {
hist_index = (central_bin - delta / resolution) - 1;
intensities[hist_index] += 1;
Expand All @@ -906,9 +906,9 @@ JOIN(stub, second_order_coherence)(shared_ring_buffer* const buf,
ringbuffer_u64_push(buffer_idler, time_of_arrival);

for (u64 s = 0; s < N; s++) {
delta =
time_of_arrival -
ringbuffer_u64_get(buffer_signal, s + buffer_idler->head);
delta = time_of_arrival -
ringbuffer_u64_get(buffer_signal,
buffer_signal->head - s - 1);
if (delta < correlation_window) {
hist_index = central_bin + delta / resolution;
intensities[hist_index] += 1;
Expand All @@ -925,6 +925,7 @@ JOIN(stub, second_order_coherence)(shared_ring_buffer* const buf,
ringbuffer_u64_deinit(buffer_idler);
}

// TODO: test this
static inline void
JOIN(stub, second_order_coherence_delays)(shared_ring_buffer* const buf,
const slice* data,
Expand All @@ -938,46 +939,78 @@ JOIN(stub, second_order_coherence_delays)(shared_ring_buffer* const buf,
u64* intensities) {

u8 channels[2] = { signal, idler };
i64 resolution_signed[2] = { resolution, -1 * resolution };
u64 offset[2] = { 0, 1 };
// i64 resolution_signed[2] = { resolution, -1 * resolution };
// u64 offset[2] = { 0, 1 };
pattern_iterator pattern =
patternIteratorInit(buf, data, 2, channels, delays, read_time);

u64 delay_bins[2] = { 0, 0 };
binsFromTimeDelays(buf, 2, delays, delay_bins);

u64 conversion_factor = srb_get_conversion_factor(buf);
u64 current_times[2] = { 0, 0 };
for (usize i = 0; i < 2; i++) {
current_times[i] =
arrivalTimeAt(data, conversion_factor, pattern.index[i]);
arrivalTimeAt(data, conversion_factor, pattern.index[i]) +
delay_bins[i];
}

const u32 N = 4096;
// ring buffer for signal and idler (in that order)
ringbuffer_u64* buffers[2] = { ringbuffer_u64_init(N), ringbuffer_u64_init(N) };
ringbuffer_u64* buffers[2] = { ringbuffer_u64_init(N),
ringbuffer_u64_init(N) };

u64 central_bin = length / 2;
u64 delta;
u64 hist_index = 0;

bool in_range = true;
while (in_range == true) {

pattern.oldest = argMin(buf, pattern.limit, current_times, 2);
u8 idx_a = pattern.oldest;
u8 idx_b = (pattern.oldest == 0) ? 1 : 0;

u64 time_of_arrival =
arrivalTimeAt(data, conversion_factor, pattern.index[idx_a]) +
delays[idx_a];
ringbuffer_u64_push(buffers[idx_a], time_of_arrival);
if (channels[pattern.oldest] == signal) {
ringbuffer_u64_push(buffers[0], current_times[0]);

u64 start = buffers[idx_b]->head;
u64 stop = start + N;
for (u64 i = 0; i < N; i++) {
delta =
current_times[0] -
ringbuffer_u64_get(buffers[1], buffers[1]->head - i - 1);

for (u64 i = start; i < stop; i += 1) {
u64 delta = time_of_arrival - ringbuffer_u64_get(buffers[idx_b], i);
if (delta < correlation_window) {
u64 hist_index =
(central_bin + delta / resolution_signed[idx_b]) - offset[idx_b];
intensities[hist_index] += 1;
if (delta < correlation_window) {
hist_index = (central_bin - delta / resolution) - 1;
intensities[hist_index] += 1;
} else {
break;
}
}

} else if (channels[pattern.oldest] == idler) {
ringbuffer_u64_push(buffers[1], current_times[1]);

for (u64 s = 0; s < N; s++) {
delta =
current_times[1] -
ringbuffer_u64_get(buffers[0], buffers[0]->head - s - 1);

if (delta < correlation_window) {
hist_index = central_bin + delta / resolution;
intensities[hist_index] += 1;
} else {
break;
}
}
}

in_range = nextForChannel(data,
&pattern.iters[pattern.oldest],
channels[pattern.oldest],
&pattern.index[pattern.oldest]);

current_times[pattern.oldest] =
arrivalTimeAt(
data, conversion_factor, pattern.index[pattern.oldest]) +
delay_bins[pattern.oldest];
}

patternIteratorDeinit(&pattern);
Expand Down

0 comments on commit 2ea20f3

Please sign in to comment.