From bddc1681b4688ccd285fda4efd3e777262df4247 Mon Sep 17 00:00:00 2001 From: David Parrini Date: Wed, 21 Oct 2026 12:51:03 -0300 Subject: [PATCH 1/4] Modernize packaging, CI and PyPI release workflow - Add pyproject.toml (hatchling) so the package is built from this repository; the version is read from comtrade.__version__ - Add numpy and pandas extras - Require Python >= 3.11 (the currently supported versions), bump to 0.2.0 - Rename tests/tests.py to tests/test_comtrade.py and run it with pytest - Replace the test workflow with CI: ruff, tests on Python 3.11-3.15 (plus Windows and macOS) and a build/twine check - Add a release workflow that publishes to PyPI via Trusted Publishing - Add Dependabot for GitHub Actions updates - Fix ruff findings (unused variables/import, type comparison) Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EPnfFrvYTuVEM9uDp81XiD --- .github/dependabot.yml | 9 +++++ .github/workflows/ci.yml | 59 ++++++++++++++++++++++++++++ .github/workflows/release.yml | 39 ++++++++++++++++++ .github/workflows/test.yml | 21 ---------- .gitignore | 3 ++ README.md | 26 +++++++++--- comtrade.py | 4 +- pyproject.toml | 53 +++++++++++++++++++++++++ tests/{tests.py => test_comtrade.py} | 5 +-- 9 files changed, 188 insertions(+), 31 deletions(-) create mode 100644 .github/dependabot.yml create mode 100644 .github/workflows/ci.yml create mode 100644 .github/workflows/release.yml delete mode 100644 .github/workflows/test.yml create mode 100644 pyproject.toml rename tests/{tests.py => test_comtrade.py} (99%) diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..478231b --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,9 @@ +version: 2 +updates: + - package-ecosystem: github-actions + directory: / + schedule: + interval: monthly + groups: + actions: + patterns: ["*"] diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..bcd03ca --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,59 @@ +name: CI + +on: + push: + pull_request: + workflow_dispatch: + # Called by release.yml so nothing is published without passing these jobs. + workflow_call: + +permissions: + contents: read + +jobs: + lint: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + with: + persist-credentials: false + - uses: astral-sh/setup-uv@1c37ad07a6a961277cf70c0d37d6f313000f5884 # v10.3.0 + - run: uv run ruff check + + test: + name: test (${{ matrix.python-version }}, ${{ matrix.os }}) + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest] + python-version: ["3.11", "3.12", "3.13", "3.14", "3.15"] + include: + - os: windows-latest + python-version: "3.15" + - os: macos-latest + python-version: "3.15" + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + with: + persist-credentials: false + - uses: astral-sh/setup-uv@1c37ad07a6a961277cf70c0d37d6f313000f5884 # v10.3.0 + with: + python-version: ${{ matrix.python-version }} + - run: uv run pytest + + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + with: + persist-credentials: false + - uses: astral-sh/setup-uv@1c37ad07a6a961277cf70c0d37d6f313000f5884 # v10.3.0 + - run: uv build + - run: uvx twine check --strict dist/* + - name: Smoke-test the wheel + run: uv run --isolated --no-project --with dist/*.whl python -c "import comtrade; print(comtrade.__version__)" + - uses: actions/upload-artifact@cf430e030ddbb5b0abf93d22962f4752f3646cd9 # v7.0.2 + with: + name: dist + path: dist/ diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml new file mode 100644 index 0000000..cdd92e9 --- /dev/null +++ b/.github/workflows/release.yml @@ -0,0 +1,39 @@ +name: Release + +# Publishing a GitHub Release (tag = version in comtrade.py, e.g. "0.1.3") +# runs CI and uploads the built sdist and wheel to PyPI via Trusted Publishing. + +on: + release: + types: [published] + +permissions: + contents: read + +jobs: + ci: + uses: ./.github/workflows/ci.yml + + publish: + needs: ci + runs-on: ubuntu-latest + environment: + name: pypi + url: https://pypi.org/project/comtrade/${{ github.event.release.tag_name }}/ + permissions: + id-token: write + steps: + - uses: actions/download-artifact@9000827ccba6bdab643e8b6fd33ac0654aef8333 # v8.0.2 + with: + name: dist + path: dist/ + - name: Check the release tag matches the package version + env: + TAG: ${{ github.event.release.tag_name }} + run: | + ls dist/ + if [ ! -f "dist/comtrade-${TAG#v}.tar.gz" ]; then + echo "::error::Release tag '$TAG' does not match __version__ in comtrade.py" + exit 1 + fi + - uses: pypa/gh-action-pypi-publish@7245906fb61fa33a13554f9d5fb856fe61758086 # v1.14.3 diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml deleted file mode 100644 index 5875d8a..0000000 --- a/.github/workflows/test.yml +++ /dev/null @@ -1,21 +0,0 @@ -name: Test - -on: [push, pull_request] - -jobs: - test: - runs-on: ubuntu-latest - strategy: - matrix: - python-version: ["3.7", "3.8", "3.9", "3.10", "3.11"] - - steps: - - uses: actions/checkout@v3 - - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v3 - with: - python-version: ${{ matrix.python-version }} - # Any linting using e.g. flake8, black, isort, mypy, prospector, etc. could be done here. - - name: Run tests - run: | - python -m unittest tests/tests.py diff --git a/.gitignore b/.gitignore index e72a96d..155fe5a 100644 --- a/.gitignore +++ b/.gitignore @@ -164,3 +164,6 @@ cython_debug/ #.idea/ # End of https://www.toptal.com/developers/gitignore/api/python + +# uv (library: dev dependencies are resolved fresh, not locked) +uv.lock diff --git a/README.md b/README.md index 727af46..0f9e390 100644 --- a/README.md +++ b/README.md @@ -21,9 +21,9 @@ Or just copy `comtrade.py` from this repository. ### Requirements -- Python 3.6 or newer -- (optional) numpy for speed improvement -- (optional) pandas for data frame support +- Python 3.11 or newer +- (optional) numpy for speed improvement: `pip install "comtrade[numpy]"` +- (optional) pandas for data frame support: `pip install "comtrade[pandas]"` ## How to Use @@ -183,12 +183,28 @@ Feel free to report any bugs you find. You are welcome to fork and submit pull r ## Development -To run tests, use Python's `unittest`. From a clone of the GitHub repository, run the command: +From a clone of the GitHub repository, run the tests and the linter with [uv](https://docs.astral.sh/uv/): ```sh -python3 -m unittest ./tests/tests.py +uv run pytest +uv run ruff check ``` +Or, with pip 25.1 or newer: + +```sh +python -m pip install -e . --group dev +python -m pytest +``` + +### Releasing + +1. Update `__version__` in `comtrade.py` and push the commit to `master`. +2. Publish a GitHub release tagged with that same version, e.g. `gh release create 0.1.3 --generate-notes`. + +The `Release` workflow then runs CI and uploads the sdist and wheel to PyPI using +[Trusted Publishing](https://docs.pypi.org/trusted-publishers/), so no API token is needed. + ## License The module is available at [GitHub](https://github.com/dparrini/python-comtrade) under the MIT license. diff --git a/comtrade.py b/comtrade.py index 3babda3..f228ef2 100644 --- a/comtrade.py +++ b/comtrade.py @@ -21,6 +21,7 @@ # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE. +__version__ = "0.2.0" import array import datetime as dt @@ -938,7 +939,6 @@ def _load_cff(self, cff_file_path: str, **kwargs): ftype = groups[0] if groups[1] is not None: fformat = groups[1] - fbytes_obj = groups[2] if ftype == "DAT" and fformat != _TYPE_ASCII: break @@ -1183,7 +1183,7 @@ def __init__(self, **kwargs): def parse(self, contents): """Parse a ASCII file contents.""" # Check if contents has been read as an io.BytesIO object, if so, decode into a string. - contents = contents.decode() if (type(contents) == bytes) else contents + contents = contents.decode() if (type(contents) is bytes) else contents analog_count = self._cfg.analog_count status_count = self._cfg.status_count diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..6ef36de --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,53 @@ +[build-system] +requires = ["hatchling>=1.27"] +build-backend = "hatchling.build" + +[project] +name = "comtrade" +dynamic = ["version"] +description = "A module designed to read Common Format for Transient Data Exchange (COMTRADE) file format." +readme = "README.md" +license = "MIT" +license-files = ["LICENSE.md"] +authors = [{ name = "David Parrini", email = "d.parrini@gmail.com" }] +requires-python = ">=3.11" +keywords = ["comtrade", "oscillography", "power systems", "IEEE C37.111", "IEC 60255-24"] +classifiers = [ + "Operating System :: OS Independent", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3 :: Only", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", + "Programming Language :: Python :: 3.14", + "Programming Language :: Python :: 3.15", + "Topic :: Scientific/Engineering", +] +dependencies = [] + +[project.optional-dependencies] +numpy = ["numpy"] +pandas = ["pandas"] + +[project.urls] +Homepage = "https://github.com/dparrini/python-comtrade" +Issues = "https://github.com/dparrini/python-comtrade/issues" +Changelog = "https://github.com/dparrini/python-comtrade/releases" + +[dependency-groups] +dev = ["pytest>=8", "ruff"] + +[tool.hatch.version] +path = "comtrade.py" + +[tool.hatch.build.targets.wheel] +only-include = ["comtrade.py"] + +[tool.hatch.build.targets.sdist] +include = ["comtrade.py", "tests/", "README.md", "LICENSE.md"] + +[tool.pytest.ini_options] +testpaths = ["tests"] + +[tool.ruff.lint] +select = ["E4", "E7", "E9", "F"] diff --git a/tests/tests.py b/tests/test_comtrade.py similarity index 99% rename from tests/tests.py rename to tests/test_comtrade.py index 2445c12..d58b7fc 100644 --- a/tests/tests.py +++ b/tests/test_comtrade.py @@ -2,7 +2,6 @@ import math import os import struct -import time import unittest import comtrade from comtrade import Comtrade @@ -355,7 +354,7 @@ def status(t: float) -> bool: # Load file file_path = os.path.abspath("tests/{}".format(self.filename)) self.comtrade = Comtrade(ignore_warnings=True) - self.comtrade.load(file_path + ".cfg".format(self.filename)) + self.comtrade.load(file_path + ".cfg") def tearDown(self): # Remove temporary files. @@ -378,7 +377,7 @@ def test_status_channels(self): def test_max_analog_value(self): tolerance = 2 - self.assertLessEqual(100 - max(self.comtrade.analog[0]), 2) + self.assertLessEqual(100 - max(self.comtrade.analog[0]), tolerance) def test_last_status_value(self): self.assertEqual(self.comtrade.status[0][-1], 1) From 3e1f089c34ee02a6871f07b489206d1194ec03f2 Mon Sep 17 00:00:00 2001 From: David Parrini Date: Wed, 21 Oct 2026 14:00:03 -0300 Subject: [PATCH 2/4] Fix multiple sample rates, skew and DAT reading issues #48 #49 #29 - Compute sample times across multiple sample rate segments: each segment continues from the last sample of the previous one, instead of (n - 1) / samp with the segment's own rate (see #45) - Accept a blank (non-critical) timestamp in ASCII DAT files and compute it from the sample rate (#48) - Fix loading binary DAT files with status channels only (#49) - Add Comtrade.analog_time with each analog channel's time plus its time skew (#29) - Read ASCII DAT files in blocks and unpack binary DAT fields for all rows at once (numpy is used when use_numpy_arrays is set), replacing the struct row-by-row reader - Update README supported features table and add tests Fixes #48, fixes #49, fixes #29 Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01VXtokXPfpX3vJri639eqXo --- README.md | 11 +- comtrade.py | 263 ++++++++++++++++++++++++----------------- tests/test_comtrade.py | 250 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 414 insertions(+), 110 deletions(-) diff --git a/README.md b/README.md index 0f9e390..ddec2e9 100644 --- a/README.md +++ b/README.md @@ -73,6 +73,8 @@ The module function `load` must be called with the `CFF` file path specified. `Comtrade.analog` and `Comtrade.status` lists stores analog and status channel sample lists respectively. These can be accessed through zero-based indexes, i.e., `Comtrade.analog[0]`. The list `Comtrade.time` stores each sample time in seconds. +Analog channels may have a time skew (in microseconds) from the start of the sample period, which is not included in `Comtrade.time`. The list `Comtrade.analog_time` stores the time values of each analog channel, with its skew, in seconds, i.e., `Comtrade.analog_time[0]` for `Comtrade.analog[0]`. + More information can be accessed through `Comtrade.cfg` object, which stores data such as detailed channel information. Data of additional sections, such as HDR and INF, can be accessed through `hdr` and `inf` properties, respectively. @@ -121,17 +123,16 @@ Feel free to pull requests implementing one of these unsupported features or fix | Timestamp multiplication factor | | x | x | x | x | | Time code and local code | | | | x | x | | Time quality of the samples | | | | x | x | -| Analog channel time skew | | x | x | x | Partial | +| Analog channel time skew | | x | x | x | x | | Analog channel primary and secondary VT or CT ratio | | x | x | x | x | | Status channel phase and monitored circuit | | x | x | x | x | -| Multiple sample rates | x | x | x | x | Partial | +| Multiple sample rates | x | x | x | x | x | | Nanoseconds scale | | | x | x | x | ### Unsupported features * Nanoseconds time base within Python's `datetime` objects (such as `start_timestamp` and `trigger_timestamp` properties). It warns the user but doesn't use it, truncating the numbers. -* Use of multiple sample rates in time calculations for binary data. ### Additional settings and features @@ -145,7 +146,9 @@ in `Comtrade` object constructor: obj = Comtrade(use_numpy_arrays=True) ``` -It may improve performance for computations after loading data. +It may improve performance for computations after loading data. It also speeds up loading binary +DAT files (`BINARY`, `BINARY32` and `FLOAT32` formats), which are converted with vectorized numpy +operations. #### Replace missing data with `float('nan')` diff --git a/comtrade.py b/comtrade.py index f228ef2..0f96d1a 100644 --- a/comtrade.py +++ b/comtrade.py @@ -27,10 +27,11 @@ import datetime as dt import errno import io +import itertools import math +import operator import os import re -import struct import sys from typing import Union import warnings @@ -71,6 +72,9 @@ # common separator character of data fields of CFG and ASCII DAT files _SEPARATOR = "," +# maximum number of ASCII DAT file values converted at once +_ASCII_BLOCK_VALUES = 1000000 + # timestamp regular expression _re_date = re.compile(r"([0-9]{1,2})/([0-9]{1,2})/([0-9]{2,4})") _re_time = re.compile(r"([0-9]{1,2}):([0-9]{2}):([0-9]{1,2})(.([0-9]{1,12}))?") @@ -103,6 +107,27 @@ def _preallocate_values(array_type, size, use_numpy_arrays): return array.array(array_type, [0]) * size +def _copy_values(target, start: int, values): + """Copy values into a preallocated values array, from the start index.""" + if isinstance(target, array.array): + values = array.array(target.typecode, values) + target[start:start + len(values)] = values + + +def _unpack_column(data, row_size: int, offset: int, typecode: str) -> array.array: + """Unpack the little-endian field of typecode at the offset (in bytes) of + every row of row_size bytes in data.""" + size = array.array(typecode).itemsize + column = bytearray(len(data) // row_size * size) + # gather each byte of the field from all rows at once + for k in range(size): + column[k::size] = data[offset + k::row_size] + values = array.array(typecode, column) + if sys.byteorder == "big": + values.byteswap() + return values + + def _prevent_null(str_value: str, value_type: type, default_value): if len(str_value.strip()) == 0: return default_value @@ -594,6 +619,9 @@ def __init__(self, **kwargs): ) self._analog_values = [] self._status_values = [] + # time values of each analog channel, with the channel time skew + # (computed when first accessed) + self._analog_time_values = None # Additional CFF data (or additional comtrade files) self._hdr = None @@ -659,6 +687,15 @@ def time(self) -> list: """Return the time values list.""" return self._time_values + @property + def analog_time(self) -> list: + """Return the time values list of each analog channel, i.e., the + sample time plus the channel time skew.""" + if self._analog_time_values is None: + self._analog_time_values = [self._get_skewed_time(channel.skew) + for channel in self._cfg.analog_channels] + return self._analog_time_values + @property def analog(self) -> list: """Return the analog channel values bidimensional list.""" @@ -794,6 +831,15 @@ def _dat_extract_data(self, dat) -> None: self._analog_values = dat.analog self._status_values = dat.status self._total_samples = dat.total_samples + self._analog_time_values = None + + def _get_skewed_time(self, skew: float): + """Return the time values shifted by a channel time skew (in microseconds).""" + skew_seconds = skew * 1E-6 + if _HAS_NUMPY and isinstance(self._time_values, numpy.ndarray): + return self._time_values + skew_seconds + return array.array(self._time_values.typecode, + [t + skew_seconds for t in self._time_values]) def to_dataframe(self, **kwargs) -> "pd.DataFrame": """Return a pandas DataFrame object with comtrade data.""" @@ -1111,6 +1157,17 @@ def _preallocate(self): steps = self._cfg.sample_rates[-1][1] # last samp field self._total_samples = steps + # first sample number and its time for each sample rate segment. + # Each sample is 1/samp after the previous one, where samp is the + # sample rate of the segment the sample belongs to. + self._rate_segments = [] + start_n, start_time = 1, 0.0 + for samp, endsamp in self._cfg.sample_rates: + self._rate_segments.append((samp, endsamp, start_n, start_time)) + if samp != 0.0: + start_time += (endsamp - start_n) / samp + start_n = endsamp + # analog and status count analog_count = self._cfg.analog_count status_count = self._cfg.status_count @@ -1134,24 +1191,24 @@ def _preallocate(self): self.status[i] = _preallocate_values("i", steps, self._use_numpy_arrays) - def _get_samp(self, n) -> float: - """Get the sampling rate for a sample n (1-based index).""" - # TODO: make tests. - last_sample_rate = 1.0 - for samp, endsamp in self._cfg.sample_rates: + def _get_rate_segment(self, n) -> tuple: + """Get the sampling rate, first sample number and its time for the + sample rate segment of a sample n (1-based index).""" + for samp, endsamp, start_n, start_time in self._rate_segments: if n <= endsamp: - return samp - return last_sample_rate + return samp, start_n, start_time + # samples beyond the last endsamp continue at the last sample rate + samp, _, start_n, start_time = self._rate_segments[-1] + return samp, start_n, start_time def _get_time(self, n: int, ts_value: float, time_base: float, time_multiplier: float): # TODO: add option to enforce dat file timestamp, when available. - # TODO: make tests. - sample_rate = self._get_samp(n) + sample_rate, start_n, start_time = self._get_rate_segment(n) if not self._cfg.timestamp_critical or ts_value == _TIMESTAMP_MISSING: # if the timestamp is missing, use calculated. if sample_rate != 0.0: - return (n - 1) / sample_rate + return start_time + (n - start_n) / sample_rate else: raise ComtradeError("Missing timestamp and no sample rate " "provided.") @@ -1159,9 +1216,6 @@ def _get_time(self, n: int, ts_value: float, time_base: float, # Use provided timestamp if it's not missing return ts_value * time_base * time_multiplier - def filter_missing(self, value) -> float: - return float(value) if value != self.DATA_MISSING else float('nan') - def parse(self, contents): """Virtual method, parse DAT file contents.""" pass @@ -1189,6 +1243,8 @@ def parse(self, contents): status_count = self._cfg.status_count time_mult = self._cfg.timemult time_base = self._cfg.time_base + missing = self.DATA_MISSING + nan = float('nan') # auxiliary vectors (channels gains and offsets) a = [x.a for x in self._cfg.analog_channels] @@ -1199,29 +1255,37 @@ def parse(self, contents): lines = contents.splitlines() else: lines = contents - - line_number = 0 - for line in lines: - line_number = line_number + 1 - if line_number > self._total_samples: + lines = iter(lines) + + # Lines are read in blocks and the values of each channel in a block + # are converted at once, which is much faster than one by one. + block_size = max(1, _ASCII_BLOCK_VALUES // (analog_count + status_count + 2)) + irow = 0 + while irow < self._total_samples: + rows = [line.strip().split(self.ASCII_SEPARATOR) for line in + itertools.islice(lines, min(block_size, self._total_samples - irow))] + if not rows: break - values = line.strip().split(self.ASCII_SEPARATOR) - - n = int(values[0]) - # Read time - ts_val = float(values[1]) - ts = self._get_time(n, ts_val, time_base, time_mult) - avalues = [value * a[i] + b[i] if not math.isnan(value) else float('nan') - for i, value in enumerate(map(self.filter_missing, values[2:analog_count + 2]))] - svalues = [int(x) for x in values[len(values) - status_count:]] + # Read time (the non-critical timestamp may be left blank) + _copy_values(self.time, irow, [ + self._get_time(int(values[0]), + float(values[1]) if values[1].strip() else _TIMESTAMP_MISSING, + time_base, time_mult) + for values in rows]) - # store - self.time[line_number - 1] = ts for i in range(analog_count): - self.analog[i][line_number - 1] = avalues[i] + ai, bi = a[i], b[i] + _copy_values(self.analog[i], irow, + [ai * float(value) + bi if value != missing else nan + for value in map(operator.itemgetter(i + 2), rows)]) + + # status values are the last ones of each line for i in range(status_count): - self.status[i][line_number - 1] = svalues[i] + _copy_values(self.status[i], irow, + list(map(int, map(operator.itemgetter(i - status_count), rows)))) + + irow += len(rows) class _BinaryDatReader(_DatReader): @@ -1235,6 +1299,10 @@ def __init__(self, **kwargs): self.TIME_BYTES = 4 self.SAMPLE_NUMBER_BYTES = 4 + # array typecodes of the analog and the 16 bits status values + self.ANALOG_TYPECODE = "h" + self.STATUS_TYPECODE = "H" + if self._rev_year == REV_1991: self.DATA_MISSING = -1 # 0xFFFF else: @@ -1242,30 +1310,24 @@ def __init__(self, **kwargs): self.read_mode = "rb" - if struct.calcsize("L") == 4: - self.STRUCT_FORMAT = "LL {acount:d}h {dcount:d}H" - self.STRUCT_FORMAT_ANALOG_ONLY = "LL {acount:d}h" - self.STRUCT_FORMAT_STATUS_ONLY = "LL {dcount:d}H" - else: - self.STRUCT_FORMAT = "II {acount:d}h {dcount:d}H" - self.STRUCT_FORMAT_ANALOG_ONLY = "II {acount:d}h" - self.STRUCT_FORMAT_STATUS_ONLY = "II {dcount:d}H" - - def get_reader_format(self, analog_channels, status_bytes): - # Number of status fields of 2 bytes based on the total number of - # bytes. - dcount = math.floor(status_bytes / 2) - - # Check the file configuration - if int(status_bytes) > 0 and int(analog_channels) > 0: - return self.STRUCT_FORMAT.format(acount=analog_channels, - dcount=dcount) - elif int(analog_channels) > 0: - # Analog channels only. - return self.STRUCT_FORMAT_ANALOG_ONLY.format(acount=analog_channels) - else: - # Status channels only. - return self.STRUCT_FORMAT_STATUS_ONLY.format(acount=dcount) + def _get_analog_values(self, yint, a: float, b: float): + """Return the analog values a * yint + b of the raw values yint, + with NaN for missing raw values.""" + if _HAS_NUMPY and self._use_numpy_arrays: + y = numpy.asarray(yint, dtype=numpy.float64) + is_missing = y == self.DATA_MISSING + y = a * y + b + y[is_missing] = numpy.nan + return y + missing = self.DATA_MISSING + nan = float('nan') + return [a * value + b if value != missing else nan for value in yint] + + def _get_status_values(self, group, chnindex: int): + """Return the status values of the bit chnindex of a group of 16 bits.""" + if _HAS_NUMPY and self._use_numpy_arrays: + return (numpy.asarray(group) >> chnindex) & 1 + return [(word >> chnindex) & 1 for word in group] def parse(self, contents): """Parse DAT binary file contents.""" @@ -1281,47 +1343,48 @@ def parse(self, contents): dbytes = self.STATUS_BYTES * math.ceil(schannel / 16.0) groups_of_16bits = math.floor(dbytes / self.STATUS_BYTES) - # Struct format. - row_reader = struct.Struct(self.get_reader_format(achannels, dbytes)) - - # Row reading function. if isinstance(contents, io.TextIOBase) or \ isinstance(contents, io.BufferedIOBase): # Read all buffer contents contents = contents.read() - for irow, values in enumerate(row_reader.iter_unpack(contents)): - # Sample number - n = values[0] - # Time stamp - ts_val = values[1] - ts = self._get_time(n, ts_val, time_base, time_mult) - - if irow >= self.total_samples: - break - self.time[irow] = ts - - # Extract analog channel values. - for ichannel in range(achannels): - yint = self.filter_missing(values[ichannel + 2]) - y = a[ichannel] * yint + b[ichannel] if not math.isnan(yint) else float('nan') - self.analog[ichannel][irow] = y - - # Extract status channel values. - for igroup in range(groups_of_16bits): - group = values[achannels + 2 + igroup] - - # for each group of 16 bits, extract the status channels - maxchn = min([(igroup + 1) * 16, schannel]) - for ichannel in range(igroup * 16, maxchn): - chnindex = ichannel - igroup * 16 - mask = int('0b01', 2) << chnindex - extract = (group & mask) >> chnindex - - self.status[ichannel][irow] = extract - - # Get the next row - irow += 1 + # Each row has the sample number, the time stamp, the analog values + # and the groups of 16 status bits. Only whole rows are read, up to + # the total samples. + row_size = self.SAMPLE_NUMBER_BYTES + self.TIME_BYTES + \ + achannels * self.ANALOG_BYTES + dbytes + rows = min(len(contents) // row_size, self.total_samples) + contents = contents[:rows * row_size] + + # Each field is unpacked from all rows at once, which is much faster + # than unpacking the rows one by one. + sample_numbers = _unpack_column(contents, row_size, 0, "I") + timestamps = _unpack_column(contents, row_size, + self.SAMPLE_NUMBER_BYTES, "I") + _copy_values(self.time, 0, [self._get_time(n, ts_val, time_base, time_mult) + for n, ts_val in zip(sample_numbers, timestamps)]) + + # Extract analog channel values. + offset = self.SAMPLE_NUMBER_BYTES + self.TIME_BYTES + for ichannel in range(achannels): + yint = _unpack_column(contents, row_size, offset, + self.ANALOG_TYPECODE) + _copy_values(self.analog[ichannel], 0, + self._get_analog_values(yint, a[ichannel], b[ichannel])) + offset += self.ANALOG_BYTES + + # Extract status channel values. + for igroup in range(groups_of_16bits): + group = _unpack_column(contents, row_size, offset, + self.STATUS_TYPECODE) + + # for each group of 16 bits, extract the status channels + maxchn = min((igroup + 1) * 16, schannel) + for ichannel in range(igroup * 16, maxchn): + chnindex = ichannel - igroup * 16 + _copy_values(self.status[ichannel], 0, + self._get_status_values(group, chnindex)) + offset += self.STATUS_BYTES class _Binary32DatReader(_BinaryDatReader): @@ -1331,13 +1394,7 @@ def __init__(self, **kwargs): # Call the initialization for the inherited class super().__init__(**kwargs) self.ANALOG_BYTES = 4 - - if struct.calcsize("L") == 4: - self.STRUCT_FORMAT = "LL {acount:d}l {dcount:d}H" - self.STRUCT_FORMAT_ANALOG_ONLY = "LL {acount:d}l" - else: - self.STRUCT_FORMAT = "II {acount:d}i {dcount:d}H" - self.STRUCT_FORMAT_ANALOG_ONLY = "II {acount:d}i" + self.ANALOG_TYPECODE = "i" # maximum negative value self.DATA_MISSING = -2147483648 # 0x80000000 @@ -1350,13 +1407,7 @@ def __init__(self, **kwargs): # Call the initialization for the inherited class super().__init__(**kwargs) self.ANALOG_BYTES = 4 - - if struct.calcsize("L") == 4: - self.STRUCT_FORMAT = "LL {acount:d}f {dcount:d}H" - self.STRUCT_FORMAT_ANALOG_ONLY = "LL {acount:d}f" - else: - self.STRUCT_FORMAT = "II {acount:d}f {dcount:d}H" - self.STRUCT_FORMAT_ANALOG_ONLY = "II {acount:d}f" + self.ANALOG_TYPECODE = "f" # Maximum negative value self.DATA_MISSING = sys.float_info.min diff --git a/tests/test_comtrade.py b/tests/test_comtrade.py index d58b7fc..28ac3b1 100644 --- a/tests/test_comtrade.py +++ b/tests/test_comtrade.py @@ -3,6 +3,7 @@ import os import struct import unittest +from unittest import mock import comtrade from comtrade import Comtrade @@ -53,6 +54,70 @@ COMTRADE_SAMPLE_4_CFG_FILE = "tests/sample_files/sample_bin.cfg" COMTRADE_SAMPLE_4_DAT_FILE = "tests/sample_files/sample_bin.dat" +# Two sample rates: 1000 Hz up to sample 4, then 500 Hz up to sample 6. +COMTRADE_SAMPLE_5_CFG = """STATION_NAME,EQUIPMENT,1999 +2,1A,1D +1, IA,,,A,1,0,0,-32768,32767,1,1,S +1, Trip,,,0 +60 +2 +1000,4 +500,6 +01/01/2000,10:30:00.000000 +01/01/2000,10:30:00.000000 +{format} +1 +""" + +COMTRADE_SAMPLE_5_DAT = "1,0,1,0\n2,1000,2,0\n3,2000,3,0\n4,3000,4,0\n5,5000,5,1\n6,7000,6,1\n" + +# Blank (non-critical) timestamps. +COMTRADE_SAMPLE_6_CFG = """STATION_NAME,EQUIPMENT,1999 +2,1A,1D +1, IA,,,A,1,0,0,-32768,32767,1,1,S +1, Trip,,,0 +60 +1 +1000,3 +01/01/2000,10:30:00.000000 +01/01/2000,10:30:00.000000 +ASCII +1 +""" + +COMTRADE_SAMPLE_6_DAT = "1,,1,0\n2,,2,0\n3,,3,1\n" + +# Status channels only. +COMTRADE_SAMPLE_7_CFG = """STATION_NAME,EQUIPMENT,2013 +2,0A,2D +1, Trip,,,0 +2, Close,,,0 +60 +1 +1000,3 +01/01/2000,10:30:00.000000 +01/01/2000,10:30:00.000000 +{format} +1 +""" + +# Analog channels with 0, 250 and -100 microseconds time skew. +COMTRADE_SAMPLE_8_CFG = """STATION_NAME,EQUIPMENT,1999 +3,3A,0D +1, IA,,,A,1,0,0,-32768,32767,1,1,S +2, IB,,,A,1,0,250,-32768,32767,1,1,S +3, IC,,,A,1,0,-100,-32768,32767,1,1,S +60 +1 +1000,3 +01/01/2000,10:30:00.000000 +01/01/2000,10:30:00.000000 +{format} +1 +""" + +COMTRADE_SAMPLE_8_DAT = "1,0,1,2,3\n2,1000,4,5,6\n3,2000,7,8,9\n" + class TestTimestamp(unittest.TestCase): """Test timestamp parsing.""" @@ -555,5 +620,190 @@ def test_should_use_double_when_specified(self): self.assertEqual(chan.typecode, "d") +class TestMultipleSampleRates(unittest.TestCase): + """ASCII DAT file with two sample rates.""" + def setUp(self): + self.comtrade = Comtrade(ignore_warnings=True) + self.comtrade.read(COMTRADE_SAMPLE_5_CFG.format(format="ASCII"), + COMTRADE_SAMPLE_5_DAT) + + def test_sample_rates(self): + self.assertEqual(self.comtrade.cfg.sample_rates, [[1000.0, 4], [500.0, 6]]) + + def test_time(self): + expected = [0.0, 0.001, 0.002, 0.003, 0.005, 0.007] + self.assertEqual(len(self.comtrade.time), len(expected)) + for t, t_expected in zip(self.comtrade.time, expected): + self.assertAlmostEqual(t, t_expected, places=6) + + +class TestMultipleSampleRatesBinary(TestMultipleSampleRates): + """Binary DAT file with two sample rates.""" + def setUp(self): + dat = b"".join(struct.pack(" list: + """Time, analog and status values lists of a record, with None for NaN + (since NaN != NaN).""" + return [[None if math.isnan(value) else value for value in values] + for values in [obj.time] + obj.analog + obj.status] + + +class TestAsciiBlockReading(unittest.TestCase): + """ASCII DAT file read in blocks of a few lines.""" + def test_block_sizes(self): + files = ("tests/sample_files/sample_ascii.cfg", + "tests/sample_files/sample_ascii_missing.dat") + expected = _channel_values(Comtrade(ignore_warnings=True).load(*files)) + # 10 values per line: 1 and 3 lines per block + for block_values in (1, 30): + with self.subTest(block_values=block_values), \ + mock.patch.object(comtrade, "_ASCII_BLOCK_VALUES", block_values): + obj = Comtrade(ignore_warnings=True).load(*files) + self.assertEqual(_channel_values(obj), expected) + + +class TestDatRowsCount(unittest.TestCase): + """DAT files with more or fewer rows than the CFG total samples.""" + use_numpy_arrays = False + + def read(self, dat_format, rows): + obj = Comtrade(ignore_warnings=True, + use_numpy_arrays=self.use_numpy_arrays) + if dat_format == "ASCII": + dat = "".join("{0},{1},{0},{2}\n".format(n, ts, status) + for n, ts, status in rows) + else: + dat = b"".join(struct.pack(" Date: Wed, 21 Oct 2026 14:00:03 -0300 Subject: [PATCH 3/4] Fix multiple sample rates, skew and DAT reading issues #48 #49 #29 - Compute sample times across multiple sample rate segments: each segment continues from the last sample of the previous one, instead of (n - 1) / samp with the segment's own rate (see #45) - Accept a blank (non-critical) timestamp in ASCII DAT files and compute it from the sample rate (#48) - Fix loading binary DAT files with status channels only (#49) - Add Comtrade.analog_time with each analog channel's time plus its time skew (#29) - Update README supported features table and add tests Fixes #48, fixes #49, fixes #29 Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01VXtokXPfpX3vJri639eqXo --- README.md | 7 +- comtrade.py | 58 ++++++++++---- tests/test_comtrade.py | 169 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 218 insertions(+), 16 deletions(-) diff --git a/README.md b/README.md index 0f9e390..b3c234e 100644 --- a/README.md +++ b/README.md @@ -73,6 +73,8 @@ The module function `load` must be called with the `CFF` file path specified. `Comtrade.analog` and `Comtrade.status` lists stores analog and status channel sample lists respectively. These can be accessed through zero-based indexes, i.e., `Comtrade.analog[0]`. The list `Comtrade.time` stores each sample time in seconds. +Analog channels may have a time skew (in microseconds) from the start of the sample period, which is not included in `Comtrade.time`. The list `Comtrade.analog_time` stores the time values of each analog channel, with its skew, in seconds, i.e., `Comtrade.analog_time[0]` for `Comtrade.analog[0]`. + More information can be accessed through `Comtrade.cfg` object, which stores data such as detailed channel information. Data of additional sections, such as HDR and INF, can be accessed through `hdr` and `inf` properties, respectively. @@ -121,17 +123,16 @@ Feel free to pull requests implementing one of these unsupported features or fix | Timestamp multiplication factor | | x | x | x | x | | Time code and local code | | | | x | x | | Time quality of the samples | | | | x | x | -| Analog channel time skew | | x | x | x | Partial | +| Analog channel time skew | | x | x | x | x | | Analog channel primary and secondary VT or CT ratio | | x | x | x | x | | Status channel phase and monitored circuit | | x | x | x | x | -| Multiple sample rates | x | x | x | x | Partial | +| Multiple sample rates | x | x | x | x | x | | Nanoseconds scale | | | x | x | x | ### Unsupported features * Nanoseconds time base within Python's `datetime` objects (such as `start_timestamp` and `trigger_timestamp` properties). It warns the user but doesn't use it, truncating the numbers. -* Use of multiple sample rates in time calculations for binary data. ### Additional settings and features diff --git a/comtrade.py b/comtrade.py index f228ef2..2462ecf 100644 --- a/comtrade.py +++ b/comtrade.py @@ -594,6 +594,9 @@ def __init__(self, **kwargs): ) self._analog_values = [] self._status_values = [] + # time values of each analog channel, with the channel time skew + # (computed when first accessed) + self._analog_time_values = None # Additional CFF data (or additional comtrade files) self._hdr = None @@ -659,6 +662,15 @@ def time(self) -> list: """Return the time values list.""" return self._time_values + @property + def analog_time(self) -> list: + """Return the time values list of each analog channel, i.e., the + sample time plus the channel time skew.""" + if self._analog_time_values is None: + self._analog_time_values = [self._get_skewed_time(channel.skew) + for channel in self._cfg.analog_channels] + return self._analog_time_values + @property def analog(self) -> list: """Return the analog channel values bidimensional list.""" @@ -794,6 +806,15 @@ def _dat_extract_data(self, dat) -> None: self._analog_values = dat.analog self._status_values = dat.status self._total_samples = dat.total_samples + self._analog_time_values = None + + def _get_skewed_time(self, skew: float): + """Return the time values shifted by a channel time skew (in microseconds).""" + skew_seconds = skew * 1E-6 + if _HAS_NUMPY and isinstance(self._time_values, numpy.ndarray): + return self._time_values + skew_seconds + return array.array(self._time_values.typecode, + [t + skew_seconds for t in self._time_values]) def to_dataframe(self, **kwargs) -> "pd.DataFrame": """Return a pandas DataFrame object with comtrade data.""" @@ -1111,6 +1132,17 @@ def _preallocate(self): steps = self._cfg.sample_rates[-1][1] # last samp field self._total_samples = steps + # first sample number and its time for each sample rate segment. + # Each sample is 1/samp after the previous one, where samp is the + # sample rate of the segment the sample belongs to. + self._rate_segments = [] + start_n, start_time = 1, 0.0 + for samp, endsamp in self._cfg.sample_rates: + self._rate_segments.append((samp, endsamp, start_n, start_time)) + if samp != 0.0: + start_time += (endsamp - start_n) / samp + start_n = endsamp + # analog and status count analog_count = self._cfg.analog_count status_count = self._cfg.status_count @@ -1134,24 +1166,24 @@ def _preallocate(self): self.status[i] = _preallocate_values("i", steps, self._use_numpy_arrays) - def _get_samp(self, n) -> float: - """Get the sampling rate for a sample n (1-based index).""" - # TODO: make tests. - last_sample_rate = 1.0 - for samp, endsamp in self._cfg.sample_rates: + def _get_rate_segment(self, n) -> tuple: + """Get the sampling rate, first sample number and its time for the + sample rate segment of a sample n (1-based index).""" + for samp, endsamp, start_n, start_time in self._rate_segments: if n <= endsamp: - return samp - return last_sample_rate + return samp, start_n, start_time + # samples beyond the last endsamp continue at the last sample rate + samp, _, start_n, start_time = self._rate_segments[-1] + return samp, start_n, start_time def _get_time(self, n: int, ts_value: float, time_base: float, time_multiplier: float): # TODO: add option to enforce dat file timestamp, when available. - # TODO: make tests. - sample_rate = self._get_samp(n) + sample_rate, start_n, start_time = self._get_rate_segment(n) if not self._cfg.timestamp_critical or ts_value == _TIMESTAMP_MISSING: # if the timestamp is missing, use calculated. if sample_rate != 0.0: - return (n - 1) / sample_rate + return start_time + (n - start_n) / sample_rate else: raise ComtradeError("Missing timestamp and no sample rate " "provided.") @@ -1208,8 +1240,8 @@ def parse(self, contents): values = line.strip().split(self.ASCII_SEPARATOR) n = int(values[0]) - # Read time - ts_val = float(values[1]) + # Read time (the non-critical timestamp may be left blank) + ts_val = float(values[1]) if values[1].strip() else _TIMESTAMP_MISSING ts = self._get_time(n, ts_val, time_base, time_mult) avalues = [value * a[i] + b[i] if not math.isnan(value) else float('nan') @@ -1265,7 +1297,7 @@ def get_reader_format(self, analog_channels, status_bytes): return self.STRUCT_FORMAT_ANALOG_ONLY.format(acount=analog_channels) else: # Status channels only. - return self.STRUCT_FORMAT_STATUS_ONLY.format(acount=dcount) + return self.STRUCT_FORMAT_STATUS_ONLY.format(dcount=dcount) def parse(self, contents): """Parse DAT binary file contents.""" diff --git a/tests/test_comtrade.py b/tests/test_comtrade.py index d58b7fc..97901c9 100644 --- a/tests/test_comtrade.py +++ b/tests/test_comtrade.py @@ -53,6 +53,70 @@ COMTRADE_SAMPLE_4_CFG_FILE = "tests/sample_files/sample_bin.cfg" COMTRADE_SAMPLE_4_DAT_FILE = "tests/sample_files/sample_bin.dat" +# Two sample rates: 1000 Hz up to sample 4, then 500 Hz up to sample 6. +COMTRADE_SAMPLE_5_CFG = """STATION_NAME,EQUIPMENT,1999 +2,1A,1D +1, IA,,,A,1,0,0,-32768,32767,1,1,S +1, Trip,,,0 +60 +2 +1000,4 +500,6 +01/01/2000,10:30:00.000000 +01/01/2000,10:30:00.000000 +{format} +1 +""" + +COMTRADE_SAMPLE_5_DAT = "1,0,1,0\n2,1000,2,0\n3,2000,3,0\n4,3000,4,0\n5,5000,5,1\n6,7000,6,1\n" + +# Blank (non-critical) timestamps. +COMTRADE_SAMPLE_6_CFG = """STATION_NAME,EQUIPMENT,1999 +2,1A,1D +1, IA,,,A,1,0,0,-32768,32767,1,1,S +1, Trip,,,0 +60 +1 +1000,3 +01/01/2000,10:30:00.000000 +01/01/2000,10:30:00.000000 +ASCII +1 +""" + +COMTRADE_SAMPLE_6_DAT = "1,,1,0\n2,,2,0\n3,,3,1\n" + +# Status channels only. +COMTRADE_SAMPLE_7_CFG = """STATION_NAME,EQUIPMENT,2013 +2,0A,2D +1, Trip,,,0 +2, Close,,,0 +60 +1 +1000,3 +01/01/2000,10:30:00.000000 +01/01/2000,10:30:00.000000 +{format} +1 +""" + +# Analog channels with 0, 250 and -100 microseconds time skew. +COMTRADE_SAMPLE_8_CFG = """STATION_NAME,EQUIPMENT,1999 +3,3A,0D +1, IA,,,A,1,0,0,-32768,32767,1,1,S +2, IB,,,A,1,0,250,-32768,32767,1,1,S +3, IC,,,A,1,0,-100,-32768,32767,1,1,S +60 +1 +1000,3 +01/01/2000,10:30:00.000000 +01/01/2000,10:30:00.000000 +{format} +1 +""" + +COMTRADE_SAMPLE_8_DAT = "1,0,1,2,3\n2,1000,4,5,6\n3,2000,7,8,9\n" + class TestTimestamp(unittest.TestCase): """Test timestamp parsing.""" @@ -555,5 +619,110 @@ def test_should_use_double_when_specified(self): self.assertEqual(chan.typecode, "d") +class TestMultipleSampleRates(unittest.TestCase): + """ASCII DAT file with two sample rates.""" + def setUp(self): + self.comtrade = Comtrade(ignore_warnings=True) + self.comtrade.read(COMTRADE_SAMPLE_5_CFG.format(format="ASCII"), + COMTRADE_SAMPLE_5_DAT) + + def test_sample_rates(self): + self.assertEqual(self.comtrade.cfg.sample_rates, [[1000.0, 4], [500.0, 6]]) + + def test_time(self): + expected = [0.0, 0.001, 0.002, 0.003, 0.005, 0.007] + self.assertEqual(len(self.comtrade.time), len(expected)) + for t, t_expected in zip(self.comtrade.time, expected): + self.assertAlmostEqual(t, t_expected, places=6) + + +class TestMultipleSampleRatesBinary(TestMultipleSampleRates): + """Binary DAT file with two sample rates.""" + def setUp(self): + dat = b"".join(struct.pack(" Date: Wed, 21 Oct 2026 14:05:16 -0300 Subject: [PATCH 4/4] Speed up reading of large DAT files #24 - Unpack each binary DAT field (sample number, timestamp, analog channel, 16-bit status group) for all rows at once with strided byte slices, replacing the struct row-by-row reader - Scale binary analog values and extract status bits with numpy when use_numpy_arrays is set - Read ASCII DAT files in blocks of lines (up to 1M values), converting each channel of a block at once - Binary DAT files ending in a partial row load their whole rows instead of failing with struct.error - Add tests for ASCII block reading, DAT files with more or fewer rows than the CFG samples and numpy arrays reading Loading the issue's sample (ASCII, 2212 analog channels, 10000 samples) drops from 3.6 s to 2.0 s. Binary files with 2212 channels load about 2x faster, and 5.5 to 10x faster with use_numpy_arrays. Fixes #24 Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01XUdNv59jhgaaRKFfZmRMwQ --- README.md | 4 +- comtrade.py | 209 ++++++++++++++++++++++------------------- tests/test_comtrade.py | 81 ++++++++++++++++ 3 files changed, 198 insertions(+), 96 deletions(-) diff --git a/README.md b/README.md index b3c234e..ddec2e9 100644 --- a/README.md +++ b/README.md @@ -146,7 +146,9 @@ in `Comtrade` object constructor: obj = Comtrade(use_numpy_arrays=True) ``` -It may improve performance for computations after loading data. +It may improve performance for computations after loading data. It also speeds up loading binary +DAT files (`BINARY`, `BINARY32` and `FLOAT32` formats), which are converted with vectorized numpy +operations. #### Replace missing data with `float('nan')` diff --git a/comtrade.py b/comtrade.py index 2462ecf..0f96d1a 100644 --- a/comtrade.py +++ b/comtrade.py @@ -27,10 +27,11 @@ import datetime as dt import errno import io +import itertools import math +import operator import os import re -import struct import sys from typing import Union import warnings @@ -71,6 +72,9 @@ # common separator character of data fields of CFG and ASCII DAT files _SEPARATOR = "," +# maximum number of ASCII DAT file values converted at once +_ASCII_BLOCK_VALUES = 1000000 + # timestamp regular expression _re_date = re.compile(r"([0-9]{1,2})/([0-9]{1,2})/([0-9]{2,4})") _re_time = re.compile(r"([0-9]{1,2}):([0-9]{2}):([0-9]{1,2})(.([0-9]{1,12}))?") @@ -103,6 +107,27 @@ def _preallocate_values(array_type, size, use_numpy_arrays): return array.array(array_type, [0]) * size +def _copy_values(target, start: int, values): + """Copy values into a preallocated values array, from the start index.""" + if isinstance(target, array.array): + values = array.array(target.typecode, values) + target[start:start + len(values)] = values + + +def _unpack_column(data, row_size: int, offset: int, typecode: str) -> array.array: + """Unpack the little-endian field of typecode at the offset (in bytes) of + every row of row_size bytes in data.""" + size = array.array(typecode).itemsize + column = bytearray(len(data) // row_size * size) + # gather each byte of the field from all rows at once + for k in range(size): + column[k::size] = data[offset + k::row_size] + values = array.array(typecode, column) + if sys.byteorder == "big": + values.byteswap() + return values + + def _prevent_null(str_value: str, value_type: type, default_value): if len(str_value.strip()) == 0: return default_value @@ -1191,9 +1216,6 @@ def _get_time(self, n: int, ts_value: float, time_base: float, # Use provided timestamp if it's not missing return ts_value * time_base * time_multiplier - def filter_missing(self, value) -> float: - return float(value) if value != self.DATA_MISSING else float('nan') - def parse(self, contents): """Virtual method, parse DAT file contents.""" pass @@ -1221,6 +1243,8 @@ def parse(self, contents): status_count = self._cfg.status_count time_mult = self._cfg.timemult time_base = self._cfg.time_base + missing = self.DATA_MISSING + nan = float('nan') # auxiliary vectors (channels gains and offsets) a = [x.a for x in self._cfg.analog_channels] @@ -1231,29 +1255,37 @@ def parse(self, contents): lines = contents.splitlines() else: lines = contents - - line_number = 0 - for line in lines: - line_number = line_number + 1 - if line_number > self._total_samples: + lines = iter(lines) + + # Lines are read in blocks and the values of each channel in a block + # are converted at once, which is much faster than one by one. + block_size = max(1, _ASCII_BLOCK_VALUES // (analog_count + status_count + 2)) + irow = 0 + while irow < self._total_samples: + rows = [line.strip().split(self.ASCII_SEPARATOR) for line in + itertools.islice(lines, min(block_size, self._total_samples - irow))] + if not rows: break - values = line.strip().split(self.ASCII_SEPARATOR) - n = int(values[0]) # Read time (the non-critical timestamp may be left blank) - ts_val = float(values[1]) if values[1].strip() else _TIMESTAMP_MISSING - ts = self._get_time(n, ts_val, time_base, time_mult) + _copy_values(self.time, irow, [ + self._get_time(int(values[0]), + float(values[1]) if values[1].strip() else _TIMESTAMP_MISSING, + time_base, time_mult) + for values in rows]) - avalues = [value * a[i] + b[i] if not math.isnan(value) else float('nan') - for i, value in enumerate(map(self.filter_missing, values[2:analog_count + 2]))] - svalues = [int(x) for x in values[len(values) - status_count:]] - - # store - self.time[line_number - 1] = ts for i in range(analog_count): - self.analog[i][line_number - 1] = avalues[i] + ai, bi = a[i], b[i] + _copy_values(self.analog[i], irow, + [ai * float(value) + bi if value != missing else nan + for value in map(operator.itemgetter(i + 2), rows)]) + + # status values are the last ones of each line for i in range(status_count): - self.status[i][line_number - 1] = svalues[i] + _copy_values(self.status[i], irow, + list(map(int, map(operator.itemgetter(i - status_count), rows)))) + + irow += len(rows) class _BinaryDatReader(_DatReader): @@ -1267,6 +1299,10 @@ def __init__(self, **kwargs): self.TIME_BYTES = 4 self.SAMPLE_NUMBER_BYTES = 4 + # array typecodes of the analog and the 16 bits status values + self.ANALOG_TYPECODE = "h" + self.STATUS_TYPECODE = "H" + if self._rev_year == REV_1991: self.DATA_MISSING = -1 # 0xFFFF else: @@ -1274,30 +1310,24 @@ def __init__(self, **kwargs): self.read_mode = "rb" - if struct.calcsize("L") == 4: - self.STRUCT_FORMAT = "LL {acount:d}h {dcount:d}H" - self.STRUCT_FORMAT_ANALOG_ONLY = "LL {acount:d}h" - self.STRUCT_FORMAT_STATUS_ONLY = "LL {dcount:d}H" - else: - self.STRUCT_FORMAT = "II {acount:d}h {dcount:d}H" - self.STRUCT_FORMAT_ANALOG_ONLY = "II {acount:d}h" - self.STRUCT_FORMAT_STATUS_ONLY = "II {dcount:d}H" - - def get_reader_format(self, analog_channels, status_bytes): - # Number of status fields of 2 bytes based on the total number of - # bytes. - dcount = math.floor(status_bytes / 2) - - # Check the file configuration - if int(status_bytes) > 0 and int(analog_channels) > 0: - return self.STRUCT_FORMAT.format(acount=analog_channels, - dcount=dcount) - elif int(analog_channels) > 0: - # Analog channels only. - return self.STRUCT_FORMAT_ANALOG_ONLY.format(acount=analog_channels) - else: - # Status channels only. - return self.STRUCT_FORMAT_STATUS_ONLY.format(dcount=dcount) + def _get_analog_values(self, yint, a: float, b: float): + """Return the analog values a * yint + b of the raw values yint, + with NaN for missing raw values.""" + if _HAS_NUMPY and self._use_numpy_arrays: + y = numpy.asarray(yint, dtype=numpy.float64) + is_missing = y == self.DATA_MISSING + y = a * y + b + y[is_missing] = numpy.nan + return y + missing = self.DATA_MISSING + nan = float('nan') + return [a * value + b if value != missing else nan for value in yint] + + def _get_status_values(self, group, chnindex: int): + """Return the status values of the bit chnindex of a group of 16 bits.""" + if _HAS_NUMPY and self._use_numpy_arrays: + return (numpy.asarray(group) >> chnindex) & 1 + return [(word >> chnindex) & 1 for word in group] def parse(self, contents): """Parse DAT binary file contents.""" @@ -1313,47 +1343,48 @@ def parse(self, contents): dbytes = self.STATUS_BYTES * math.ceil(schannel / 16.0) groups_of_16bits = math.floor(dbytes / self.STATUS_BYTES) - # Struct format. - row_reader = struct.Struct(self.get_reader_format(achannels, dbytes)) - - # Row reading function. if isinstance(contents, io.TextIOBase) or \ isinstance(contents, io.BufferedIOBase): # Read all buffer contents contents = contents.read() - for irow, values in enumerate(row_reader.iter_unpack(contents)): - # Sample number - n = values[0] - # Time stamp - ts_val = values[1] - ts = self._get_time(n, ts_val, time_base, time_mult) - - if irow >= self.total_samples: - break - self.time[irow] = ts - - # Extract analog channel values. - for ichannel in range(achannels): - yint = self.filter_missing(values[ichannel + 2]) - y = a[ichannel] * yint + b[ichannel] if not math.isnan(yint) else float('nan') - self.analog[ichannel][irow] = y - - # Extract status channel values. - for igroup in range(groups_of_16bits): - group = values[achannels + 2 + igroup] - - # for each group of 16 bits, extract the status channels - maxchn = min([(igroup + 1) * 16, schannel]) - for ichannel in range(igroup * 16, maxchn): - chnindex = ichannel - igroup * 16 - mask = int('0b01', 2) << chnindex - extract = (group & mask) >> chnindex - - self.status[ichannel][irow] = extract - - # Get the next row - irow += 1 + # Each row has the sample number, the time stamp, the analog values + # and the groups of 16 status bits. Only whole rows are read, up to + # the total samples. + row_size = self.SAMPLE_NUMBER_BYTES + self.TIME_BYTES + \ + achannels * self.ANALOG_BYTES + dbytes + rows = min(len(contents) // row_size, self.total_samples) + contents = contents[:rows * row_size] + + # Each field is unpacked from all rows at once, which is much faster + # than unpacking the rows one by one. + sample_numbers = _unpack_column(contents, row_size, 0, "I") + timestamps = _unpack_column(contents, row_size, + self.SAMPLE_NUMBER_BYTES, "I") + _copy_values(self.time, 0, [self._get_time(n, ts_val, time_base, time_mult) + for n, ts_val in zip(sample_numbers, timestamps)]) + + # Extract analog channel values. + offset = self.SAMPLE_NUMBER_BYTES + self.TIME_BYTES + for ichannel in range(achannels): + yint = _unpack_column(contents, row_size, offset, + self.ANALOG_TYPECODE) + _copy_values(self.analog[ichannel], 0, + self._get_analog_values(yint, a[ichannel], b[ichannel])) + offset += self.ANALOG_BYTES + + # Extract status channel values. + for igroup in range(groups_of_16bits): + group = _unpack_column(contents, row_size, offset, + self.STATUS_TYPECODE) + + # for each group of 16 bits, extract the status channels + maxchn = min((igroup + 1) * 16, schannel) + for ichannel in range(igroup * 16, maxchn): + chnindex = ichannel - igroup * 16 + _copy_values(self.status[ichannel], 0, + self._get_status_values(group, chnindex)) + offset += self.STATUS_BYTES class _Binary32DatReader(_BinaryDatReader): @@ -1363,13 +1394,7 @@ def __init__(self, **kwargs): # Call the initialization for the inherited class super().__init__(**kwargs) self.ANALOG_BYTES = 4 - - if struct.calcsize("L") == 4: - self.STRUCT_FORMAT = "LL {acount:d}l {dcount:d}H" - self.STRUCT_FORMAT_ANALOG_ONLY = "LL {acount:d}l" - else: - self.STRUCT_FORMAT = "II {acount:d}i {dcount:d}H" - self.STRUCT_FORMAT_ANALOG_ONLY = "II {acount:d}i" + self.ANALOG_TYPECODE = "i" # maximum negative value self.DATA_MISSING = -2147483648 # 0x80000000 @@ -1382,13 +1407,7 @@ def __init__(self, **kwargs): # Call the initialization for the inherited class super().__init__(**kwargs) self.ANALOG_BYTES = 4 - - if struct.calcsize("L") == 4: - self.STRUCT_FORMAT = "LL {acount:d}f {dcount:d}H" - self.STRUCT_FORMAT_ANALOG_ONLY = "LL {acount:d}f" - else: - self.STRUCT_FORMAT = "II {acount:d}f {dcount:d}H" - self.STRUCT_FORMAT_ANALOG_ONLY = "II {acount:d}f" + self.ANALOG_TYPECODE = "f" # Maximum negative value self.DATA_MISSING = sys.float_info.min diff --git a/tests/test_comtrade.py b/tests/test_comtrade.py index 97901c9..28ac3b1 100644 --- a/tests/test_comtrade.py +++ b/tests/test_comtrade.py @@ -3,6 +3,7 @@ import os import struct import unittest +from unittest import mock import comtrade from comtrade import Comtrade @@ -724,5 +725,85 @@ class TestAnalogChannelSkewNumpy(TestAnalogChannelSkew): use_numpy_arrays = True +def _channel_values(obj) -> list: + """Time, analog and status values lists of a record, with None for NaN + (since NaN != NaN).""" + return [[None if math.isnan(value) else value for value in values] + for values in [obj.time] + obj.analog + obj.status] + + +class TestAsciiBlockReading(unittest.TestCase): + """ASCII DAT file read in blocks of a few lines.""" + def test_block_sizes(self): + files = ("tests/sample_files/sample_ascii.cfg", + "tests/sample_files/sample_ascii_missing.dat") + expected = _channel_values(Comtrade(ignore_warnings=True).load(*files)) + # 10 values per line: 1 and 3 lines per block + for block_values in (1, 30): + with self.subTest(block_values=block_values), \ + mock.patch.object(comtrade, "_ASCII_BLOCK_VALUES", block_values): + obj = Comtrade(ignore_warnings=True).load(*files) + self.assertEqual(_channel_values(obj), expected) + + +class TestDatRowsCount(unittest.TestCase): + """DAT files with more or fewer rows than the CFG total samples.""" + use_numpy_arrays = False + + def read(self, dat_format, rows): + obj = Comtrade(ignore_warnings=True, + use_numpy_arrays=self.use_numpy_arrays) + if dat_format == "ASCII": + dat = "".join("{0},{1},{0},{2}\n".format(n, ts, status) + for n, ts, status in rows) + else: + dat = b"".join(struct.pack("