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