Skip to content

dyned.py

DyNED encoder + DyNEDc compressor - 1,739 lines. View on GitHub (dyned.py).

"""
DyNED - Dynamic Neural Encoding/Decoding

DyNED: Spike encoder using sigma-delta quantization of frequency-domain signals.
DyNEDc: Compressor optimized for spike trains (v1: RLE+Huffman, v2: Golomb coding).

Both techniques form a universal encoding pipeline for temporal (speech) and
static (image) signals.

Usage:
    from dyned import sigma_delta_quantisation, generate_step_signal
    from dyned import DyNEDcCompressor, compression_summary, print_compression_table
"""

import heapq
import math
from collections import Counter, defaultdict

import numpy as np
import torch

try:
    from numba import njit

    @njit(cache=True)
    def _sigma_delta_loop(normalized, levels):
        """JIT-compiled sigma-delta feedback loop (~100x faster than Python)."""
        n = len(normalized)
        quantized = np.zeros(n, dtype=np.float64)
        errors = np.zeros(n, dtype=np.float64)
        error = 0.0
        for i in range(n):
            sample_with_error = normalized[i] + error
            q = round(sample_with_error * (levels - 1))
            if q < 0:
                q = 0
            elif q > levels - 1:
                q = levels - 1
            q_norm = q / (levels - 1)
            quantized[i] = q_norm
            error = sample_with_error - q_norm
            errors[i] = error
        return quantized, errors

    @njit(cache=True)
    def _dyned_encode_2d(log_dense, levels):
        """Fully JIT-compiled 2D DyNED encode: sigma-delta + step signal.

        Processes all time bins in compiled code  -  no Python loop overhead.
        Produces identical output to the per-bin sigma_delta_quantisation +
        generate_step_signal loop but ~50-100x faster.
        """
        n_channels, n_bins = log_dense.shape
        spike_train = np.zeros((n_channels, n_bins), dtype=np.float32)

        for t in range(n_bins):
            # Find min/max for this time bin
            data_min = log_dense[0, t]
            data_max = log_dense[0, t]
            for i in range(1, n_channels):
                val = log_dense[i, t]
                if val < data_min:
                    data_min = val
                if val > data_max:
                    data_max = val

            data_range = data_max - data_min
            if data_range == 0.0:
                continue  # All same value -> no transitions -> all zeros

            # Sigma-delta quantisation (inlined)
            prev_quantized = 0.0
            error = 0.0
            for i in range(n_channels):
                normalized = (log_dense[i, t] - data_min) / data_range
                sample_with_error = normalized + error
                q = round(sample_with_error * (levels - 1))
                if q < 0:
                    q = 0
                elif q > levels - 1:
                    q = levels - 1
                q_norm = q / (levels - 1)
                quantized_val = q_norm * data_range + data_min
                error = sample_with_error - q_norm

                # Step signal: 1 where quantized value changes
                if i > 0 and quantized_val != prev_quantized:
                    spike_train[i, t] = 1.0
                prev_quantized = quantized_val

        return spike_train

    @njit(cache=True)
    def _dyned_quantise_2d(data, levels):
        """Fully JIT-compiled 2D sigma-delta quantisation returning quantised values + step signal.

        Used by speech scripts that need the quantised output (not just the step signal).
        """
        n_feat, n_time = data.shape
        quantised_out = np.zeros((n_feat, n_time), dtype=np.float32)
        step_out = np.zeros((n_feat, n_time), dtype=np.float32)

        for t in range(n_time):
            # Find min/max
            data_min = data[0, t]
            data_max = data[0, t]
            for i in range(1, n_feat):
                val = data[i, t]
                if val < data_min:
                    data_min = val
                if val > data_max:
                    data_max = val

            data_range = data_max - data_min
            if data_range == 0.0:
                # All same value  -  quantised = original, step = 0
                for i in range(n_feat):
                    quantised_out[i, t] = np.float32(data[i, t])
                continue

            prev_quantized = 0.0
            error = 0.0
            for i in range(n_feat):
                normalized = (data[i, t] - data_min) / data_range
                sample_with_error = normalized + error
                q = round(sample_with_error * (levels - 1))
                if q < 0:
                    q = 0
                elif q > levels - 1:
                    q = levels - 1
                q_norm = q / (levels - 1)
                quantized_val = q_norm * data_range + data_min
                error = sample_with_error - q_norm

                quantised_out[i, t] = np.float32(quantized_val)
                if i > 0 and quantized_val != prev_quantized:
                    step_out[i, t] = 1.0
                prev_quantized = quantized_val

        return quantised_out, step_out

    HAS_NUMBA = True
except ImportError:
    HAS_NUMBA = False

try:
    import plotly.graph_objects as go

    HAS_PLOTLY = True
except ImportError:
    HAS_PLOTLY = False


# =============================================================================
# DyNED Encoder - Sigma-Delta Spike Encoding
# =============================================================================


def sigma_delta_quantisation(data, levels=256):
    """
    Sigma-delta quantization of signal data.

    Converts a continuous signal into a quantized representation using
    feedback-driven error diffusion. Uses numba JIT if available.

    Args:
        data: Input signal (numpy array or torch Tensor).
        levels: Number of quantization levels (default 256).

    Returns:
        quantized: Quantized signal (numpy array, same scale as input).
        errors: Quantization error at each sample.
    """
    if hasattr(data, "cpu"):
        data = data.cpu()
    data = np.asarray(data, dtype=np.float64)

    data_min, data_max = data.min(), data.max()
    data_range = data_max - data_min
    if data_range == 0:
        return data.copy(), np.zeros_like(data)

    normalized = (data - data_min) / data_range

    if HAS_NUMBA:
        quantized, errors = _sigma_delta_loop(normalized, levels)
    else:
        quantized = np.zeros_like(normalized)
        error = 0.0
        errors = np.zeros_like(normalized)
        for i in range(len(normalized)):
            sample_with_error = normalized[i] + error
            q = np.clip(np.round(sample_with_error * (levels - 1)), 0, levels - 1)
            q_norm = q / (levels - 1)
            quantized[i] = q_norm
            error = sample_with_error - q_norm
            errors[i] = error

    # Denormalize back to original scale
    quantized = quantized * data_range + data_min
    return quantized, errors


def generate_step_signal(quantised_data):
    """Generate binary step signal from quantized data (1 where value changes)."""
    step_signal = np.zeros_like(quantised_data, dtype=bool)
    step_signal[1:] = np.diff(quantised_data) != 0
    return step_signal


def dyned_encode_dense(dense, dyned_levels=256):
    """Apply per-time-bin DyNED sigma-delta quantisation to a dense 2D array.

    Numba-accelerated: processes the entire [n_channels, n_bins] matrix in one
    compiled call instead of looping through Python per time bin.

    Args:
        dense: float32 array [n_channels, n_bins]
        dyned_levels: number of sigma-delta quantisation levels

    Returns:
        spike_train: float32 tensor [n_channels, n_bins]
    """
    # log1p on float32 first (matches original per-bin behaviour), then promote
    log_dense = np.log1p(np.asarray(dense, dtype=np.float32)).astype(np.float64)

    if HAS_NUMBA:
        spike_train = _dyned_encode_2d(log_dense, dyned_levels)
    else:
        # Pure-Python fallback (same logic, slower)
        n_channels, n_bins = log_dense.shape
        spike_train = np.zeros((n_channels, n_bins), dtype=np.float32)
        for t in range(n_bins):
            frame = log_dense[:, t]
            quantised, _ = sigma_delta_quantisation(frame, levels=dyned_levels)
            step_signal = generate_step_signal(quantised)
            spike_train[:, t] = step_signal.astype(np.float32)

    return torch.from_numpy(spike_train)


def dyned_encode_features(features, levels=256):
    """Sigma-delta encode pre-computed features, returning step signal only.

    Unlike dyned_encode_dense(), no log1p is applied  -  input should already
    be in feature space (e.g. log STFT magnitude, mel spectrogram, etc.).

    Args:
        features: float32 array [n_feat, n_time]
        levels: number of sigma-delta quantisation levels

    Returns:
        spike_train: float32 tensor [n_feat, n_time]
    """
    data = np.asarray(features, dtype=np.float64)

    if HAS_NUMBA:
        spike_train = _dyned_encode_2d(data, levels)
    else:
        n_feat, n_time = data.shape
        spike_train = np.zeros((n_feat, n_time), dtype=np.float32)
        for t in range(n_time):
            q, _ = sigma_delta_quantisation(data[:, t], levels=levels)
            s = generate_step_signal(q)
            spike_train[:, t] = s.astype(np.float32)

    return torch.from_numpy(spike_train)


def dyned_quantise_2d(data, levels=256, boost=False):
    """Sigma-delta quantise a 2D feature array, returning quantised values.

    Numba-accelerated. Used by speech scripts that need the quantised output
    (not just the step signal). Optionally applies boost: quantised * (1 + step).

    Args:
        data: float32 array [n_feat, n_time] (already in feature space, e.g. mel/MFCC)
        levels: number of sigma-delta quantisation levels
        boost: if True, return quantised * (1 + step_signal)

    Returns:
        result: float32 tensor [n_feat, n_time]
    """
    data_f64 = np.asarray(data, dtype=np.float64)

    if HAS_NUMBA:
        quantised, step = _dyned_quantise_2d(data_f64, levels)
    else:
        n_feat, n_time = data_f64.shape
        quantised = np.zeros((n_feat, n_time), dtype=np.float32)
        step = np.zeros((n_feat, n_time), dtype=np.float32)
        for t in range(n_time):
            frame = data_f64[:, t]
            q, _ = sigma_delta_quantisation(frame, levels=levels)
            s = generate_step_signal(q)
            quantised[:, t] = q.astype(np.float32)
            step[:, t] = s.astype(np.float32)

    if boost:
        return torch.from_numpy(quantised * (1.0 + step))
    return torch.from_numpy(quantised)


def identify_peaks_and_troughs(data):
    """Find peaks and troughs in data using scipy."""
    from scipy.signal import find_peaks

    peaks, _ = find_peaks(data)
    troughs, _ = find_peaks(-data)
    return peaks, troughs


# =============================================================================
# Compression Algorithms (for comparison benchmarks)
# =============================================================================


class RunLengthEncoding:
    """Run-Length Encoding with variable-length binary count encoding."""

    def compress(self, binary_data: str) -> tuple[str, dict]:
        if not binary_data:
            return "", {"original_length": 0, "compressed_length": 0}

        compressed = []
        count = 1
        current = binary_data[0]

        for i in range(1, len(binary_data)):
            if binary_data[i] == current:
                count += 1
            else:
                count_binary = format(count, "b")
                count_size = format(len(count_binary) - 1, "03b")
                compressed.append(count_size + count_binary + current)
                current = binary_data[i]
                count = 1

        count_binary = format(count, "b")
        count_size = format(len(count_binary) - 1, "03b")
        compressed.append(count_size + count_binary + current)

        compressed_str = "".join(compressed)
        compression_info = {
            "original_length": len(binary_data),
            "compressed_length": len(compressed_str),
            "compression_ratio": len(compressed_str) / len(binary_data),
        }
        return compressed_str, compression_info

    def decompress(self, compressed_data: str) -> str:
        if not compressed_data:
            return ""

        decompressed = []
        i = 0
        while i < len(compressed_data):
            count_size = int(compressed_data[i : i + 3], 2) + 1
            i += 3
            count = int(compressed_data[i : i + count_size], 2)
            i += count_size
            value = compressed_data[i]
            i += 1
            decompressed.append(value * count)

        return "".join(decompressed)


class LZWCompressor:
    """LZW compression for binary strings."""

    def compress(self, data: str) -> list[int]:
        dictionary = {str(i): i for i in range(2)}
        next_code = len(dictionary)
        current = ""
        compressed = []

        for char in data:
            current_plus_char = current + char
            if current_plus_char in dictionary:
                current = current_plus_char
            else:
                compressed.append(dictionary[current])
                dictionary[current_plus_char] = next_code
                next_code += 1
                current = char

        if current:
            compressed.append(dictionary[current])

        return compressed

    def decompress(self, compressed: list[int]) -> str:
        dictionary = {i: str(i) for i in range(2)}
        next_code = len(dictionary)

        result = [str(compressed[0])]
        current = result[0]

        for code in compressed[1:]:
            if code in dictionary:
                entry = dictionary[code]
            elif code == next_code:
                entry = current + current[0]
            else:
                raise ValueError("Invalid compressed data")

            result.append(entry)
            dictionary[next_code] = current + entry[0]
            next_code += 1
            current = entry

        return "".join(result)


class HuffmanNode:
    """Node in a Huffman tree."""

    def __init__(self, char, freq):
        self.char = char
        self.freq = freq
        self.left = None
        self.right = None

    def __lt__(self, other):
        return self.freq < other.freq


class HuffmanCompressor:
    """Huffman coding compression with chunked input."""

    def __init__(self, chunk_size=4):
        self.codes = {}
        self.chunk_size = chunk_size

    def _make_codes(self, node, code=""):
        if node is None:
            return
        if node.char is not None:
            self.codes[node.char] = code or "0"
            return
        self._make_codes(node.left, code + "0")
        self._make_codes(node.right, code + "1")

    def compress(self, binary_data: str) -> tuple[str, dict]:
        if not binary_data:
            return "", {}

        chunks = [binary_data[i : i + self.chunk_size] for i in range(0, len(binary_data), self.chunk_size)]
        freq = Counter(chunks)
        heap = [HuffmanNode(chunk, count) for chunk, count in freq.items()]
        heapq.heapify(heap)

        while len(heap) > 1:
            left = heapq.heappop(heap)
            right = heapq.heappop(heap)
            internal = HuffmanNode(None, left.freq + right.freq)
            internal.left = left
            internal.right = right
            heapq.heappush(heap, internal)

        self.codes = {}
        if heap:
            self._make_codes(heap[0])

        compressed = "".join(self.codes[chunk] for chunk in chunks)
        compression_info = {
            "codes": self.codes,
            "original_length": len(binary_data),
            "compressed_length": len(compressed),
            "compression_ratio": len(compressed) / len(binary_data),
        }
        return compressed, compression_info

    def decompress(self, compressed_data: str, codes: dict) -> str:
        if not compressed_data or not codes:
            return ""

        reverse_codes = {v: k for k, v in codes.items()}
        decoded = []
        buffer = ""

        for bit in compressed_data:
            buffer += bit
            if buffer in reverse_codes:
                decoded.append(reverse_codes[buffer])
                buffer = ""

        return "".join(decoded)


class DeltaCompressor:
    """Delta encoding compression for integer sequences (e.g. spike positions)."""

    def compress(self, indices: list[int]) -> list[int]:
        if not indices:
            return []
        deltas = [indices[0]]
        for i in range(1, len(indices)):
            deltas.append(indices[i] - indices[i - 1])
        return deltas

    def decompress(self, deltas: list[int]) -> list[int]:
        if not deltas:
            return []
        values = [deltas[0]]
        for i in range(1, len(deltas)):
            values.append(values[-1] + deltas[i])
        return values


class BWTransform:
    """Burrows-Wheeler Transform (preprocessing for compression)."""

    # BWT is O(n^2) memory; limit to prevent OOM on long signals
    MAX_LENGTH = 4096

    def transform(self, data: str) -> tuple[str, int]:
        if not data:
            return "", 0
        if len(data) > self.MAX_LENGTH:
            data = data[: self.MAX_LENGTH]
        rotations = [data[i:] + data[:i] for i in range(len(data))]
        sorted_rot = sorted(rotations)
        I = sorted_rot.index(data)
        return "".join(r[-1] for r in sorted_rot), I


class LZ77Compressor:
    """LZ77 sliding-window compression."""

    def __init__(self, window_size=8):
        self.window_size = window_size

    def compress(self, data: str) -> list:
        compressed = []
        i = 0
        while i < len(data):
            longest_match = ""
            match_pos = -1
            for j in range(max(0, i - self.window_size), i):
                for length in range(len(data[j : min(j + self.window_size, len(data))]), 0, -1):
                    pattern = data[j : j + length]
                    if data[i : i + length] == pattern:
                        if length > len(longest_match):
                            longest_match = pattern
                            match_pos = j
                        break

            if match_pos >= 0 and len(longest_match) > 2:
                compressed.append((match_pos, len(longest_match)))
                i += len(longest_match)
            else:
                compressed.append(data[i])
                i += 1

        return compressed


# =============================================================================
# DyNEDc Compressor - Hybrid RLE + Huffman for Spike Trains
# =============================================================================


class DyNEDcCompressor:
    """
    DyNEDc (Dynamic Neural Encoding/Decoding Compression).

    Two-stage compression for binary spike trains:
      1. Extract alternating run lengths (no binary overhead per run)
      2. Huffman-encode the run-length values (common lengths get short codes)

    Designed for sparse spike trains from DyNED encoding where most of the
    signal is 0 with occasional 1s.
    """

    def __init__(self):
        self.huffman_codes = {}
        self._start_value = "0"

    def _extract_runs(self, binary_data):
        """Extract alternating run lengths and the starting value."""
        runs = []
        current = binary_data[0]
        count = 1
        for c in binary_data[1:]:
            if c == current:
                count += 1
            else:
                runs.append(count)
                current = c
                count = 1
        runs.append(count)
        return runs, binary_data[0]

    def _build_huffman_codes(self, symbols):
        """Build Huffman codes for a list of string symbols."""
        freq = Counter(symbols)
        if len(freq) == 1:
            self.huffman_codes = {list(freq.keys())[0]: "0"}
            return

        heap = [[count, [sym, ""]] for sym, count in freq.items()]
        heapq.heapify(heap)

        while len(heap) > 1:
            lo = heapq.heappop(heap)
            hi = heapq.heappop(heap)
            for pair in lo[1:]:
                pair[1] = "0" + pair[1]
            for pair in hi[1:]:
                pair[1] = "1" + pair[1]
            heapq.heappush(heap, [lo[0] + hi[0]] + lo[1:] + hi[1:])

        if heap:
            self.huffman_codes = {sym: code for sym, code in heap[0][1:]}
        else:
            self.huffman_codes = {}

    def compress(self, binary_data: str) -> tuple[str, dict]:
        """Compress binary spike train."""
        runs, self._start_value = self._extract_runs(binary_data)

        # Huffman-encode run lengths (as string symbols)
        symbols = [str(r) for r in runs]
        self._build_huffman_codes(symbols)
        compressed = "".join(self.huffman_codes[s] for s in symbols)

        compression_info = {
            "original_length": len(binary_data),
            "num_runs": len(runs),
            "compressed_length": len(compressed),
            "compression_ratio": len(compressed) / len(binary_data),
        }
        return compressed, compression_info

    def decompress(self, compressed_data: str) -> str:
        """Decompress DyNEDc data back to binary spike train."""
        reverse = {v: k for k, v in self.huffman_codes.items()}

        # Decode Huffman to get run lengths
        runs = []
        buffer = ""
        for bit in compressed_data:
            buffer += bit
            if buffer in reverse:
                runs.append(int(reverse[buffer]))
                buffer = ""

        # Reconstruct binary string from alternating runs
        result = []
        current = self._start_value
        for length in runs:
            result.append(current * length)
            current = "1" if current == "0" else "0"

        return "".join(result)


# =============================================================================
# Golomb Coding + DyNEDc V2 (Position-Based Golomb)
# =============================================================================


class GolombCoder:
    """Golomb coding for non-negative integers.

    Provably optimal prefix-free code for geometric distributions,
    which is the distribution of inter-spike gaps in sparse spike trains.
    """

    @staticmethod
    def optimal_m(p):
        """Optimal Golomb parameter for geometric distribution with spike density p."""
        if p <= 0 or p >= 1:
            return 1
        return max(1, int(math.ceil(-1.0 / math.log2(1.0 - p))))

    @staticmethod
    def encode(n, m):
        """Encode non-negative integer n with Golomb parameter m."""
        q, r = divmod(n, m)
        # Quotient in unary: q zeros + "1"
        code = "0" * q + "1"
        if m == 1:
            return code
        b = int(math.ceil(math.log2(m)))
        cutoff = (1 << b) - m
        if r < cutoff:
            if b > 1:
                code += format(r, f"0{b - 1}b")
        else:
            code += format(r + cutoff, f"0{b}b")
        return code

    @staticmethod
    def decode(bits, pos, m):
        """Decode one Golomb-coded integer. Returns (value, new_pos)."""
        # Read unary quotient (count zeros until "1")
        q = 0
        while pos < len(bits) and bits[pos] == "0":
            q += 1
            pos += 1
        pos += 1  # skip the "1"
        if m == 1:
            return q, pos
        b = int(math.ceil(math.log2(m)))
        cutoff = (1 << b) - m
        r_bits = b - 1 if b > 1 else 0
        r = int(bits[pos : pos + r_bits], 2) if r_bits > 0 else 0
        pos += r_bits
        if r >= cutoff:
            r = (r << 1) | int(bits[pos], 2)
            pos += 1
            r -= cutoff
        return q * m + r, pos


class DyNEDcCompressorV2:
    """DyNEDc v2: Adaptive position-based Golomb coding for binary spike trains.

    Instead of RLE + Huffman (v1):
      1. Pick the minority symbol (whichever of 0s or 1s is fewer)
      2. Extract positions of minority symbols, delta-encode to gaps
      3. Golomb-code each gap (optimal for geometric distributions)

    Advantages over v1:
      - Works well for both sparse AND dense spike trains (adaptive)
      - No codebook overhead (just parameter M, auto-computed from density)
      - Provably optimal for memoryless spike processes
      - Self-contained compressed stream (57-bit header)
    """

    def compress(self, binary_data: str) -> tuple[str, dict]:
        """Compress binary spike train using adaptive position-based Golomb coding."""
        if not binary_data:
            return "", {"original_length": 0, "compressed_length": 0}

        n = len(binary_data)
        ones = [i for i, c in enumerate(binary_data) if c == "1"]
        k1 = len(ones)
        k0 = n - k1

        # Encode whichever set of positions is smaller
        if k1 <= k0:
            positions, k, invert = ones, k1, 0
        else:
            positions = [i for i, c in enumerate(binary_data) if c == "0"]
            k, invert = k0, 1

        if k == 0:
            # All same symbol  -  header only
            header = str(invert) + format(n, "020b") + format(0, "020b") + format(1, "016b")
            return header, {
                "original_length": n,
                "num_positions": 0,
                "inverted": bool(invert),
                "golomb_m": 1,
                "compressed_length": len(header),
                "compression_ratio": len(header) / n,
            }

        # Delta-encode positions
        gaps = [positions[0]]
        for i in range(1, k):
            gaps.append(positions[i] - positions[i - 1])

        # Optimal Golomb parameter from minority density
        p = k / n
        m = GolombCoder.optimal_m(p)

        # Header: invert (1b) + length (20b) + count (20b) + M (16b) = 57 bits
        header = str(invert) + format(n, "020b") + format(k, "020b") + format(m, "016b")
        body = "".join(GolombCoder.encode(g, m) for g in gaps)
        compressed = header + body

        return compressed, {
            "original_length": n,
            "num_positions": k,
            "inverted": bool(invert),
            "golomb_m": m,
            "compressed_length": len(compressed),
            "compression_ratio": len(compressed) / n,
        }

    def decompress(self, compressed_data: str) -> str:
        """Decompress Golomb-coded spike train back to binary string."""
        if not compressed_data:
            return ""

        invert = int(compressed_data[0])
        n = int(compressed_data[1:21], 2)
        k = int(compressed_data[21:41], 2)
        m = int(compressed_data[41:57], 2)

        # Fill with majority symbol, mark minority positions
        fill = "1" if invert else "0"
        mark = "0" if invert else "1"

        if k == 0:
            return fill * n

        # Decode Golomb-coded gaps
        pos = 57
        gaps = []
        for _ in range(k):
            gap, pos = GolombCoder.decode(compressed_data, pos, m)
            gaps.append(gap)

        # Delta-decode back to absolute positions
        positions = [gaps[0]]
        for i in range(1, len(gaps)):
            positions.append(positions[-1] + gaps[i])

        # Reconstruct binary string
        result = [fill] * n
        for idx in positions:
            result[idx] = mark

        return "".join(result)


class DyNEDcCompressorV3:
    """DyNEDc v3: Full hybrid compression picking the best strategy.

    Tries four modes and picks whichever gives the smallest output:
      1. Variable-length RLE (3-bit count_size + variable count + symbol per run)
      2. Chunk-based Huffman (groups bits into fixed-size chunks before Huffman)
      3. RLE + chunk Huffman (hybrid)
      4. Alternating-run RLE + Huffman on run-length values (same as V1/DyNEDc)

    Note: Uses a safe RLE that splits runs > 255 to stay within the 3-bit
    count_size field (max 8-bit count = 255).
    """

    _MAX_RUN = 255  # max representable with 3-bit count_size (8-bit count)

    def __init__(self, chunk_size=4):
        self.chunk_size = chunk_size
        self._mode = None
        self._huff_codes = {}
        self._alt_start = "0"

    @staticmethod
    def _flush_run(compressed, count, value):
        """Encode one RLE entry: 3-bit count_size + variable count + 1-bit value."""
        count_binary = format(count, "b")
        count_size = format(len(count_binary) - 1, "03b")
        compressed.append(count_size + count_binary + value)

    def _rle_compress(self, binary_data):
        """Variable-length RLE, safe for arbitrary run lengths."""
        if not binary_data:
            return ""
        compressed = []
        count = 1
        current = binary_data[0]

        for i in range(1, len(binary_data)):
            if binary_data[i] == current:
                count += 1
                if count == self._MAX_RUN:
                    self._flush_run(compressed, count, current)
                    count = 0
            else:
                if count > 0:
                    self._flush_run(compressed, count, current)
                current = binary_data[i]
                count = 1

        if count > 0:
            self._flush_run(compressed, count, current)

        return "".join(compressed)

    @staticmethod
    def _rle_decompress(compressed_data):
        """Decompress variable-length RLE (same format as RunLengthEncoding)."""
        if not compressed_data:
            return ""
        decompressed = []
        i = 0
        while i < len(compressed_data):
            count_size = int(compressed_data[i : i + 3], 2) + 1
            i += 3
            count = int(compressed_data[i : i + count_size], 2)
            i += count_size
            value = compressed_data[i]
            i += 1
            decompressed.append(value * count)
        return "".join(decompressed)

    @staticmethod
    def _alt_rle_compress(binary_data):
        """Alternating-run RLE: extract run lengths + starting value."""
        runs = []
        current = binary_data[0]
        count = 1
        for c in binary_data[1:]:
            if c == current:
                count += 1
            else:
                runs.append(count)
                current = c
                count = 1
        runs.append(count)
        return runs, binary_data[0]

    def _alt_huffman_compress(self, binary_data):
        """V1-style: alternating runs + Huffman on run-length values."""
        runs, start = self._alt_rle_compress(binary_data)
        symbols = [str(r) for r in runs]

        freq = Counter(symbols)
        if len(freq) == 1:
            codes = {list(freq.keys())[0]: "0"}
        else:
            heap = [HuffmanNode(s, c) for s, c in freq.items()]
            heapq.heapify(heap)
            while len(heap) > 1:
                left = heapq.heappop(heap)
                right = heapq.heappop(heap)
                internal = HuffmanNode(None, left.freq + right.freq)
                internal.left = left
                internal.right = right
                heapq.heappush(heap, internal)
            codes = {}
            self._build_codes(heap[0], "", codes)

        compressed = "".join(codes[s] for s in symbols)
        return compressed, codes, start

    @staticmethod
    def _build_codes(node, prefix, codes):
        if node is None:
            return
        if node.char is not None:
            codes[node.char] = prefix or "0"
            return
        DyNEDcCompressorV3._build_codes(node.left, prefix + "0", codes)
        DyNEDcCompressorV3._build_codes(node.right, prefix + "1", codes)

    def compress(self, binary_data: str) -> tuple[str, dict]:
        """Compress using hybrid mode selection  -  picks the best ratio."""
        if not binary_data:
            return "", {"original_length": 0, "compressed_length": 0}

        n = len(binary_data)

        # Mode 1: Variable-length RLE only
        rle_compressed = self._rle_compress(binary_data)

        # Mode 2: Chunk-based Huffman only
        huff_only = HuffmanCompressor(chunk_size=self.chunk_size)
        huff_compressed, _ = huff_only.compress(binary_data)

        # Mode 3: RLE + chunk Huffman (hybrid)
        huff_hybrid = HuffmanCompressor(chunk_size=self.chunk_size)
        hybrid_compressed, _ = huff_hybrid.compress(rle_compressed)

        # Mode 4: Alternating-run RLE + Huffman on values (V1-style)
        alt_compressed, alt_codes, alt_start = self._alt_huffman_compress(binary_data)

        # Pick the best mode
        candidates = {
            "rle": (rle_compressed, len(rle_compressed), {}, ""),
            "huffman": (huff_compressed, len(huff_compressed), dict(huff_only.codes), ""),
            "hybrid": (hybrid_compressed, len(hybrid_compressed), dict(huff_hybrid.codes), ""),
            "alt_huffman": (alt_compressed, len(alt_compressed), alt_codes, alt_start),
        }

        best_mode = min(candidates, key=lambda m: candidates[m][1])
        best_compressed, best_len, best_codes, best_start = candidates[best_mode]

        self._mode = best_mode
        self._huff_codes = best_codes
        self._alt_start = best_start

        return best_compressed, {
            "original_length": n,
            "compressed_length": best_len,
            "compression_ratio": best_len / n,
            "mode": best_mode,
        }

    def decompress(self, compressed_data: str) -> str:
        """Decompress based on the mode used during compression."""
        if not compressed_data:
            return ""

        if self._mode == "rle":
            return self._rle_decompress(compressed_data)

        if self._mode == "alt_huffman":
            # Decode Huffman to get run lengths, then reconstruct
            reverse = {v: k for k, v in self._huff_codes.items()}
            runs = []
            buffer = ""
            for bit in compressed_data:
                buffer += bit
                if buffer in reverse:
                    runs.append(int(reverse[buffer]))
                    buffer = ""
            result = []
            current = self._alt_start
            for length in runs:
                result.append(current * length)
                current = "1" if current == "0" else "0"
            return "".join(result)

        # Decode chunk-based Huffman
        reverse = {v: k for k, v in self._huff_codes.items()}
        decoded = []
        buffer = ""
        for bit in compressed_data:
            buffer += bit
            if buffer in reverse:
                decoded.append(reverse[buffer])
                buffer = ""
        huff_output = "".join(decoded)

        if self._mode == "huffman":
            return huff_output
        # hybrid: Huffman output is the RLE-compressed data
        return self._rle_decompress(huff_output)


# =============================================================================
# DyNEDc V4  -  Numba-Accelerated Hybrid Compression
# =============================================================================

if HAS_NUMBA:

    @njit(cache=True)
    def _extract_runs_jit(arr):
        """Extract alternating run lengths from uint8 array (~100x faster)."""
        n = len(arr)
        if n == 0:
            return np.zeros(0, dtype=np.int64), np.int8(0)
        # Worst case: every element is a new run
        runs = np.empty(n, dtype=np.int64)
        num_runs = 0
        current = arr[0]
        count = 1
        for i in range(1, n):
            if arr[i] == current:
                count += 1
            else:
                runs[num_runs] = count
                num_runs += 1
                current = arr[i]
                count = 1
        runs[num_runs] = count
        num_runs += 1
        return runs[:num_runs], np.int8(arr[0])

    @njit(cache=True)
    def _rle_size_jit(runs, max_run):
        """Compute variable-length RLE compressed size in bits."""
        size = np.int64(0)
        for i in range(len(runs)):
            count = runs[i]
            while count > max_run:
                # 3-bit count_size + bits_in(max_run) + 1-bit value
                cb_len = 0
                tmp = max_run
                while tmp > 0:
                    cb_len += 1
                    tmp >>= 1
                size += 3 + cb_len + 1
                count -= max_run
            if count > 0:
                cb_len = 0
                tmp = count
                while tmp > 0:
                    cb_len += 1
                    tmp >>= 1
                size += 3 + cb_len + 1
        return size

    @njit(cache=True)
    def _pack_chunks_jit(arr, chunk_size):
        """Pack bit array into integer chunk keys."""
        n = len(arr)
        n_full = n - n % chunk_size
        n_chunks = n_full // chunk_size + (1 if n % chunk_size > 0 else 0)
        result = np.empty(n_chunks, dtype=np.int64)
        idx = 0
        for i in range(0, n_full, chunk_size):
            val = np.int64(0)
            for j in range(chunk_size):
                val = val * 2 + np.int64(arr[i + j])
            result[idx] = val
            idx += 1
        if n % chunk_size > 0:
            val = np.int64(0)
            for j in range(n_full, n):
                val = val * 2 + np.int64(arr[j])
            result[idx] = val
            idx += 1
        return result[:idx]


class DyNEDcCompressorV4:
    """DyNEDc v4: Numba-accelerated hybrid compression for spike trains.

    Same four-mode strategy as V3 but accepts tensors/arrays directly
    and uses numba JIT-compiled loops for RLE extraction, size
    computation, and chunk packing (~100x faster than Python iteration).

    Modes:
      1. Variable-length RLE
      2. Chunk-based Huffman
      3. RLE + chunk Huffman (hybrid)
      4. Alternating-run RLE + Huffman on run-length values
    """

    _MAX_RUN = 255

    def __init__(self, chunk_size=4):
        self.chunk_size = chunk_size
        self._mode = None
        self._huff_codes = {}
        self._alt_start = "0"

    # -- Primitives ------------------------------------------------------

    @staticmethod
    def _to_array(spike_data):
        """Convert any input to a flat uint8 numpy array."""
        if isinstance(spike_data, str):
            return np.frombuffer(spike_data.encode("ascii"), dtype=np.uint8) - ord("0")
        if hasattr(spike_data, "numpy"):
            return spike_data.detach().cpu().flatten().numpy().astype(np.uint8)
        return np.asarray(spike_data, dtype=np.uint8).flatten()

    @staticmethod
    def _extract_runs(arr):
        """Extract alternating run lengths. Uses numba if available."""
        if len(arr) == 0:
            return np.zeros(0, dtype=np.int64), "0"
        if HAS_NUMBA:
            runs, start_val = _extract_runs_jit(arr)
            return runs, str(int(start_val))
        # Fallback: numpy vectorised
        changes = np.diff(arr) != 0
        change_idx = np.flatnonzero(changes) + 1
        boundaries = np.concatenate(([0], change_idx, [len(arr)]))
        return np.diff(boundaries), str(int(arr[0]))

    @staticmethod
    def _arr_to_chunk_ints(arr, chunk_size):
        """Pack bit array into integer chunk keys. Uses numba if available."""
        if HAS_NUMBA:
            return _pack_chunks_jit(arr, chunk_size)
        # Fallback: numpy
        n = len(arr) - len(arr) % chunk_size
        powers = (2 ** np.arange(chunk_size - 1, -1, -1)).astype(np.int64)
        if n > 0:
            chunk_ints = (arr[:n].reshape(-1, chunk_size).astype(np.int64) * powers).sum(axis=1)
        else:
            chunk_ints = np.array([], dtype=np.int64)
        if len(arr) - n > 0:
            tail = arr[n:]
            tp = (2 ** np.arange(len(tail) - 1, -1, -1)).astype(np.int64)
            chunk_ints = np.append(chunk_ints, int((tail.astype(np.int64) * tp).sum()))
        return chunk_ints

    # -- RLE encoding from pre-extracted runs ----------------------------

    @classmethod
    def _rle_from_runs(cls, runs, start_value):
        """Encode pre-extracted runs to variable-length RLE binary format."""
        compressed = []
        current = start_value
        for count in runs:
            while count > cls._MAX_RUN:
                cb = format(cls._MAX_RUN, "b")
                compressed.append(format(len(cb) - 1, "03b") + cb + current)
                count -= cls._MAX_RUN
            if count > 0:
                cb = format(count, "b")
                compressed.append(format(len(cb) - 1, "03b") + cb + current)
            current = "1" if current == "0" else "0"
        return "".join(compressed)

    @staticmethod
    def _rle_size_from_runs(runs):
        """Compute compressed size of variable-length RLE without building the string."""
        if HAS_NUMBA:
            return int(_rle_size_jit(runs, np.int64(255)))
        size = 0
        for count in runs:
            while count > 255:
                size += 3 + len(format(255, "b")) + 1
                count -= 255
            if count > 0:
                size += 3 + len(format(count, "b")) + 1
        return size

    # -- Huffman helpers -------------------------------------------------

    @staticmethod
    def _build_huffman_codes(freq):
        """Build Huffman codes from a frequency dict."""
        if len(freq) == 1:
            return {list(freq.keys())[0]: "0"}
        heap = [HuffmanNode(sym, c) for sym, c in freq.items()]
        heapq.heapify(heap)
        while len(heap) > 1:
            left = heapq.heappop(heap)
            right = heapq.heappop(heap)
            internal = HuffmanNode(None, left.freq + right.freq)
            internal.left = left
            internal.right = right
            heapq.heappush(heap, internal)
        codes = {}
        DyNEDcCompressorV4._walk(heap[0], "", codes)
        return codes

    @staticmethod
    def _walk(node, prefix, codes):
        if node is None:
            return
        if node.char is not None:
            codes[node.char] = prefix or "0"
            return
        DyNEDcCompressorV4._walk(node.left, prefix + "0", codes)
        DyNEDcCompressorV4._walk(node.right, prefix + "1", codes)

    @staticmethod
    def _huffman_size(freq, codes):
        """Compute Huffman compressed size from frequency dict and codes."""
        return sum(count * len(codes[sym]) for sym, count in freq.items())

    # -- Main compress interface -----------------------------------------

    def compress(self, spike_data) -> tuple[str, dict]:
        """Compress a spike train using hybrid mode selection.

        Args:
            spike_data: tensor, numpy array, or binary string of 0s and 1s.

        Returns:
            (compressed_str, info_dict) same format as V3.
        """
        arr = self._to_array(spike_data)
        n = len(arr)
        self._original_length = n  # needed for correct decompression
        if n == 0:
            return "", {"original_length": 0, "compressed_length": 0}

        # Extract runs once (numpy vectorised)  -  reused by modes 1, 3, 4
        runs, start = self._extract_runs(arr)

        # -- Compute sizes for all modes (without building full strings
        #    for modes we won't use) -------------------------------------

        # Mode 1: Variable-length RLE  -  compute size from runs
        rle_size = self._rle_size_from_runs(runs)

        # Mode 2: Chunk-based Huffman  -  compute size from chunk frequencies
        chunk_ints = self._arr_to_chunk_ints(arr, self.chunk_size)
        chunk_freq = Counter(chunk_ints.tolist())
        huff_codes = self._build_huffman_codes(chunk_freq)
        huff_size = self._huffman_size(chunk_freq, huff_codes)

        # Mode 4: Alternating-run RLE + Huffman on run-length values
        runs_list = runs.tolist() if isinstance(runs, np.ndarray) else runs
        sym_freq = Counter(runs_list)
        alt_codes = self._build_huffman_codes(sym_freq)
        alt_size = self._huffman_size(sym_freq, alt_codes)

        # Mode 3: RLE + chunk Huffman  -  only compute if potentially best
        # (skip the expensive RLE string -> re-chunk -> Huffman pipeline
        #  if another mode is already clearly better)
        min_so_far = min(rle_size, huff_size, alt_size)
        hybrid_size = min_so_far + 1  # assume worse unless computed
        hybrid_codes = {}
        rle_compressed = None

        if min_so_far > n * 0.05:  # only try hybrid if no mode is already <5%
            rle_compressed = self._rle_from_runs(runs, start)
            rle_arr = np.frombuffer(rle_compressed.encode("ascii"), dtype=np.uint8) - ord("0")
            rle_chunks = self._arr_to_chunk_ints(rle_arr, self.chunk_size)
            rle_chunk_freq = Counter(rle_chunks.tolist())
            hybrid_codes = self._build_huffman_codes(rle_chunk_freq)
            hybrid_size = self._huffman_size(rle_chunk_freq, hybrid_codes)

        # Pick the best mode
        sizes = {
            "rle": rle_size,
            "huffman": huff_size,
            "hybrid": hybrid_size,
            "alt_huffman": alt_size,
        }
        best_mode = min(sizes, key=sizes.get)
        best_size = sizes[best_mode]

        # Only build the compressed string for the winning mode
        if best_mode == "rle":
            if rle_compressed is None:
                rle_compressed = self._rle_from_runs(runs, start)
            best_compressed = rle_compressed
            best_codes = {}
            best_start = ""
        elif best_mode == "huffman":
            chunk_list = chunk_ints.tolist()
            best_compressed = "".join(huff_codes[k] for k in chunk_list)
            best_codes = huff_codes
            best_start = ""
        elif best_mode == "hybrid":
            if rle_compressed is None:
                rle_compressed = self._rle_from_runs(runs, start)
            rle_arr = np.frombuffer(rle_compressed.encode("ascii"), dtype=np.uint8) - ord("0")
            rle_chunks = self._arr_to_chunk_ints(rle_arr, self.chunk_size)
            best_compressed = "".join(hybrid_codes[k] for k in rle_chunks.tolist())
            best_codes = hybrid_codes
            best_start = ""
        else:  # alt_huffman
            best_compressed = "".join(alt_codes[r] for r in runs_list)
            best_codes = alt_codes
            best_start = start

        self._mode = best_mode
        self._huff_codes = best_codes
        self._alt_start = best_start

        return best_compressed, {
            "original_length": n,
            "compressed_length": best_size,
            "compression_ratio": best_size / n,
            "mode": best_mode,
        }

    def decompress(self, compressed_data: str) -> str:
        """Decompress based on the mode used during compression."""
        if not compressed_data:
            return ""

        if self._mode == "rle":
            return DyNEDcCompressorV3._rle_decompress(compressed_data)

        if self._mode == "alt_huffman":
            reverse = {v: k for k, v in self._huff_codes.items()}
            runs = []
            buffer = ""
            for bit in compressed_data:
                buffer += bit
                if buffer in reverse:
                    runs.append(int(reverse[buffer]))
                    buffer = ""
            result = []
            current = self._alt_start
            for length in runs:
                result.append(current * length)
                current = "1" if current == "0" else "0"
            return "".join(result)

        # Huffman or hybrid  -  decode chunk-based Huffman
        reverse = {v: k for k, v in self._huff_codes.items()}
        decoded = []
        buffer = ""
        for bit in compressed_data:
            buffer += bit
            if buffer in reverse:
                decoded.append(reverse[buffer])
                buffer = ""

        # Convert chunk integers back to zero-padded binary strings.
        cs = self.chunk_size
        huff_output = "".join(format(int(k), f"0{cs}b") for k in decoded)

        # The last chunk may have been shorter than chunk_size during
        # compression, so truncate to the original length.
        orig_len = getattr(self, "_original_length", None)
        if orig_len is not None and len(huff_output) > orig_len:
            huff_output = huff_output[:orig_len]

        if self._mode == "huffman":
            return huff_output
        return DyNEDcCompressorV3._rle_decompress(huff_output)


# =============================================================================
# Visualization Utilities
# =============================================================================


def _spike_train_to_binary(step_signal):
    """Convert boolean step signal to binary string for compression."""
    return "".join(str(int(s)) for s in step_signal)


def _print_metrics(name, original_size, compressed_size):
    """Print compression statistics and return formatted ratio/saving strings."""
    ratio = compressed_size / original_size if original_size > 0 else 0
    saving = (1 - ratio) * 100 if original_size > 0 else 0
    print(f"{name}:")
    print(f"    Original data size: {original_size} bits")
    print(f"    Compressed data size: {compressed_size} bits")
    print(f"    Compression ratio: {ratio:.3f}")
    print(f"    Space saving: {saving:.1f}%")
    return f"{ratio:.3f}", f"{saving:.1f}%"


def _make_spike_trace(freqs, step_signal):
    """Create plotly spike trace showing spikes at unit height."""
    spike_indices = np.where(step_signal == 1)[0]
    spike_freqs = freqs[spike_indices]
    x_values = np.repeat(spike_freqs, 3)
    y_values = np.tile([0, 1, None], len(spike_freqs))
    return go.Scatter(
        x=x_values,
        y=y_values,
        mode="lines",
        line=dict(color="blue", width=1),
        hoverinfo="skip",
    )


def create_spike_plot(freqs, step_signal):
    """Create a plotly spike plot for the original (uncompressed) spike train."""
    spike_indices = np.where(step_signal == 1)[0]
    print(f"Number of spikes: {len(spike_indices)}")
    return _make_spike_trace(freqs, step_signal)


def create_rle_spike_plot(freqs, step_signal):
    """Compress spike train with RLE and create visualization."""
    spike_binary = _spike_train_to_binary(step_signal)
    n_spikes = int(np.sum(step_signal))
    print(f"Number of spikes: {n_spikes}")

    compressor = RunLengthEncoding()
    compressed, _ = compressor.compress(spike_binary)

    # Verify lossless round-trip
    assert compressor.decompress(compressed) == spike_binary, "RLE round-trip failed!"

    cr, sa = _print_metrics("Run-Length Encoding", len(spike_binary), len(compressed))
    return _make_spike_trace(freqs, step_signal), cr, sa


def create_lzw_spike_plot(freqs, step_signal):
    """Compress spike train with LZW and create visualization."""
    spike_binary = _spike_train_to_binary(step_signal)
    n_spikes = int(np.sum(step_signal))
    print(f"Number of spikes: {n_spikes}")

    compressor = LZWCompressor()
    compressed = compressor.compress(spike_binary)

    compressed_size = int(len(compressed) * np.ceil(np.log2(max(compressed) + 1)))

    # Verify lossless round-trip
    assert compressor.decompress(compressed) == spike_binary, "LZW round-trip failed!"

    cr, sa = _print_metrics("LZW Compression", len(spike_binary), compressed_size)
    return _make_spike_trace(freqs, step_signal), cr, sa


def create_huffman_spike_plot(freqs, step_signal):
    """Compress spike train with Huffman coding and create visualization."""
    spike_binary = _spike_train_to_binary(step_signal)
    n_spikes = int(np.sum(step_signal))
    print(f"Number of spikes: {n_spikes}")

    compressor = HuffmanCompressor()
    compressed, info = compressor.compress(spike_binary)

    # Verify lossless round-trip
    assert compressor.decompress(compressed, info["codes"]) == spike_binary, "Huffman round-trip failed!"

    cr, sa = _print_metrics("Huffman Compression", len(spike_binary), len(compressed))
    return _make_spike_trace(freqs, step_signal), cr, sa


def create_delta_spike_plot(freqs, step_signal):
    """Compress spike positions with delta encoding and create visualization."""
    spike_indices = np.where(step_signal == 1)[0]
    n_spikes = len(spike_indices)
    print(f"Number of spikes: {n_spikes}")

    compressor = DeltaCompressor()
    compressed = compressor.compress(spike_indices.tolist())

    compressed_size = sum(len(bin(abs(x))[2:]) + 1 for x in compressed)
    original_size = len(step_signal)

    # Verify lossless round-trip
    assert compressor.decompress(compressed) == spike_indices.tolist(), "Delta round-trip failed!"

    cr, sa = _print_metrics("Delta Compression", original_size, compressed_size)
    return _make_spike_trace(freqs, step_signal), cr, sa


def create_bwt_spike_plot(freqs, step_signal):
    """Apply BWT to spike train and report metrics."""
    spike_binary = _spike_train_to_binary(step_signal)
    n_spikes = int(np.sum(step_signal))
    print(f"Number of spikes: {n_spikes}")

    bwt = BWTransform()
    transformed, index = bwt.transform(spike_binary)

    compressed_size = len(transformed) + len(bin(index)[2:])

    cr, sa = _print_metrics("Burrows-Wheeler Transform", len(spike_binary), compressed_size)
    return _make_spike_trace(freqs, step_signal), cr, sa


def create_lz77_spike_plot(freqs, step_signal):
    """Compress spike train with LZ77 and create visualization."""
    spike_binary = _spike_train_to_binary(step_signal)
    n_spikes = int(np.sum(step_signal))
    print(f"Number of spikes: {n_spikes}")

    compressor = LZ77Compressor(window_size=8)
    compressed = compressor.compress(spike_binary)

    compressed_size = sum(1 if isinstance(x, str) else (len(bin(x[0])[2:]) + len(bin(x[1])[2:])) for x in compressed)

    cr, sa = _print_metrics("LZ77", len(spike_binary), compressed_size)
    return _make_spike_trace(freqs, step_signal), cr, sa


def create_dynedc_spike_plot(freqs, step_signal):
    """Compress spike train with DyNEDc and create visualization."""
    spike_binary = _spike_train_to_binary(step_signal)
    n_spikes = int(np.sum(step_signal))
    print(f"Number of spikes: {n_spikes}")

    compressor = DyNEDcCompressor()
    compressed, _ = compressor.compress(spike_binary)

    # Verify lossless round-trip
    assert compressor.decompress(compressed) == spike_binary, "DyNEDc round-trip failed!"

    cr, sa = _print_metrics("DyNEDc Compression", len(spike_binary), len(compressed))
    return _make_spike_trace(freqs, step_signal), cr, sa


def create_dynedc_v2_spike_plot(freqs, step_signal):
    """Compress spike train with DyNEDc V2 (Golomb) and create visualization."""
    spike_binary = _spike_train_to_binary(step_signal)
    n_spikes = int(np.sum(step_signal))
    print(f"Number of spikes: {n_spikes}")

    compressor = DyNEDcCompressorV2()
    compressed, info = compressor.compress(spike_binary)

    # Verify lossless round-trip
    assert compressor.decompress(compressed) == spike_binary, "DyNEDc V2 round-trip failed!"

    print(f"    Golomb parameter M: {info.get('golomb_m', 'N/A')}")
    cr, sa = _print_metrics("DyNEDc V2 (Golomb)", len(spike_binary), len(compressed))
    return _make_spike_trace(freqs, step_signal), cr, sa


def create_dynedc_v3_spike_plot(freqs, step_signal):
    """Compress spike train with DyNEDc V3 (hybrid) and create visualization."""
    spike_binary = _spike_train_to_binary(step_signal)
    n_spikes = int(np.sum(step_signal))
    print(f"Number of spikes: {n_spikes}")

    compressor = DyNEDcCompressorV3()
    compressed, info = compressor.compress(spike_binary)

    # Verify lossless round-trip
    assert compressor.decompress(compressed) == spike_binary, "DyNEDc V3 round-trip failed!"

    print(f"    Mode selected: {info.get('mode', 'N/A')}")
    cr, sa = _print_metrics("DyNEDc V3 (Hybrid)", len(spike_binary), len(compressed))
    return _make_spike_trace(freqs, step_signal), cr, sa


def create_dynedc_v4_spike_plot(freqs, step_signal):
    """Compress spike train with DyNEDc V4 (numba-accelerated hybrid) and create visualization."""
    spike_binary = _spike_train_to_binary(step_signal)
    n_spikes = int(np.sum(step_signal))
    print(f"Number of spikes: {n_spikes}")

    compressor = DyNEDcCompressorV4()
    compressed, info = compressor.compress(spike_binary)

    print(f"    Mode selected: {info.get('mode', 'N/A')}")
    cr, sa = _print_metrics("DyNEDc Compression", len(spike_binary), len(compressed))
    return _make_spike_trace(freqs, step_signal), cr, sa


# =============================================================================
# Benchmarking - Compression Summary
# =============================================================================


def compression_summary(step_signal):
    """
    Run all compression methods on a spike train and return a summary.

    Args:
        step_signal: Boolean array from generate_step_signal().

    Returns:
        dict with original_bits, num_spikes, and per-method metrics including
        compressed_bits, bits_per_spike, compression_ratio, space_saving_pct.
    """
    spike_binary = _spike_train_to_binary(step_signal)
    spike_indices = np.where(step_signal == 1)[0]
    n_spikes = len(spike_indices)
    original_bits = len(spike_binary)

    methods = {}

    # RLE
    rle = RunLengthEncoding()
    rle_compressed, _ = rle.compress(spike_binary)
    methods["RLE"] = len(rle_compressed)

    # LZW
    lzw = LZWCompressor()
    lzw_compressed = lzw.compress(spike_binary)
    methods["LZW"] = int(len(lzw_compressed) * np.ceil(np.log2(max(lzw_compressed) + 1))) if lzw_compressed else 0

    # Huffman
    huffman = HuffmanCompressor()
    huffman_compressed, _ = huffman.compress(spike_binary)
    methods["Huffman"] = len(huffman_compressed)

    # Delta (on spike positions)
    delta = DeltaCompressor()
    delta_compressed = delta.compress(spike_indices.tolist())
    methods["Delta"] = sum(len(bin(abs(x))[2:]) + 1 for x in delta_compressed) if delta_compressed else 0

    # BWT (skip if signal too long)
    if len(spike_binary) <= BWTransform.MAX_LENGTH:
        bwt = BWTransform()
        bwt_transformed, bwt_index = bwt.transform(spike_binary)
        methods["BWT"] = len(bwt_transformed) + len(bin(bwt_index)[2:])
    else:
        methods["BWT"] = None  # skipped

    # LZ77
    lz77 = LZ77Compressor()
    lz77_compressed = lz77.compress(spike_binary)
    methods["LZ77"] = sum(
        1 if isinstance(x, str) else (len(bin(x[0])[2:]) + len(bin(x[1])[2:])) for x in lz77_compressed
    )

    # DyNEDc
    dynedc = DyNEDcCompressor()
    dynedc_compressed, _ = dynedc.compress(spike_binary)
    methods["DyNEDc"] = len(dynedc_compressed)

    # DyNEDc V2 (Golomb)
    dynedc_v2 = DyNEDcCompressorV2()
    dynedc_v2_compressed, _ = dynedc_v2.compress(spike_binary)
    methods["DyNEDc V2"] = len(dynedc_v2_compressed)

    # DyNEDc V3 (Hybrid: variable-length RLE + chunk Huffman + mode selection)
    dynedc_v3 = DyNEDcCompressorV3()
    dynedc_v3_compressed, _ = dynedc_v3.compress(spike_binary)
    methods["DyNEDc V3"] = len(dynedc_v3_compressed)

    # DyNEDc V4 (Numba-accelerated hybrid)
    dynedc_v4 = DyNEDcCompressorV4()
    dynedc_v4_compressed, _ = dynedc_v4.compress(spike_binary)
    methods["DyNEDc V4"] = len(dynedc_v4_compressed)

    results = {}
    for name, bits in methods.items():
        if bits is None:
            results[name] = {
                "compressed_bits": None,
                "bits_per_spike": None,
                "compression_ratio": None,
                "space_saving_pct": None,
            }
            continue
        ratio = bits / original_bits if original_bits > 0 else 0
        bps = bits / n_spikes if n_spikes > 0 else 0
        results[name] = {
            "compressed_bits": bits,
            "bits_per_spike": round(bps, 2),
            "compression_ratio": round(ratio, 4),
            "space_saving_pct": round((1 - ratio) * 100, 1),
        }

    return {"original_bits": original_bits, "num_spikes": n_spikes, "methods": results}


def print_compression_table(summary):
    """Print a formatted comparison table from compression_summary output."""
    print(f"\nSpike Train: {summary['original_bits']} bits, {summary['num_spikes']} spikes")
    print(f"{'Method':<12} {'Compressed':>12} {'Ratio':>8} {'Bits/Spike':>12} {'Saving':>8}")
    print("-" * 56)
    for name, data in summary["methods"].items():
        if data["compressed_bits"] is None:
            print(f"{name:<12} {'(skipped)':>12} {'N/A':>8} {'N/A':>12} {'N/A':>8}")
            continue
        print(
            f"{name:<12} {data['compressed_bits']:>12} "
            f"{data['compression_ratio']:>8.4f} "
            f"{data['bits_per_spike']:>12.2f} "
            f"{data['space_saving_pct']:>7.1f}%"
        )


def batch_compression_stats(step_signals):
    """
    Compute compression statistics across multiple samples.

    Args:
        step_signals: List of boolean step signal arrays.

    Returns:
        dict with mean/std/min/max compression ratio per method.
    """
    all_results = defaultdict(list)

    for step_signal in step_signals:
        summary = compression_summary(step_signal)
        for method, data in summary["methods"].items():
            if data["compression_ratio"] is not None:
                all_results[method].append(data["compression_ratio"])

    stats = {}
    for method, ratios in all_results.items():
        ratios = np.array(ratios)
        stats[method] = {
            "mean_ratio": round(float(ratios.mean()), 4),
            "std_ratio": round(float(ratios.std()), 4),
            "min_ratio": round(float(ratios.min()), 4),
            "max_ratio": round(float(ratios.max()), 4),
            "n_samples": len(ratios),
        }

    return stats