initial
15 changed files+526−168
.gitignoremodified+22−0View file
@@ -1,6 +1,8 @@
1+# Python
12 __pycache__/
23 *.py[cod]
34 *$py.class
5+*.so
46 .Python
57 build/
68 develop-eggs/
@@ -17,3 +19,23 @@ wheels/
1719 *.egg-info/
1820 .installed.cfg
1921 *.egg
22+
23+# Virtual Environment
24+venv/
25+env/
26+ENV/
27+
28+# IDE
29+.idea/
30+.vscode/
31+*.swp
32+*.swo
33+
34+# Testing
35+.coverage
36+htmlcov/
37+.pytest_cache/
38+.mypy_cache/
39+
40+# Misc
41+.DS_Store
LICENSEadded+201−0View file
@@ -0,0 +1,201 @@
1+ Apache License
2+ Version 2.0, January 2004
3+ http://www.apache.org/licenses/
4+
5+ TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
6+
7+ 1. Definitions.
8+
9+ "License" shall mean the terms and conditions for use, reproduction,
10+ and distribution as defined by Sections 1 through 9 of this document.
11+
12+ "Licensor" shall mean the copyright owner or entity authorized by
13+ the copyright owner that is granting the License.
14+
15+ "Legal Entity" shall mean the union of the acting entity and all
16+ other entities that control, are controlled by, or are under common
17+ control with that entity. For the purposes of this definition,
18+ "control" means (i) the power, direct or indirect, to cause the
19+ direction or management of such entity, whether by contract or
20+ otherwise, or (ii) ownership of fifty percent (50%) or more of the
21+ outstanding shares, or (iii) beneficial ownership of such entity.
22+
23+ "You" (or "Your") shall mean an individual or Legal Entity
24+ exercising permissions granted by this License.
25+
26+ "Source" form shall mean the preferred form for making modifications,
27+ including but not limited to software source code, documentation
28+ source, and configuration files.
29+
30+ "Object" form shall mean any form resulting from mechanical
31+ transformation or translation of a Source form, including but
32+ not limited to compiled object code, generated documentation,
33+ and conversions to other media types.
34+
35+ "Work" shall mean the work of authorship, whether in Source or
36+ Object form, made available under the License, as indicated by a
37+ copyright notice that is included in or attached to the work
38+ (an example is provided in the Appendix below).
39+
40+ "Derivative Works" shall mean any work, whether in Source or Object
41+ form, that is based on (or derived from) the Work and for which the
42+ editorial revisions, annotations, elaborations, or other modifications
43+ represent, as a whole, an original work of authorship. For the purposes
44+ of this License, Derivative Works shall not include works that remain
45+ separable from, or merely link (or bind by name) to the interfaces of,
46+ the Work and Derivative Works thereof.
47+
48+ "Contribution" shall mean any work of authorship, including
49+ the original version of the Work and any modifications or additions
50+ to that Work or Derivative Works thereof, that is intentionally
51+ submitted to Licensor for inclusion in the Work by the copyright owner
52+ or by an individual or Legal Entity authorized to submit on behalf of
53+ the copyright owner. For the purposes of this definition, "submitted"
54+ means any form of electronic, verbal, or written communication sent
55+ to the Licensor or its representatives, including but not limited to
56+ communication on electronic mailing lists, source code control systems,
57+ and issue tracking systems that are managed by, or on behalf of, the
58+ Licensor for the purpose of discussing and improving the Work, but
59+ excluding communication that is conspicuously marked or otherwise
60+ designated in writing by the copyright owner as "Not a Contribution."
61+
62+ "Contributor" shall mean Licensor and any individual or Legal Entity
63+ on behalf of whom a Contribution has been received by Licensor and
64+ subsequently incorporated within the Work.
65+
66+ 2. Grant of Copyright License. Subject to the terms and conditions of
67+ this License, each Contributor hereby grants to You a perpetual,
68+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable
69+ copyright license to reproduce, prepare Derivative Works of,
70+ publicly display, publicly perform, sublicense, and distribute the
71+ Work and such Derivative Works in Source or Object form.
72+
73+ 3. Grant of Patent License. Subject to the terms and conditions of
74+ this License, each Contributor hereby grants to You a perpetual,
75+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable
76+ (except as stated in this section) patent license to make, have made,
77+ use, offer to sell, sell, import, and otherwise transfer the Work,
78+ where such license applies only to those patent claims licensable
79+ by such Contributor that are necessarily infringed by their
80+ Contribution(s) alone or by combination of their Contribution(s)
81+ with the Work to which such Contribution(s) was submitted. If You
82+ institute patent litigation against any entity (including a
83+ cross-claim or counterclaim in a lawsuit) alleging that the Work
84+ or a Contribution incorporated within the Work constitutes direct
85+ or contributory patent infringement, then any patent licenses
86+ granted to You under this License for that Work shall terminate
87+ as of the date such litigation is filed.
88+
89+ 4. Redistribution. You may reproduce and distribute copies of the
90+ Work or Derivative Works thereof in any medium, with or without
91+ modifications, and in Source or Object form, provided that You
92+ meet the following conditions:
93+
94+ (a) You must give any other recipients of the Work or
95+ Derivative Works a copy of this License; and
96+
97+ (b) You must cause any modified files to carry prominent notices
98+ stating that You changed the files; and
99+
100+ (c) You must retain, in the Source form of any Derivative Works
101+ that You distribute, all copyright, patent, trademark, and
102+ attribution notices from the Source form of the Work,
103+ excluding those notices that do not pertain to any part of
104+ the Derivative Works; and
105+
106+ (d) If the Work includes a "NOTICE" text file as part of its
107+ distribution, then any Derivative Works that You distribute must
108+ include a readable copy of the attribution notices contained
109+ within such NOTICE file, excluding those notices that do not
110+ pertain to any part of the Derivative Works, in at least one
111+ of the following places: within a NOTICE text file distributed
112+ as part of the Derivative Works; within the Source form or
113+ documentation, if provided along with the Derivative Works; or,
114+ within a display generated by the Derivative Works, if and
115+ wherever such third-party notices normally appear. The contents
116+ of the NOTICE file are for informational purposes only and
117+ do not modify the License. You may add Your own attribution
118+ notices within Derivative Works that You distribute, alongside
119+ or as an addendum to the NOTICE text from the Work, provided
120+ that such additional attribution notices cannot be construed
121+ as modifying the License.
122+
123+ You may add Your own copyright statement to Your modifications and
124+ may provide additional or different license terms and conditions
125+ for use, reproduction, or distribution of Your modifications, or
126+ for any such Derivative Works as a whole, provided Your use,
127+ reproduction, and distribution of the Work otherwise complies with
128+ the conditions stated in this License.
129+
130+ 5. Submission of Contributions. Unless You explicitly state otherwise,
131+ any Contribution intentionally submitted for inclusion in the Work
132+ by You to the Licensor shall be under the terms and conditions of
133+ this License, without any additional terms or conditions.
134+ Notwithstanding the above, nothing herein shall supersede or modify
135+ the terms of any separate license agreement you may have executed
136+ with Licensor regarding such Contributions.
137+
138+ 6. Trademarks. This License does not grant permission to use the trade
139+ names, trademarks, service marks, or product names of the Licensor,
140+ except as required for reasonable and customary use in describing the
141+ origin of the Work and reproducing the content of the NOTICE file.
142+
143+ 7. Disclaimer of Warranty. Unless required by applicable law or
144+ agreed to in writing, Licensor provides the Work (and each
145+ Contributor provides its Contributions) on an "AS IS" BASIS,
146+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
147+ implied, including, without limitation, any warranties or conditions
148+ of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A
149+ PARTICULAR PURPOSE. You are solely responsible for determining the
150+ appropriateness of using or redistributing the Work and assume any
151+ risks associated with Your exercise of permissions under this License.
152+
153+ 8. Limitation of Liability. In no event and under no legal theory,
154+ whether in tort (including negligence), contract, or otherwise,
155+ unless required by applicable law (such as deliberate and grossly
156+ negligent acts) or agreed to in writing, shall any Contributor be
157+ liable to You for damages, including any direct, indirect, special,
158+ incidental, or consequential damages of any character arising as a
159+ result of this License or out of the use or inability to use the
160+ Work (including but not limited to damages for loss of goodwill,
161+ work stoppage, computer failure or malfunction, or any and all
162+ other commercial damages or losses), even if such Contributor
163+ has been advised of the possibility of such damages.
164+
165+ 9. Accepting Warranty or Additional Liability. While redistributing
166+ the Work or Derivative Works thereof, You may choose to offer,
167+ and charge a fee for, acceptance of support, warranty, indemnity,
168+ or other liability obligations and/or rights consistent with this
169+ License. However, in accepting such obligations, You may act only
170+ on Your own behalf and on Your sole responsibility, not on behalf
171+ of any other Contributor, and only if You agree to indemnify,
172+ defend, and hold each Contributor harmless for any liability
173+ incurred by, or claims asserted against, such Contributor by reason
174+ of your accepting any such warranty or additional liability.
175+
176+ END OF TERMS AND CONDITIONS
177+
178+ APPENDIX: How to apply the Apache License to your work.
179+
180+ To apply the Apache License to your work, attach the following
181+ boilerplate notice, with the fields enclosed by brackets "[]"
182+ replaced with your own identifying information. (Don't include
183+ the brackets!) The text should be enclosed in the appropriate
184+ comment syntax for the file format. We also recommend that a
185+ file or class name and description of purpose be included on the
186+ same "printed page" as the copyright notice for easier
187+ identification within third-party archives.
188+
189+ Copyright [yyyy] [name of copyright owner]
190+
191+ Licensed under the Apache License, Version 2.0 (the "License");
192+ you may not use this file except in compliance with the License.
193+ You may obtain a copy of the License at
194+
195+ http://www.apache.org/licenses/LICENSE-2.0
196+
197+ Unless required by applicable law or agreed to in writing, software
198+ distributed under the License is distributed on an "AS IS" BASIS,
199+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
200+ See the License for the specific language governing permissions and
201+ limitations under the License.
README.mdadded+3−0View file
@@ -0,0 +1,3 @@
1+# zia
2+
3+Compression of integer arrays
helpers.pydeleted+0−44View file
@@ -1,44 +0,0 @@
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.pydeleted+0−38View file
@@ -1,38 +0,0 @@
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
pyproject.tomladded+8−0View file
@@ -0,0 +1,8 @@
1+[build-system]
2+requires = ["setuptools>=45", "wheel", "setuptools_scm>=6.2"]
3+build-backend = "setuptools.build_meta"
4+
5+[tool.pytest.ini_options]
6+testpaths = ["tests"]
7+python_files = ["test_*.py"]
8+addopts = "-ra -q"
setup.pyadded+31−0View file
@@ -0,0 +1,31 @@
1+from setuptools import setup, find_packages
2+
3+setup(
4+ name="zia",
5+ version="0.1.0",
6+ packages=find_packages(where="src"),
7+ package_dir={"": "src"},
8+ install_requires=[
9+ "numpy"
10+ ],
11+ extras_require={
12+ "dev": [
13+ "pytest",
14+ ],
15+ },
16+ python_requires=">=3.8",
17+ author="Jeremy Magland",
18+ description="Compression of integer arrays",
19+ long_description=open("README.md").read(),
20+ long_description_content_type="text/markdown",
21+ keywords="neuroscience, compression, signal processing",
22+ classifiers=[
23+ "Development Status :: 3 - Alpha",
24+ "Intended Audience :: Science/Research",
25+ "Programming Language :: Python :: 3",
26+ "Programming Language :: Python :: 3.8",
27+ "Programming Language :: Python :: 3.9",
28+ "Programming Language :: Python :: 3.10",
29+ "Programming Language :: Python :: 3.11"
30+ ],
31+)
src/zia/__init__.pyadded+0−0View file
No changes to the file's content.
src/zia/_analysis.pyadded+61−0View file
@@ -0,0 +1,61 @@
1+import numpy as np
2+from ._filters import highpass_filter
3+
4+
5+def estimate_noise_level(array: np.ndarray, *, sampling_frequency: float) -> float:
6+ """Estimate the noise level of a signal using the median absolute deviation.
7+
8+ Args:
9+ array: Input signal array
10+ sampling_frequency: Sampling frequency in Hz
11+
12+ Returns:
13+ Estimated noise level
14+ """
15+ array_filtered = highpass_filter(
16+ array, sampling_frequency=sampling_frequency, lowcut=300
17+ )
18+ MAD = float(
19+ np.median(np.abs(array_filtered.ravel() - np.median(array_filtered.ravel())))
20+ / 0.6745
21+ )
22+ return MAD
23+
24+
25+def compute_entropy_per_sample(array: np.ndarray) -> float:
26+ """Compute the entropy per sample of a signal.
27+
28+ Args:
29+ array: Input signal array
30+
31+ Returns:
32+ Entropy per sample in bits
33+ """
34+ _, counts = np.unique(array, return_counts=True)
35+ p = counts / len(array)
36+ return float(-np.sum(p * np.log2(p)))
37+
38+
39+from typing import Callable
40+
41+def linear_fit(x: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, Callable[[np.ndarray], np.ndarray]]:
42+ """Perform linear fit with constant term.
43+
44+ Args:
45+ x: Input array of shape (N, M-1) containing M-1 predictors for N samples
46+ y: Target array of shape (N,) containing values to predict
47+
48+ Returns:
49+ Tuple containing:
50+ - coefficients array of shape (M,)
51+ - prediction function that takes x_new and returns predictions
52+ """
53+ from numpy.linalg import lstsq
54+ X = np.column_stack([x, np.ones(len(x))])
55+ coeffs = lstsq(X, y, rcond=None)[0]
56+
57+ def predict(x_new: np.ndarray) -> np.ndarray:
58+ X_new = np.column_stack([x_new, np.ones(len(x_new))])
59+ return np.dot(X_new, coeffs)
60+
61+ return coeffs, predict
compress_ints_lossless.py →src/zia/_compress_ints_lossless.pyrenamed+14−2View file
@@ -1,7 +1,19 @@
11 def compress_ints_lossless(x, *, method: str = "zstd") -> bytes:
2- import zstandard as zstd
2+ """Compress integer data using various compression methods.
33
4+ Args:
5+ x: NumPy array of integers to compress
6+ method: Compression method to use. One of:
7+ - "zstd": Zstandard compression (default)
8+ - "zlib": zlib compression
9+ - "lzma": LZMA compression
10+ - "simple_ans": Simple Asymmetric Numeral Systems
11+
12+ Returns:
13+ Compressed bytes
14+ """
415 if method == "zstd":
16+ import zstandard as zstd
517 cctx = zstd.ZstdCompressor(level=22)
618 return cctx.compress(x.tobytes())
719 elif method == "zlib":
@@ -15,4 +27,4 @@ def compress_ints_lossless(x, *, method: str = "zstd") -> bytes:
1527 encoding = ans_encode(x)
1628 return encoding.bitstream + encoding.symbol_counts.tobytes() + encoding.symbol_values.tobytes()
1729 else:
18- raise ValueError(f"Unknown method: {method}")
30+ raise ValueError(f"Unknown method: {method}")
src/zia/_data_loaders.pyadded+60−0View file
@@ -0,0 +1,60 @@
1+import numpy as np
2+import lindi
3+from typing import cast
4+
5+
6+def load_real_000876(*, num_samples: int, num_channels: int, start_channel: int) -> np.ndarray:
7+ """Load data from DANDI dataset 000876.
8+
9+ Args:
10+ num_samples: Number of samples to load
11+ num_channels: Number of channels to load
12+ start_channel: Starting channel index
13+
14+ Returns:
15+ Array of shape (num_samples, num_channels) containing the loaded data
16+ """
17+ nwb_url = "https://api.dandiarchive.org/api/assets/7e1de06d-d478-40e2-9b64-9dd04eafaa4c/download/"
18+ h5f = lindi.LindiH5pyFile.from_hdf5_file(nwb_url)
19+ ds = h5f["/acquisition/ElectricalSeriesAP/data"]
20+ assert isinstance(ds, lindi.LindiH5pyDataset)
21+ ret = ds[:num_samples, start_channel:start_channel + num_channels]
22+ return cast(np.ndarray, ret)
23+
24+
25+def load_real_000409(*, num_samples: int, num_channels: int, start_channel: int) -> np.ndarray:
26+ """Load data from DANDI dataset 000409.
27+
28+ Args:
29+ num_samples: Number of samples to load
30+ num_channels: Number of channels to load
31+ start_channel: Starting channel index
32+
33+ Returns:
34+ Array of shape (num_samples, num_channels) containing the loaded data
35+ """
36+ nwb_url = "https://api.dandiarchive.org/api/assets/c04f6b30-82bf-40e1-9210-34f0bcd8be24/download/"
37+ h5f = lindi.LindiH5pyFile.from_hdf5_file(nwb_url)
38+ ds = h5f['/acquisition/ElectricalSeriesAp/data']
39+ assert isinstance(ds, lindi.LindiH5pyDataset)
40+ ret = ds[:num_samples, start_channel:start_channel + num_channels]
41+ return cast(np.ndarray, ret)
42+
43+
44+def load_real_001290(*, num_samples: int, num_channels: int, start_channel: int) -> np.ndarray:
45+ """Load data from DANDI dataset 001290.
46+
47+ Args:
48+ num_samples: Number of samples to load
49+ num_channels: Number of channels to load
50+ start_channel: Starting channel index
51+
52+ Returns:
53+ Array of shape (num_samples, num_channels) containing the loaded data
54+ """
55+ nwb_url = "https://api.dandiarchive.org/api/assets/78c99d23-da88-4ecd-9086-c488a126eac5/download/"
56+ h5f = lindi.LindiH5pyFile.from_hdf5_file(nwb_url)
57+ ds = h5f['/acquisition/ElectricalSeriesAPImec/data']
58+ assert isinstance(ds, lindi.LindiH5pyDataset)
59+ ret = ds[:num_samples, start_channel:start_channel + num_channels]
60+ return cast(np.ndarray, ret)
src/zia/_filters.pyadded+56−0View file
@@ -0,0 +1,56 @@
1+from typing import cast
2+import numpy as np
3+from scipy.signal import butter, lfilter
4+
5+
6+def bandpass_filter(array: np.ndarray, *, sampling_frequency: float, lowcut: float, highcut: float) -> np.ndarray:
7+ """Apply a bandpass filter to the input array.
8+
9+ Args:
10+ array: Input signal array
11+ sampling_frequency: Sampling frequency in Hz
12+ lowcut: Lower cutoff frequency in Hz
13+ highcut: Higher cutoff frequency in Hz
14+
15+ Returns:
16+ Filtered signal array
17+ """
18+ nyquist = 0.5 * sampling_frequency
19+ low = lowcut / nyquist
20+ high = highcut / nyquist
21+ b, a = butter(5, [low, high], btype="band")
22+ return cast(np.ndarray, lfilter(b, a, array, axis=0))
23+
24+
25+def lowpass_filter(array: np.ndarray, *, sampling_frequency: float, highcut: float) -> np.ndarray:
26+ """Apply a lowpass filter to the input array.
27+
28+ Args:
29+ array: Input signal array
30+ sampling_frequency: Sampling frequency in Hz
31+ highcut: Cutoff frequency in Hz
32+
33+ Returns:
34+ Filtered signal array
35+ """
36+ nyquist = 0.5 * sampling_frequency
37+ high = highcut / nyquist
38+ b, a = butter(5, high, btype="low")
39+ return cast(np.ndarray, lfilter(b, a, array, axis=0))
40+
41+
42+def highpass_filter(array: np.ndarray, *, sampling_frequency: float, lowcut: float) -> np.ndarray:
43+ """Apply a highpass filter to the input array.
44+
45+ Args:
46+ array: Input signal array
47+ sampling_frequency: Sampling frequency in Hz
48+ lowcut: Cutoff frequency in Hz
49+
50+ Returns:
51+ Filtered signal array
52+ """
53+ nyquist = 0.5 * sampling_frequency
54+ low = lowcut / nyquist
55+ b, a = butter(5, low, btype="high")
56+ return cast(np.ndarray, lfilter(b, a, array, axis=0))
test1.pymodified+63−84View file
@@ -1,119 +1,98 @@
11 # %%
22 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
3+from zia._filters import bandpass_filter
4+from zia._data_loaders import load_real_000876, load_real_000409, load_real_001290
5+from zia._compress_ints_lossless import compress_ints_lossless
6+from zia._analysis import linear_fit, compute_entropy_per_sample, estimate_noise_level
197
208 # %%
219 N = 500_000
22-# X = np.round(np.random.randn(N) * 500)
2310
2411 # X = load_real_001290(num_samples=N, num_channels=1, start_channel=0)
2512 X = load_real_000409(num_samples=N, num_channels=1, start_channel=101)
2613 # X = load_real_000876(num_samples=N, num_channels=1, start_channel=45)
27-# X = np.random.randn(len(X)) * 100
2814
2915 X = X.astype(np.int16)
3016 X = X.flatten()
3117 # %%
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}')
18+def print_ideal_compression_ratio(X):
19+ ee = compute_entropy_per_sample(X)
20+ print(f'Ideal compression ratio: {X.itemsize * 8 / ee:.2f} ({ee:.2f} bits per sample)')
21+
22+def print_actual_compression_ratios(X):
23+ buf_zstd = compress_ints_lossless(X, method='zstd')
24+ buf_zlib = compress_ints_lossless(X, method='zlib')
25+ buf_lzma = compress_ints_lossless(X, method='lzma')
26+ buf_ans = compress_ints_lossless(X, method='simple_ans')
27+ print(f'Zstd compression ratio: {len(X) * X.itemsize / len(buf_zstd):.2f}')
28+ print(f'Zlib compression ratio: {len(X) * X.itemsize / len(buf_zlib):.2f}')
29+ print(f'Lzma compression ratio: {len(X) * X.itemsize / len(buf_lzma):.2f}')
30+ print(f'simple_ans compression ratio: {len(X) * X.itemsize / len(buf_ans):.2f}')
31+
32+def get_marcovian_prediction_residual(X, M):
33+ sequences = np.array([X[i:i+M] for i in range(len(X) - 2 * M + 1)])
34+ predictors = sequences[:, :M - 1]
35+ target = sequences[:, M - 1]
36+
37+ coeffs, predict = linear_fit(predictors, target)
38+ predictions = predict(predictors)
39+ predictions = np.round(predictions)
40+ residuals = target - predictions
41+ residuals = residuals.astype(np.int16)
42+ return residuals
3543
3644 # %%
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}')
45+print('RAW')
46+print_ideal_compression_ratio(X)
4047
4148 # %%
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}')
49+print('RAW DELTA ENCODING')
50+print_ideal_compression_ratio(np.diff(X))
5151
5252 # %%
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}')
53+print('RAW DELTA ENCODING - actual compression ratios')
54+print_actual_compression_ratios(np.diff(X))
55+print_ideal_compression_ratio(np.diff(X))
6956
7057 # %%
71-v = 5
58+X_mr = get_marcovian_prediction_residual(X, 20)
59+print('RAW MARCOVIAN')
60+print_ideal_compression_ratio(X_mr)
61+
62+# %%
63+v = 0.25 # step size for quantization
7264 lowcut = 300
7365 highcut = 6000
7466 X2 = bandpass_filter(X - np.median(X), sampling_frequency=30000, lowcut=lowcut, highcut=highcut)
7567 noise_level = estimate_noise_level(X2, sampling_frequency=30000)
76-X2b = X2 / noise_level * v
68+X2b = X2 / noise_level / v
7769 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}')
8170
8271 # %%
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}')
72+print('FILTERED (and quantized)')
73+print_ideal_compression_ratio(X2)
8674
8775 # %%
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}')
76+print('FILTERED DELTA ENCODING')
77+print_ideal_compression_ratio(np.diff(X2))
78+
10279 # %%
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}')
80+residuals2 = get_marcovian_prediction_residual(X2, 20)
81+print('FILTERED MARCOVIAN')
82+print_ideal_compression_ratio(residuals2)
11283
11384 # %%
11485 import matplotlib.pyplot as plt
115-plt.figure(figsize=(10, 6))
116-plt.plot(X2[:600])
117-# %%
118-plt.plot(coeffs)
86+plt.figure(figsize=(12, 4))
87+plt.plot(X[:800])
88+plt.title('RAW')
89+
90+plt.figure(figsize=(12, 4))
91+plt.plot(X2[:800])
92+plt.title('FILTERED')
93+
94+plt.figure(figsize=(12, 4))
95+plt.plot(residuals2[:800])
96+plt.title('FILTERED MARCOVIAN')
97+
11998 # %%
tests/conftest.pyadded+7−0View file
@@ -0,0 +1,7 @@
1+"""Configuration for pytest."""
2+
3+import os
4+import sys
5+
6+# Add the src directory to the Python path
7+sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "../src")))
tests/test_basic.pyadded+0−0View file
No changes to the file's content.