1# %%
2import numpy as np
3from zia_benchmark._filters import bandpass_filter, highpass_filter
4from zia_benchmark._data_loaders import load_real_000876, load_real_000409, load_real_001290
5from zia_benchmark._compress_ints_lossless import compress_ints_lossless
6from zia_benchmark._analysis import linear_fit, compute_entropy_per_sample, estimate_noise_level
7import matplotlib.pyplot as plt
9# %%
10N = 500_000
12channel_number = 101
13X = load_real_000409(num_samples=N, num_channels=1, start_channel=channel_number).flatten()
15# X = load_real_001290(num_samples=N, num_channels=1, start_channel=0).flatten()
16# X = load_real_000876(num_samples=N, num_channels=1, start_channel=45).flatten()
18X = X.astype(np.int16)
20# %%
21plt.figure(figsize=(12, 4))
22plt.plot(X[:2400])
23# %%
24def print_ideal_compression_ratio(X):
25 ee = compute_entropy_per_sample(X)
26 print(f'Ideal compression ratio: {X.itemsize * 8 / ee:.2f} ({ee:.2f} bits per sample)')
28def print_actual_compression_ratios(X):
29 buf_zstd = compress_ints_lossless(X, method='zstd')
30 buf_zlib = compress_ints_lossless(X, method='zlib')
31 buf_lzma = compress_ints_lossless(X, method='lzma')
32 buf_ans = compress_ints_lossless(X, method='simple_ans')
33 print(f'Zstd compression ratio: {len(X) * X.itemsize / len(buf_zstd):.2f}')
34 print(f'Zlib compression ratio: {len(X) * X.itemsize / len(buf_zlib):.2f}')
35 print(f'Lzma compression ratio: {len(X) * X.itemsize / len(buf_lzma):.2f}')
36 print(f'simple_ans compression ratio: {len(X) * X.itemsize / len(buf_ans):.2f}')
38def get_marcovian_prediction_residual(X, M):
39 sequences = np.array([X[i:i+M] for i in range(len(X) - 2 * M + 1)])
40 predictors = sequences[:, :M - 1]
41 target = sequences[:, M - 1]
43 coeffs, predict = linear_fit(predictors, target)
44 predictions = predict(predictors)
45 predictions = np.round(predictions)
46 residuals = target - predictions
47 residuals = residuals.astype(np.int16)
48 return residuals
50# %%
51print('RAW')
52print_ideal_compression_ratio(X)
54# %%
55print('RAW DELTA ENCODING')
56print_ideal_compression_ratio(np.diff(X))
58# %%
59print('RAW DELTA ENCODING - actual compression ratios')
60print_actual_compression_ratios(np.diff(X))
61print_ideal_compression_ratio(np.diff(X))
63# %%
64X_mr = get_marcovian_prediction_residual(X, 20)
65print('RAW MARCOVIAN')
66print_ideal_compression_ratio(X_mr)
68# %%
69v = 0.25 # step size for quantization
70lowcut = 300
71highcut = 6000
72X_filt = bandpass_filter(X - np.median(X), sampling_frequency=30000, lowcut=lowcut, highcut=highcut)
73noise_level = estimate_noise_level(X_filt, sampling_frequency=30000)
74X_filt_normalized = X_filt / noise_level
75X2b = X_filt_normalized / v
76X2 = np.round(X2b).astype(np.int16)
78# %%
79plt.figure(figsize=(12, 4))
80plt.plot(X[:2400])
82# %%
83print('FILTERED (and quantized)')
84print_ideal_compression_ratio(X2)
86# %%
87print('FILTERED DELTA ENCODING')
88print_ideal_compression_ratio(np.diff(X2))
90# %%
91residuals2 = get_marcovian_prediction_residual(X2, 20)
92print('FILTERED MARCOVIAN')
93print_ideal_compression_ratio(residuals2)
95# %%
96import matplotlib.pyplot as plt
97plt.figure(figsize=(12, 4))
98plt.plot(X[:2400])
99plt.title('RAW')
101plt.figure(figsize=(12, 4))
102plt.plot(X2[:2400])
103plt.title('FILTERED')
105plt.figure(figsize=(12, 4))
106plt.plot(residuals2[:2400])
107plt.title('FILTERED MARCOVIAN')
109# %%
110def sliding_max(x, delta):
111 y = np.zeros_like(x)
112 for i in range(len(x)):
113 y[i] = np.max(x[max(0, i - delta):min(len(x), i + delta + 1)])
114 return y
116def smoothed(x, delta):
117 y = np.zeros_like(x)
118 for i in range(len(x)):
119 y[i] = np.mean(x[max(0, i - delta):min(len(x), i + delta + 1)])
120 return y
122cc = [3, 6]
123Y = sliding_max(np.abs(X_filt_normalized), 50)
124Y = smoothed(Y, 20)
125Y = np.minimum(1, np.maximum(0, (Y - cc[0]) / (cc[1] - cc[0])))
126# Y = highpass_filter(Y, sampling_frequency=30000, lowcut=3)
127Y_scaled = X_filt_normalized * Y
128X3b = Y_scaled / v
129X3 = np.round(X3b).astype(np.int16)
131plt.figure(figsize=(12, 4))
132plt.plot(X2[:2400], color='lightgray')
133# plt.plot(Y[4000:5000])
134plt.plot(X3b[:2400])
135# %%
136print('FILTERED (and quantized) with suppression')
137print_ideal_compression_ratio(X3)
138print('')
139print('FILTERED DELTA ENCODING with suppression')
140print_ideal_compression_ratio(np.diff(X3))
141print('')
142print('FILTERED MARCOVIAN with suppression')
143residuals3 = get_marcovian_prediction_residual(X3, 20)
144print_ideal_compression_ratio(residuals3)
145print('ACTUAL FILTERED MARCOVIAN with suppression')
146print_actual_compression_ratios(residuals3)
147# %%
148def get_run_lengths(x):
149 runs = []
150 i = 0
151 current_nonzero_run_length = 0
152 while i < len(x):
153 if np.all(x[i:i+10] == 0):
154 runs.append(current_nonzero_run_length)
155 current_nonzero_run_length = 0
156 j = i
157 while j < len(x) and x[j] == 0:
158 j += 1
159 runs.append(j - i)
160 i = j
161 else:
162 current_nonzero_run_length += 1
163 i += 1
164 if np.max(runs) < 256:
165 return np.array(runs, dtype=np.uint8)
166 if np.max(runs) < 2 ** 16:
167 return np.array(runs, dtype=np.uint16)
168 return np.array(runs, dtype=np.uint32)
170AA = residuals3[residuals3 != 0]
171run_lengths = get_run_lengths(residuals3)
172print(run_lengths, run_lengths.itemsize)
173ee = compute_entropy_per_sample(AA)
174theoretical_size = (len(AA) * ee / 8) + run_lengths.nbytes
175theoretical_compression_ratio = len(residuals3) * X.itemsize / theoretical_size
176print(f'Theoretical compression ratio: {theoretical_compression_ratio:.2f}')