Initial commit
5 changed files+238−0
.gitignoreadded+19−0View file
@@ -0,0 +1,19 @@
1+__pycache__/
2+*.py[cod]
3+*$py.class
4+.Python
5+build/
6+develop-eggs/
7+dist/
8+downloads/
9+eggs/
10+.eggs/
11+lib/
12+lib64/
13+parts/
14+sdist/
15+var/
16+wheels/
17+*.egg-info/
18+.installed.cfg
19+*.egg
compress_ints_lossless.pyadded+18−0View file
@@ -0,0 +1,18 @@
1+def compress_ints_lossless(x, *, method: str = "zstd") -> bytes:
2+ import zstandard as zstd
3+
4+ if method == "zstd":
5+ cctx = zstd.ZstdCompressor(level=22)
6+ return cctx.compress(x.tobytes())
7+ elif method == "zlib":
8+ import zlib
9+ return zlib.compress(x.tobytes(), level=9)
10+ elif method == "lzma":
11+ import lzma
12+ return lzma.compress(x.tobytes(), preset=9)
13+ elif method == "simple_ans":
14+ from simple_ans import ans_encode
15+ encoding = ans_encode(x)
16+ return encoding.bitstream + encoding.symbol_counts.tobytes() + encoding.symbol_values.tobytes()
17+ else:
18+ raise ValueError(f"Unknown method: {method}")
helpers.pyadded+44−0View file
@@ -0,0 +1,44 @@
1+import numpy as np
2+
3+def bandpass_filter(array, *, sampling_frequency, lowcut, highcut) -> np.ndarray:
4+ from scipy.signal import butter, lfilter
5+
6+ nyquist = 0.5 * sampling_frequency
7+ low = lowcut / nyquist
8+ high = highcut / nyquist
9+ b, a = butter(5, [low, high], btype="band")
10+ return lfilter(b, a, array, axis=0) # type: ignore
11+
12+def lowpass_filter(array, *, sampling_frequency, highcut) -> np.ndarray:
13+ from scipy.signal import butter, lfilter
14+
15+ nyquist = 0.5 * sampling_frequency
16+ high = highcut / nyquist
17+ b, a = butter(5, high, btype="low")
18+ return lfilter(b, a, array, axis=0) # type: ignore
19+
20+
21+def highpass_filter(array, *, sampling_frequency, lowcut) -> np.ndarray:
22+ from scipy.signal import butter, lfilter
23+
24+ nyquist = 0.5 * sampling_frequency
25+ low = lowcut / nyquist
26+ b, a = butter(5, low, btype="high")
27+ return lfilter(b, a, array, axis=0) # type: ignore
28+
29+
30+def estimate_noise_level(array: np.ndarray, *, sampling_frequency: float) -> float:
31+ array_filtered = highpass_filter(
32+ array, sampling_frequency=sampling_frequency, lowcut=300
33+ )
34+ MAD = float(
35+ np.median(np.abs(array_filtered.ravel() - np.median(array_filtered.ravel())))
36+ / 0.6745
37+ )
38+ return MAD
39+
40+
41+def compute_entropy_per_sample(a):
42+ _, counts = np.unique(a, return_counts=True)
43+ p = counts / len(a)
44+ return -np.sum(p * np.log2(p))
load_real.pyadded+38−0View file
@@ -0,0 +1,38 @@
1+import numpy as np
2+
3+
4+def load_real_000876(
5+ *, num_samples: int, num_channels: int, start_channel: int
6+) -> np.ndarray:
7+ import lindi
8+ # https://neurosift.app/?p=/nwb&url=https://api.dandiarchive.org/api/assets/7e1de06d-d478-40e2-9b64-9dd04eafaa4c/download/&dandisetId=000876&dandisetVersion=draft
9+ nwb_url = "https://api.dandiarchive.org/api/assets/7e1de06d-d478-40e2-9b64-9dd04eafaa4c/download/"
10+ h5f = lindi.LindiH5pyFile.from_hdf5_file(nwb_url)
11+ ds = h5f["/acquisition/ElectricalSeriesAP/data"]
12+ assert isinstance(ds, lindi.LindiH5pyDataset)
13+ ret = ds[:num_samples, start_channel:start_channel + num_channels]
14+ return ret
15+
16+def load_real_000409(
17+ *, num_samples: int, num_channels: int, start_channel: int
18+) -> np.ndarray:
19+ import lindi
20+ # https://neurosift.app/?p=/nwb&url=https://api.dandiarchive.org/api/assets/c04f6b30-82bf-40e1-9210-34f0bcd8be24/download/&dandisetId=000409&dandisetVersion=draft
21+ nwb_url = "https://api.dandiarchive.org/api/assets/c04f6b30-82bf-40e1-9210-34f0bcd8be24/download/"
22+ h5f = lindi.LindiH5pyFile.from_hdf5_file(nwb_url)
23+ ds = h5f['/acquisition/ElectricalSeriesAp/data']
24+ assert isinstance(ds, lindi.LindiH5pyDataset)
25+ ret = ds[:num_samples, start_channel:start_channel + num_channels]
26+ return ret
27+
28+def load_real_001290(
29+ *, num_samples: int, num_channels: int, start_channel: int
30+):
31+ import lindi
32+ # https://neurosift.app/?p=/nwb&url=https://api.dandiarchive.org/api/assets/78c99d23-da88-4ecd-9086-c488a126eac5/download/&dandisetId=001290&dandisetVersion=draft
33+ nwb_url = "https://api.dandiarchive.org/api/assets/78c99d23-da88-4ecd-9086-c488a126eac5/download/"
34+ h5f = lindi.LindiH5pyFile.from_hdf5_file(nwb_url)
35+ ds = h5f['/acquisition/ElectricalSeriesAPImec/data']
36+ assert isinstance(ds, lindi.LindiH5pyDataset)
37+ ret = ds[:num_samples, start_channel:start_channel + num_channels]
38+ return ret
test1.pyadded+119−0View file
@@ -0,0 +1,119 @@
1+# %%
2+import numpy as np
3+from helpers import bandpass_filter, estimate_noise_level, compute_entropy_per_sample
4+from load_real import load_real_000876, load_real_000409, load_real_001290
5+from compress_ints_lossless import compress_ints_lossless
6+
7+def linear_fit(x, y):
8+ """Perform linear fit with constant term.
9+ Returns coefficients and prediction function."""
10+ from numpy.linalg import lstsq
11+ X = np.column_stack([x, np.ones(len(x))])
12+ coeffs = lstsq(X, y, rcond=None)[0]
13+
14+ def predict(x_new):
15+ X_new = np.column_stack([x_new, np.ones(len(x_new))])
16+ return np.dot(X_new, coeffs)
17+
18+ return coeffs, predict
19+
20+# %%
21+N = 500_000
22+# X = np.round(np.random.randn(N) * 500)
23+
24+# X = load_real_001290(num_samples=N, num_channels=1, start_channel=0)
25+X = load_real_000409(num_samples=N, num_channels=1, start_channel=101)
26+# X = load_real_000876(num_samples=N, num_channels=1, start_channel=45)
27+# X = np.random.randn(len(X)) * 100
28+
29+X = X.astype(np.int16)
30+X = X.flatten()
31+# %%
32+e1 = compute_entropy_per_sample(X)
33+print(f'(raw) Bits per sample: {e1:.2f}')
34+print(f'Ideal compression ratio: {X.itemsize * 8 / e1:.2f}')
35+
36+# %%
37+e1 = compute_entropy_per_sample(np.diff(X))
38+print(f'(raw diff) Bits per sample: {e1:.2f}')
39+print(f'Ideal compression ratio: {X.itemsize * 8 / e1:.2f}')
40+
41+# %%
42+# Actual compression ratio
43+buf_zstd = compress_ints_lossless(np.diff(X), method='zstd')
44+buf_zlib = compress_ints_lossless(np.diff(X), method='zlib')
45+buf_lzma = compress_ints_lossless(np.diff(X), method='lzma')
46+buf_ans = compress_ints_lossless(np.diff(X), method='simple_ans')
47+print(f'Zstd compression ratio: {len(X) * X.itemsize / len(buf_zstd):.2f}')
48+print(f'Zlib compression ratio: {len(X) * X.itemsize / len(buf_zlib):.2f}')
49+print(f'Lzma compression ratio: {len(X) * X.itemsize / len(buf_lzma):.2f}')
50+print(f'Simple ANS compression ratio: {len(X) * X.itemsize / len(buf_ans):.2f}')
51+
52+# %%
53+M = 20
54+# N - M + 1 x M
55+sequences = np.array([X[i:i+M] for i in range(len(X) - 2 * M + 1)])
56+predictors = sequences[:, :M - 1]
57+target = sequences[:, M - 1]
58+
59+# Can choose either linear or quadratic fit
60+# coeffs, predict = quadratic_fit(predictors, y)
61+coeffs, predict = linear_fit(predictors, target)
62+predictions = predict(predictors)
63+predictions = np.round(predictions)
64+residuals = target - predictions
65+residuals = residuals.astype(np.int16)
66+e3 = compute_entropy_per_sample(residuals)
67+print(f'(raw adjusted) Bits per sample: {e3:.2f}')
68+print(f'Ideal compression ratio: {X.itemsize * 8 / e3:.2f}')
69+
70+# %%
71+v = 5
72+lowcut = 300
73+highcut = 6000
74+X2 = bandpass_filter(X - np.median(X), sampling_frequency=30000, lowcut=lowcut, highcut=highcut)
75+noise_level = estimate_noise_level(X2, sampling_frequency=30000)
76+X2b = X2 / noise_level * v
77+X2 = np.round(X2b).astype(np.int16)
78+e2 = compute_entropy_per_sample(X2)
79+print(f'(filtered) Bits per sample: {e2:.2f}')
80+print(f'Ideal compression ratio: {X.itemsize * 8 / e2:.2f}')
81+
82+# %%
83+e2 = compute_entropy_per_sample(np.diff(X2))
84+print(f'(filtered diff) Bits per sample: {e2:.2f}')
85+print(f'Ideal compression ratio: {X.itemsize * 8 / e2:.2f}')
86+
87+# %%
88+M = 20
89+# N - M + 1 x M
90+sequences = np.array([X2[i:i+M] for i in range(len(X) - 2 * M + 1)])
91+predictors = sequences[:, :M - 1]
92+target = sequences[:, M - 1]
93+
94+coeffs, predict = linear_fit(predictors, target)
95+predictions = predict(predictors)
96+predictions = np.round(predictions)
97+residuals = target - predictions
98+residuals = residuals.astype(np.int16)
99+e3 = compute_entropy_per_sample(residuals)
100+print(f'(filtered adjusted) Bits per sample: {e3:.2f}')
101+print(f'Ideal compression ratio: {X.itemsize * 8 / e3:.2f}')
102+# %%
103+# Get the actual compression ratio
104+buf_zstd = compress_ints_lossless(residuals, method='zstd')
105+buf_zlib = compress_ints_lossless(residuals, method='zlib')
106+buf_lzma = compress_ints_lossless(residuals, method='lzma')
107+buf_ans = compress_ints_lossless(residuals, method='simple_ans')
108+print(f'Zstd compression ratio: {len(residuals) * residuals.itemsize / len(buf_zstd):.2f}')
109+print(f'Zlib compression ratio: {len(residuals) * residuals.itemsize / len(buf_zlib):.2f}')
110+print(f'Lzma compression ratio: {len(residuals) * residuals.itemsize / len(buf_lzma):.2f}')
111+print(f'Simple ANS compression ratio: {len(residuals) * residuals.itemsize / len(buf_ans):.2f}')
112+
113+# %%
114+import matplotlib.pyplot as plt
115+plt.figure(figsize=(10, 6))
116+plt.plot(X2[:600])
117+# %%
118+plt.plot(coeffs)
119+# %%