Comparing Codecs

This example compresses a DAS recording of an event with each ComDAS codec at several settings, then compares the waterfall plots, the compression ratios, and the reconstruction error.

import matplotlib.pyplot as plt
import numpy as np
import dascore as dc

from comdas import SVDCodec, WaveletCodec, compress_patch

patch = dc.get_example_patch("example_event_2")
print(patch.dims, patch.shape, patch.dtype)
print(f"uncompressed size: {patch.data.nbytes / 1e6:.2f} MB")

Downloading file 'example_dasdae_event_1.h5' from 'https://github.com/dasdae/test_data/raw/master/das/example_dasdae_event_1.h5' to '/home/runner/.cache/dascore/0.0.0'.

('distance', 'time') (601, 1001) float64
uncompressed size: 4.81 MB

The original data

SCALE = 0.5

fig, ax = plt.subplots(figsize=(7, 4))
patch.viz.waterfall(ax=ax, scale=SCALE)
ax.set_title("Original")
plt.show()

Compress with several codecs

Each compressed patch is an ordinary DASCore patch. The codec and compressed payload are available from patch.data, which lets us compute the compression ratio. We also measure the relative error

\[\varepsilon = \frac{\lVert A - \hat A \rVert_F}{\lVert A \rVert_F}.\]

codecs = {
    "SVD, rank=5": SVDCodec(rank=5),
    "SVD, rank=25": SVDCodec(rank=25),
    "SVD, energy=0.95": SVDCodec(energy=0.95),
    "Wavelet db4, keep 2%": WaveletCodec(wavelet="db4", keep_fraction=0.02),
    "Wavelet db4, keep 10%": WaveletCodec(wavelet="db4", keep_fraction=0.10),
    "Wavelet sym8, keep 5%": WaveletCodec(wavelet="sym8", keep_fraction=0.05),
}


def summarize(original, compressed):
    """Return (compression ratio, relative error) for a compressed patch."""
    ratio = compressed.data.codec.compression_ratio(compressed.data.payload)
    a = np.asarray(original.data)
    error = np.linalg.norm(a - np.asarray(compressed.data)) / np.linalg.norm(a)
    return ratio, error


compressed = {
    label: compress_patch(patch, codec) for label, codec in codecs.items()
}

print(f"{'codec':<24}{'ratio':>10}{'rel. error':>12}")
for label, cpatch in compressed.items():
    ratio, error = summarize(patch, cpatch)
    print(f"{label:<24}{ratio:>9.1f}x{error:>12.3f}")

codec                        ratio  rel. error
SVD, rank=5                  75.1x       0.687
SVD, rank=25                 15.0x       0.312
SVD, energy=0.95             10.7x       0.221
Wavelet db4, keep 2%         32.2x       0.308
Wavelet db4, keep 10%         6.4x       0.129
Wavelet sym8, keep 5%        12.4x       0.204

Waterfall plots

Compressed patches plot with the same patch.viz.waterfall call as the original.

fig, axes = plt.subplots(3, 2, figsize=(12, 11), sharex=True, sharey=True)
for ax, (label, cpatch) in zip(axes.flat, compressed.items()):
    ratio, error = summarize(patch, cpatch)
    cpatch.viz.waterfall(ax=ax, scale=SCALE)
    ax.set_title(f"{label}\nratio {ratio:.1f}x, error {error:.3f}")
fig.tight_layout()
plt.show()

What was lost

Plotting the residual (original minus reconstruction) on the same colour scale as the original shows which features each codec discards.

vmax = SCALE * np.abs(patch.data).max()

fig, axes = plt.subplots(3, 2, figsize=(12, 11), sharex=True, sharey=True)
for ax, (label, cpatch) in zip(axes.flat, compressed.items()):
    residual = patch.new(data=np.asarray(patch.data) - np.asarray(cpatch.data))
    residual.viz.waterfall(ax=ax, scale=vmax, scale_type="absolute")
    ax.set_title(f"Residual: {label}")
fig.tight_layout()
plt.show()

Compression ratio vs. error

Sweeping each codec’s main setting shows the trade-off between size and fidelity for this recording.

sweeps = {
    "SVD": [SVDCodec(rank=k) for k in (1, 2, 5, 10, 20, 40, 80)],
    "Wavelet db4": [
        WaveletCodec(wavelet="db4", keep_fraction=f)
        for f in (0.005, 0.01, 0.02, 0.05, 0.1, 0.2)
    ],
    "Wavelet sym8": [
        WaveletCodec(wavelet="sym8", keep_fraction=f)
        for f in (0.005, 0.01, 0.02, 0.05, 0.1, 0.2)
    ],
}

fig, ax = plt.subplots(figsize=(7, 4.5))
for name, sweep in sweeps.items():
    points = [summarize(patch, compress_patch(patch, codec)) for codec in sweep]
    ratios, errors = zip(*points)
    ax.plot(ratios, errors, marker="o", label=name)
ax.set_xscale("log")
ax.set_xlabel("compression ratio")
ax.set_ylabel("relative error")
ax.grid(True, which="both", alpha=0.3)
ax.legend()
plt.show()

Saving to disk

The same codecs are used when writing files. Here each version is written to its own file and read back with plain DASCore.

import tempfile
from pathlib import Path

from comdas import write_compressed

with tempfile.TemporaryDirectory() as tmp:
    raw_path = Path(tmp) / "original.h5"
    dc.write(patch, raw_path, "DASDAE")
    print(f"{'original (DASDAE)':<24}{raw_path.stat().st_size / 1e6:>8.2f} MB")

    for i, (label, codec) in enumerate(codecs.items()):
        path = Path(tmp) / f"compressed_{i}.h5"
        write_compressed(patch, path, codec)
        restored = dc.spool(path)[0]
        assert restored.shape == patch.shape
        print(f"{label:<24}{path.stat().st_size / 1e6:>8.2f} MB")

/opt/hostedtoolcache/Python/3.12.14/x64/lib/python3.12/site-packages/tables/path.py:146: NaturalNameWarning: object name is not a valid Python identifier: 'DAS_________0.0__0.10000000000000286'; it does not match the pattern ``^[a-zA-Z_][a-zA-Z0-9_]*$``; you will not be able to use natural naming to access this object; using ``getattr()`` will still work, though
  check_attribute_name(name)

original (DASDAE)           4.83 MB
SVD, rank=5                 0.09 MB

SVD, rank=25                0.35 MB
SVD, energy=0.95            0.47 MB
Wavelet db4, keep 2%        0.17 MB

Wavelet db4, keep 10%       0.77 MB
Wavelet sym8, keep 5%       0.41 MB