How to Compress Time-Series Metrics in Python With Gorilla Delta-of-Delta and XOR Encoding

Every monitoring system stores the same two numbers over and over: a time and a measurement. Stored plainly, each pair takes 16 bytes, and a busy system records millions of them per second. In this tutorial you will build, in plain Python with only the standard library, the compression scheme from Facebook’s 2015 Gorilla paper, which Prometheus adapted for the chunks that store its metrics. You will test it on a real week of river measurements from the U.S. Geological Survey (USGS). The 32,272 bytes of one series shrink to 1,678 bytes (19.2 times smaller) and decode back to exactly the same numbers.

You will also see where the scheme falls short, which is just as useful. The same code compresses a CPU-usage gauge only 2.3 times, and plain xz compression does about twice as well on it. A one-line change then makes the gauge 8.6 times smaller. Understanding why that happens will help you more than memorizing the algorithm.

Why monitoring data is worth compressing

A time series is a list of samples, and a sample is one timestamp plus one value. For example, at 04:00 UTC on October 14, 2026 the Potomac River was carrying 15,400 cubic feet of water per second. Monitoring tools collect many series like this every few seconds. The obvious storage layout, an 8-byte integer for the time and an 8-byte floating-point number for the value, costs 16 bytes per sample.

The Gorilla paper (Pelkonen et al., VLDB 2015) puts the scale in perspective: “At 16 bytes per point, the resulting 16TB of RAM would be too resource intensive for practical deployment.” Its authors, working at Facebook, repurposed an existing scheme to “compress time series to an average of 1.37 bytes per point, a 12x reduction in size.” Prometheus reports a similar outcome for its own storage: “Prometheus stores an average of only 1-2 bytes per sample.”

Two properties of real measurements make this possible. Samples usually arrive at a steady interval, so the time between neighbors barely changes. And a value is usually the same as, or close to, the one before it. Gorilla stores only what changed, in as few bits as possible. You will build the two halves separately, prove that each one works, and then measure the whole thing.

Prerequisites

• Python 3.13 or newer. I ran everything on Python 3.13.14 on Windows 11 and did not try other versions or operating systems; the code is plain Python and should run the same elsewhere. Where I write python, use python3 if that is your command name.
• pytest 9.1.1 for the tests (pip install pytest==9.1.1). Nothing else is needed: the compressor itself uses only the standard library.
• Internet access once, in Step 1, to download about 320 KB of public data from USGS.
• Basic Python: functions, lists, tuples, integers and floats. You do not need to know bit operations; each one is explained where it first appears.
• About an hour to read and run everything.

Create an empty folder called gorilla_lab and work inside it. Every file in this tutorial goes in that folder, and every command assumes you run it from there.

How Gorilla works, in plain words

Five words appear all the time, so here they are once:

• Bit: a 0 or a 1. A byte is 8 bits.
• Delta: the difference between a number and the one before it.
• Delta of delta: the difference between neighboring deltas. For a steady clock it is zero.
• XOR (exclusive or): combines two bit patterns and gives 1 wherever they differ and 0 wherever they match. XOR of two identical patterns is all zeros.
• Chunk: a run of consecutive samples compressed together as one unit.

Idea 1: for timestamps, store the change in the change

If a scraper reads a gauge every 15 seconds, the timestamps in milliseconds go 0, 15000, 30000, 45000, and so on. The deltas are 15000, 15000, 15000, so the deltas of deltas are all 0. The paper’s own example has deltas of “60, 60, 59 and 61” seconds, which gives deltas of deltas of “0, -1 and 2”. Small numbers need few bits, so Gorilla writes each delta of delta with a short prefix that says how many bits follow:

The paper says the limits “were selected by sampling a set of real time series from the production system”, and it reports that “about 96% of all time stamps can be compressed to a single bit”.

Idea 2: for values, store only the bits that flipped

A Python float is a 64-bit IEEE 754 number: one sign bit, 11 exponent bits and 52 fraction bits. Two measurements that are close together usually share the sign, the exponent and the first few fraction bits. If you XOR the two 64-bit patterns, the shared bits turn into zeros. What remains is a pattern with a run of zeros on the left (leading zeros) and often another run on the right (trailing zeros). Only the bits in between, the meaningful bits, carry information. If the value did not change at all, the XOR is all zeros and costs a single bit.

The window idea matters because successive XORs tend to sit in the same place, so reusing the previous window skips the 11 bits of leading-zero and length fields. The paper observes that “Integer values compress especially well because the location of the one bits after the XOR operation is often the same for the whole time series, meaning most values have the same number of trailing zeros.” Keep that sentence in mind; Step 8 turns it into a practical trick.

The plan

1. Download a real week of sensor data and price the raw format.
2. Write a tool that writes and reads single bits.
3. Encode timestamps.
4. Encode values.
5. Combine both into encode() and decode().
6. Test it, then plant bugs to see whether the tests notice.
7. Measure it against general-purpose compressors.
8. Fix the weak spot with integers.
9. Choose the timestamp table that matches your clock.
10. Cut the stream into chunks.

Step 1: Download a real week of sensor data and price the raw format

Download the data

Synthetic numbers hide the messy parts, so start with real ones. USGS publishes readings from its stream gauges, and its copyright page says that “USGS-authored or produced data and information are considered to be in the U.S. Public Domain”, so they are safe to download and reuse. The script below asks the USGS Instantaneous Values service for site 01646500 (the Potomac River near Washington, D.C., at the Little Falls pump station) for seven days, from 00:00 on October 14, 2026 to 00:00 on October 21, 2026, US Eastern daylight time (UTC-4). Parameter code 00060 is streamflow in cubic feet per second, and 00065 is gage height, the water level in feet.

Save this as fetch_usgs.py:

"""fetch_usgs.py: download 7 days of 5-minute Potomac River readings from USGS and save them as a CSV."""
import csv
import datetime as dt
import hashlib
import json
import urllib.request
from pathlib import Path

URL = ("https://waterservices.usgs.gov/nwis/iv/?format=json&sites=01646500&parameterCd=00060,00065"
       "&startDT=2026-10-14T00:00:00-04:00&endDT=2026-10-21T00:00:00-04:00")
OUT = Path(__file__).parent / "data" / "potomac_7d.csv"

request = urllib.request.Request(URL, headers={"User-Agent": "gorilla-tutorial/1.0"})
with urllib.request.urlopen(request, timeout=60) as response:
    payload = json.load(response)

columns = {}   # parameter code -> {epoch_ms: reading as the text USGS sent}
for series in payload["value"]["timeSeries"]:
    code = series["variable"]["variableCode"][0]["value"]
    readings = {}
    for row in series["values"][0]["value"]:
        when = dt.datetime.fromisoformat(row["dateTime"])
        readings[int(when.timestamp() * 1000)] = row["value"]
    columns[code] = readings

times = sorted(set(columns["00060"]) & set(columns["00065"]))
OUT.parent.mkdir(exist_ok=True)
with open(OUT, "w", newline="") as f:
    writer = csv.writer(f, lineterminator="\n")
    writer.writerow(["epoch_ms", "flow_cfs", "stage_ft"])
    for ms in times:
        writer.writerow([ms, columns["00060"][ms], columns["00065"][ms]])

print(f"saved {len(times)} rows to {OUT.name}")
print("sha256:", hashlib.sha256(OUT.read_bytes()).hexdigest())

The service answers with JSON that holds one list of readings per parameter. Each reading carries an ISO date with a UTC offset and the value as text. The script converts each date to milliseconds since 1970-01-01 UTC, the unit Prometheus uses for timestamps, keeps only the moments present in both series, and writes a CSV that keeps USGS’s value text unchanged. Saving to a file means you download once and everything later runs offline on the same bytes. The file is written with lineterminator="\n" so that Windows does not add carriage returns; that makes the sha256 fingerprint the same on every operating system.

python fetch_usgs.py

saved 2017 rows to potomac_7d.csv
sha256: 9434d0fc6bdf1e44cfecbc07d43c9d0754e167bc2f2419f2bcdabfd27677fb37

Check it worked

You should get 2,017 rows: seven days of one reading every five minutes is 7 × 288 = 2,016 readings, plus the closing one at the end of the range. USGS labels this data provisional (its response carries the qualifier “Provisional data subject to revision.”). If USGS revised a reading after I downloaded the file, your sha256 will differ from mine and a few numbers below will differ slightly. The code and the lessons do not change.

Common mistake

Dropping the UTC offset and treating the local clock time as UTC. Differences between neighbors would still compress the same way, but your timestamps would be four hours off from real time.

Load it and price the raw format

Now a small module that loads the CSV, and a script that prices the raw format. Save these as datasets.py and step01_baseline.py:

"""datasets.py: the series used in this tutorial. The synthetic ones are added in Step 7."""
import csv
from pathlib import Path

DATA = Path(__file__).parent / "data"


def usgs_series():
    """Real readings: {"flow": [(epoch_ms, cubic feet per second)], "stage": [(epoch_ms, feet)]}."""
    flow, stage = [], []
    with open(DATA / "potomac_7d.csv", newline="") as f:
        for row in csv.DictReader(f):
            ms = int(row["epoch_ms"])
            flow.append((ms, float(row["flow_cfs"])))
            stage.append((ms, float(row["stage_ft"])))
    return {"flow": flow, "stage": stage}

"""step01_baseline.py: what two real sensor series cost before any compression."""
import hashlib
from pathlib import Path

from datasets import usgs_series

csv_file = Path(__file__).parent / "data" / "potomac_7d.csv"
print("csv sha256:", hashlib.sha256(csv_file.read_bytes()).hexdigest()[:16], "...")
for name, samples in usgs_series().items():
    print(f"{name}: {len(samples)} samples, first {samples[0]}, last {samples[-1]}")
    print(f"  raw size: {len(samples) * 16:,} bytes (8-byte timestamp + 8-byte float per sample)")

python step01_baseline.py

csv sha256: 9434d0fc6bdf1e44 ...
flow: 2017 samples, first (1791950400000, 15400.0), last (1792555200000, 5290.0)
  raw size: 32,272 bytes (8-byte timestamp + 8-byte float per sample)
stage: 2017 samples, first (1791950400000, 4.51), last (1792555200000, 3.4)
  raw size: 32,272 bytes (8-byte timestamp + 8-byte float per sample)

Each sample is an 8-byte timestamp plus an 8-byte float, so 2,017 samples cost 2,017 × 16 = 32,272 bytes per series. That is the number to beat. The flow starts at 15,400 cubic feet per second and falls to 5,290 over the week, while the gage height drops from 4.51 to 3.4 feet.

Step 2: Write and read single bits

Python can store whole bytes, but Gorilla needs fields that are 1, 2, 5, 6, 7, 9, 12 or 14 bits wide placed back to back. bits.py is a small tool for that. Save it:

"""bits.py: write and read fields that are not a whole number of bytes."""


class BitWriter:
    """Collects bits, most significant bit first, and packs them into bytes."""

    def __init__(self):
        self._buf = bytearray()
        self._acc = 0       # bits waiting to fill a byte (always fewer than 8)
        self._pending = 0   # how many bits are waiting in _acc
        self.nbits = 0      # total bits written so far

    def write_bits(self, value, width):
        if width < 0 or not 0 <= value < (1 << width):
            raise ValueError(f"{value} does not fit in {width} bits")
        self.nbits += width
        self._acc = (self._acc << width) | value
        self._pending += width
        while self._pending >= 8:
            self._pending -= 8
            self._buf.append((self._acc >> self._pending) & 0xFF)
        self._acc &= (1 << self._pending) - 1

    def getvalue(self):
        """All bits written so far, padded with zero bits up to a whole byte."""
        out = bytes(self._buf)
        if self._pending:
            out += bytes([(self._acc << (8 - self._pending)) & 0xFF])
        return out


class BitReader:
    """Reads fields back in the order they were written."""

    def __init__(self, data):
        self._data = bytes(data)
        self._pos = 0   # position in bits

    def read_bits(self, width):
        if width < 0 or self._pos + width > len(self._data) * 8:
            raise EOFError("ran out of bits")
        if width == 0:
            return 0
        first = self._pos >> 3
        last = (self._pos + width - 1) >> 3
        chunk = int.from_bytes(self._data[first:last + 1], "big")
        shift = (last - first + 1) * 8 - (self._pos & 7) - width
        self._pos += width
        return (chunk >> shift) & ((1 << width) - 1)

BitWriter keeps finished bytes in a bytearray and the leftover bits (always fewer than 8) in a small integer. write_bits(value, width) appends width bits, most significant bit first, and refuses a value that does not fit. The nbits counter is what lets us measure what every sample costs later. getvalue() pads the last partial byte with zero bits. BitReader does the reverse: it finds the bytes that hold the requested bits, turns them into an integer, and shifts and masks out the field. Asking for bits beyond the end raises EOFError.

Try it with step02_bits.py, which writes fields of 1, 2, 7 and 3 bits:

"""step02_bits.py: write fields of 1, 2, 7 and 3 bits, then read them back."""
from bits import BitReader, BitWriter

w = BitWriter()
for value, width in [(0b1, 1), (0b10, 2), (0b0111110, 7), (0b101, 3)]:
    w.write_bits(value, width)
data = w.getvalue()
print(w.nbits, "bits written ->", len(data), "bytes:", " ".join(f"{b:08b}" for b in data))
r = BitReader(data)
print("read back:", [r.read_bits(width) for width in (1, 2, 7, 3)])
try:
    w.write_bits(8, 3)
except ValueError as error:
    print("ValueError:", error)
try:
    BitReader(b"\xff").read_bits(9)
except EOFError as error:
    print("EOFError:", error)

python step02_bits.py

13 bits written -> 2 bytes: 11001111 10101000
read back: [1, 2, 62, 5]
ValueError: 8 does not fit in 3 bits
EOFError: ran out of bits

The four fields (1, 10, 0111110 and 101) add up to 13 bits. They were packed as 1100111110101 and padded with three zero bits to fill two bytes. The read-back list [1, 2, 62, 5] is exactly what went in (0111110 is 62 and 101 is 5). The two error lines show that both classes fail loudly rather than silently truncating, which the decoder will rely on in Step 6.

Check it worked

Your output should match the block above exactly: 13 bits in 2 bytes, the read-back list [1, 2, 62, 5] and the two error lines. If the list differs, recheck write_bits and read_bits in bits.py.

Common mistake

Forgetting that the padding looks like data. A reader that keeps reading until the bytes run out would decode up to seven phantom values from the padding bits. That is why our format stores the sample count first (Step 5).

Step 3: Encode timestamps as a delta of a delta

The timestamp half of the codec is a table of buckets and two small functions. Save this as gorilla.py; later steps append to the same file:

"""gorilla.py: Gorilla-style compression for one chunk of (timestamp, value) samples."""
import struct

from bits import BitReader, BitWriter

# Each bucket is (prefix bits, prefix length, payload width in bits).
PAPER_BUCKETS = ((0b10, 2, 7), (0b110, 3, 9), (0b1110, 4, 12), (0b1111, 4, 32))
PROMETHEUS_BUCKETS = ((0b10, 2, 14), (0b110, 3, 17), (0b1110, 4, 20), (0b1111, 4, 64))
MASK64 = (1 << 64) - 1


def fits(x, width):
    """True if x fits a payload of `width` bits. The paper's ranges run from
    -(2**(width-1) - 1) to 2**(width-1), so 7 bits hold -63 through 64."""
    return -((1 << (width - 1)) - 1) <= x <= (1 << (width - 1))


def write_dod(w, dod, buckets):
    """Write one delta of delta and return which bucket it used ("0" or the payload width)."""
    if dod == 0:
        w.write_bits(0, 1)
        return "0"
    for prefix, prefix_len, width in buckets:
        if fits(dod, width):
            w.write_bits(prefix, prefix_len)
            w.write_bits(dod & ((1 << width) - 1), width)
            return str(width)
    raise ValueError(f"delta of delta {dod} does not fit in {buckets[-1][2]} bits")


def read_dod(r, buckets):
    if r.read_bits(1) == 0:
        return 0
    ones = 1
    while ones < len(buckets) and r.read_bits(1) == 1:
        ones += 1
    width = buckets[ones - 1][2]
    raw = r.read_bits(width)
    return raw - (1 << width) if raw > (1 << (width - 1)) else raw

Each bucket is a prefix, the prefix length, and the payload width. PAPER_BUCKETS is the table from the paper (7, 9, 12 and 32 bits of payload). write_dod tries the buckets from the smallest to the largest and writes the first one that fits, after a single 0 bit for the common case of zero. read_dod counts the leading 1 bits to learn which bucket was used (10, 110, 1110 or 1111), reads the payload, and turns it back into a signed number.

The odd-looking range deserves a minute. Seven bits give 128 patterns, and the paper’s range of -63 to 64 contains exactly 128 values. The writer stores the value modulo 128 and the reader treats anything above 64 as negative: raw 0 to 64 mean 0 to 64, and raw 65 to 127 mean -63 to -1. So 64 is included and -64 is not. Prometheus’s source uses the same range: its bitRange function tests -((1<<(nbits-1))-1) <= x && x <= 1<<(nbits-1).

One deliberate simplification: before the second sample there is no earlier delta, so we start from a delta of 0 and the second sample pays for its full interval once. The paper instead stores the first delta against a block header aligned to a two-hour window, and Prometheus special-cases the second sample. Avoiding special cases keeps the code short, and it costs one extra 36-bit code per chunk (you will see it in Step 5).

Run step03_timestamps.py to price the paper’s example and the edges of every bucket:

"""step03_timestamps.py: price delta-of-delta values under the paper's bucket table."""
import gorilla
from bits import BitReader, BitWriter

BUCKETS = gorilla.PAPER_BUCKETS

print("The paper's example: deltas of 60, 60, 59 and 61 seconds")
times = [0, 60, 120, 179, 240]
prev_delta = 0   # before the second sample there is no earlier delta, so we start from 0
for earlier, later in zip(times, times[1:]):
    delta = later - earlier
    w = BitWriter()
    kind = gorilla.write_dod(w, delta - prev_delta, BUCKETS)
    print(f"  delta {delta:>3}  delta of delta {delta - prev_delta:>3}  -> {w.nbits:>2} bits (bucket {kind})")
    prev_delta = delta

print()
print("Bucket edges")
for dod in (0, 1, 64, 65, -63, -64, 256, 257, 2048, 2049, 2**31, 2**31 + 1):
    w = BitWriter()
    try:
        kind = gorilla.write_dod(w, dod, BUCKETS)
    except ValueError as error:
        print(f"  {dod:>11} -> ValueError: {error}")
        continue
    back = gorilla.read_dod(BitReader(w.getvalue()), BUCKETS)
    print(f"  {dod:>11} -> {w.nbits:>2} bits (bucket {kind:>2}), read back {back}")

python step03_timestamps.py

The paper's example: deltas of 60, 60, 59 and 61 seconds
  delta  60  delta of delta  60  ->  9 bits (bucket 7)
  delta  60  delta of delta   0  ->  1 bits (bucket 0)
  delta  59  delta of delta  -1  ->  9 bits (bucket 7)
  delta  61  delta of delta   2  ->  9 bits (bucket 7)

Bucket edges
            0 ->  1 bits (bucket  0), read back 0
            1 ->  9 bits (bucket  7), read back 1
           64 ->  9 bits (bucket  7), read back 64
           65 -> 12 bits (bucket  9), read back 65
          -63 ->  9 bits (bucket  7), read back -63
          -64 -> 12 bits (bucket  9), read back -64
          256 -> 12 bits (bucket  9), read back 256
          257 -> 16 bits (bucket 12), read back 257
         2048 -> 16 bits (bucket 12), read back 2048
         2049 -> 36 bits (bucket 32), read back 2049
   2147483648 -> 36 bits (bucket 32), read back 2147483648
   2147483649 -> ValueError: delta of delta 2147483649 does not fit in 32 bits

The first block reproduces the paper’s idea. Our first delta of delta is 60 because we start from zero (9 bits), then the steady 60 second delta costs 1 bit, and the wobbles of -1 and +2 cost 9 bits each. In the second block, 64 still fits the 7-bit bucket but 65 needs the 9-bit one, -63 fits but -64 does not, and each bucket edge behaves as the table promised. The last line is the important safety net: a delta of delta that cannot fit in 32 bits raises ValueError. Silently wrapping it around would corrupt every later timestamp in the chunk.

Check it worked

Compare your output with the block above. The rows at the bucket edges (64 and 65, -63 and -64) are the ones that go wrong first if a range limit is off by one.

Common mistake

Using the ordinary two’s complement range (-64 to 63) for the 7-bit field. It looks natural, but it disagrees with the reader above, and the value -64 would come back as +64. Step 6 plants exactly this bug to show that the tests catch it.

Step 4: Encode values by XOR-ing them with the previous value

To XOR two floats you first need their raw bits. struct.pack(">d", value) gives the eight bytes of a float. The Python documentation says that for the d format “the packed representation uses the IEEE 754 binary32, binary64 or binary16 format”, “regardless of the floating-point format used by the platform”, so the result is the same on every machine. Wikipedia’s description of binary64 lists a sign bit, an 11-bit exponent and a significand precision of “53 bits (52 explicitly stored)”.

Append this to gorilla.py:

# gorilla.py (continued)
def float_to_bits(value):
    return struct.unpack(">Q", struct.pack(">d", value))[0]


def bits_to_float(bits):
    return struct.unpack(">d", struct.pack(">Q", bits))[0]


class Window:
    """Leading and trailing zero counts of the last XOR that opened a new window."""

    def __init__(self):
        self.leading = None
        self.trailing = None


def write_xor(w, prev_bits, bits, window):
    """Write one value as its XOR with the previous value; return "same", "reuse" or "new"."""
    x = prev_bits ^ bits
    if x == 0:
        w.write_bits(0, 1)
        return "same"
    w.write_bits(1, 1)
    leading = min(64 - x.bit_length(), 31)    # the field is only 5 bits wide
    trailing = (x & -x).bit_length() - 1
    if window.leading is not None and leading >= window.leading and trailing >= window.trailing:
        w.write_bits(0, 1)
        w.write_bits(x >> window.trailing, 64 - window.leading - window.trailing)
        return "reuse"
    window.leading, window.trailing = leading, trailing
    meaningful = 64 - leading - trailing      # 1 to 64
    w.write_bits(1, 1)
    w.write_bits(leading, 5)
    w.write_bits(meaningful & 0x3F, 6)        # 64 does not fit in 6 bits, so it is stored as 0
    w.write_bits(x >> trailing, meaningful)
    return "new"


def read_xor(r, prev_bits, window):
    if r.read_bits(1) == 0:
        return prev_bits
    if r.read_bits(1) == 1:
        leading = r.read_bits(5)
        meaningful = r.read_bits(6) or 64     # a stored 0 means 64
        if leading + meaningful > 64:
            raise ValueError("corrupt stream: window is wider than 64 bits")
        window.leading, window.trailing = leading, 64 - leading - meaningful
    elif window.leading is None:
        raise ValueError("corrupt stream: window reused before it was defined")
    meaningful = 64 - window.leading - window.trailing
    return prev_bits ^ (r.read_bits(meaningful) << window.trailing)

write_xor handles the three cases from the table. It computes the XOR, counts the leading and trailing zeros, and checks whether the new meaningful bits fit inside the Window remembered from the last time a window was written. Two details are easy to miss, and both come from the field widths:

First, the leading-zero count is stored in 5 bits, so the largest value it can hold is 31. If the XOR has 32 or more leading zeros, the code pretends it has 31 and stores a few extra zero bits. Prometheus’s source has the same cap: “Clamp number of leading zeros to avoid overflow when encoding.” Second, the length of the meaningful part is stored in 6 bits, which holds 0 to 63, yet a length of 64 is possible. Length 0 can never occur (a zero XOR is handled earlier), so 64 is written as 0 and the reader turns a 0 back into 64. The Prometheus source explains the same trick: “But that value doesn’t actually fit into the 6 bits we have.”

See it on real bit patterns with step04_values.py:

"""step04_values.py: XOR two floats and price the result."""
import gorilla
from bits import BitWriter


def show(label, bits):
    text = f"{bits:064b}"
    print(f"  {label:<5} {text[0]} {text[1:12]} {text[12:]}")


a, b = 12.0, 24.0
xa, xb = gorilla.float_to_bits(a), gorilla.float_to_bits(b)
x = xa ^ xb
print("columns: sign, exponent (11 bits), fraction (52 bits)")
show("12.0", xa)
show("24.0", xb)
show("xor", x)
leading, trailing = 64 - x.bit_length(), (x & -x).bit_length() - 1
print(f"leading zeros {leading}, trailing zeros {trailing}, meaningful bits {64 - leading - trailing}")

print()
print("value   case   bits")
values = [12.0, 24.0, 24.0, 12.0, 24.0, 12.5, 12.5]
w, window, prev = BitWriter(), gorilla.Window(), gorilla.float_to_bits(values[0])
for v in values[1:]:
    bits = gorilla.float_to_bits(v)
    before = w.nbits
    kind = gorilla.write_xor(w, prev, bits, window)
    print(f"{v:>5}  {kind:<5} {w.nbits - before:>4}")
    prev = bits

python step04_values.py

columns: sign, exponent (11 bits), fraction (52 bits)
  12.0  0 10000000010 1000000000000000000000000000000000000000000000000000
  24.0  0 10000000011 1000000000000000000000000000000000000000000000000000
  xor   0 00000000001 0000000000000000000000000000000000000000000000000000
leading zeros 11, trailing zeros 52, meaningful bits 1

value   case   bits
 24.0  new     14
 24.0  same     1
 12.0  reuse    3
 24.0  reuse    3
 12.5  new     18
 12.5  same     1

The floats 12.0 and 24.0 have the same sign and fraction, and their exponents differ by one, so the XOR is a single 1 bit in the last exponent position: 11 leading zeros, 52 trailing zeros and one meaningful bit. Writing it as a new window costs 1 + 1 + 5 + 6 + 1 = 14 bits, which is the 14 bits the paper counts for the same shape of XOR. The second 24.0 is unchanged and costs 1 bit. Returning to 12.0 produces the same XOR pattern as before, so the window is reused and the cost is 1 + 1 + 1 = 3 bits. The value 12.5 needs a 1 bit at a lower position, so the old window no longer fits: a new window with 5 meaningful bits costs 18 bits. One more 12.5 costs a single bit.

Check it worked

The last column of costs should read 14, 1, 3, 3, 18 and 1. If the first one is not 14, the window bookkeeping in write_xor is the likely culprit.

Step 5: Put the pieces together into encode() and decode()

A compressed chunk has a simple layout: the sample count (32 bits), the first timestamp (64 bits, in two’s complement so negative times survive), the first value (64 bits, with no compression because there is nothing to compare it with), and then, for every later sample, a delta-of-delta code followed by an XOR code. Append the last part to gorilla.py:

# gorilla.py (continued)
def to_signed64(raw):
    return raw - (1 << 64) if raw >= (1 << 63) else raw


def encode(samples, buckets=PAPER_BUCKETS, ledger=None):
    """Compress [(timestamp_ms, value), ...] into bytes. If `ledger` is a list, it receives one
    (timestamp bucket, timestamp bits, value case, value bits) row per sample after the first."""
    w = BitWriter()
    w.write_bits(len(samples), 32)
    if not samples:
        return w.getvalue()
    t, v = samples[0]
    if not -(1 << 63) <= t < (1 << 63):
        raise ValueError("the first timestamp must fit in a signed 64-bit integer")
    w.write_bits(t & MASK64, 64)              # two's complement, so negative times survive
    prev_bits = float_to_bits(v)
    w.write_bits(prev_bits, 64)
    prev_t, prev_delta, window = t, 0, Window()
    for t, v in samples[1:]:
        delta = t - prev_t
        start = w.nbits
        ts_kind = write_dod(w, delta - prev_delta, buckets)
        middle = w.nbits
        bits = float_to_bits(v)
        val_kind = write_xor(w, prev_bits, bits, window)
        if ledger is not None:
            ledger.append((ts_kind, middle - start, val_kind, w.nbits - middle))
        prev_t, prev_delta, prev_bits = t, delta, bits
    return w.getvalue()


def decode(data, buckets=PAPER_BUCKETS):
    r = BitReader(data)
    count = r.read_bits(32)
    if count == 0:
        return []
    t = to_signed64(r.read_bits(64))
    bits = r.read_bits(64)
    out = [(t, bits_to_float(bits))]
    delta, window = 0, Window()
    for _ in range(count - 1):
        delta += read_dod(r, buckets)
        t += delta
        bits = read_xor(r, bits, window)
        out.append((t, bits_to_float(bits)))
    return out

The optional ledger argument is a learning and measuring aid: when you pass a list, encode adds one row per sample after the first, with the timestamp bucket and the number of bits it used and the value case and its bits. Run step05_ledger.py. It encodes the three samples from the paper’s Figure 2 and then eight samples of its own:

"""step05_ledger.py: encode a few samples and print what every sample cost."""
import gorilla

T0 = 1_760_000_000_000


def print_ledger(samples):
    ledger = []
    data = gorilla.encode(samples, ledger=ledger)
    print(f"{'sample':>6} {'timestamp bucket':>16} {'bits':>5}   {'value case':<10} {'bits':>5}")
    for number, (ts_kind, ts_bits, val_kind, val_bits) in enumerate(ledger, start=2):
        print(f"{number:>6} {ts_kind:>16} {ts_bits:>5}   {val_kind:<10} {val_bits:>5}")
    body_bits = sum(ts_bits + val_bits for _, ts_bits, _, val_bits in ledger)
    print(f"count and first sample: {32 + 64 + 64} bits; the rest: {body_bits} bits")
    print(f"total: {len(data)} bytes, against {len(samples) * 16} bytes raw")
    print("decoded equals input:", gorilla.decode(data) == samples)
    print()


print("Figure 2 of the paper: the values 12, 12 and 24 at 62 s, 122 s and 182 s")
print_ledger([(62, 12.0), (122, 12.0), (182, 24.0)])

print("Eight samples, 15 seconds apart, with 1 ms of drift in the last four")
offsets_and_values = [(0, 12.0), (15_000, 12.0), (30_000, 12.5), (45_000, 12.5),
                      (60_001, 13.0), (75_001, 13.0), (90_001, 12.5), (105_001, 12.5)]
print_ledger([(T0 + offset, value) for offset, value in offsets_and_values])

python step05_ledger.py

Figure 2 of the paper: the values 12, 12 and 24 at 62 s, 122 s and 182 s
sample timestamp bucket  bits   value case  bits
     2                7     9   same           1
     3                0     1   new           14
count and first sample: 160 bits; the rest: 25 bits
total: 24 bytes, against 48 bytes raw
decoded equals input: True

Eight samples, 15 seconds apart, with 1 ms of drift in the last four
sample timestamp bucket  bits   value case  bits
     2               32    36   same           1
     3                0     1   new           14
     4                0     1   same           1
     5                7     9   new           15
     6                7     9   same           1
     7                0     1   reuse          4
     8                0     1   same           1
count and first sample: 160 bits; the rest: 95 bits
total: 32 bytes, against 128 bytes raw
decoded equals input: True

The paper’s text says that the delta of delta in its example “is stored with a two bit header (‘10’), and the value is stored in seven bits, for a total size of just 9 bits”, and that the XOR example “is stored in fourteen total bits”. Its Figure 2 uses the values 12, 12 and 24. Our encoder, given the same values, spends 9 bits on the second timestamp and 1 on its value, then 1 bit and 14 bits on the third sample, which agrees with the paper. (Our second delta of delta is 60 where the paper’s is -2, because we have no block header to measure the first delta against, but both land in the 7-bit bucket.)

The second table shows the codec in action. Sample 2 pays 36 bits for its timestamp, the price of not yet knowing the interval. After that, steady timestamps cost 1 bit each. At sample 5 the clock drifts by one millisecond: the drift costs 9 bits when it appears and another 9 when the clock returns to normal at sample 6, so 1 ms of jitter costs 18 bits. The value column shows all three cases: same, new and reuse. Of the 32 bytes in total, 20 are the count and first sample, and the other seven samples take 95 bits, which is why tiny chunks compress poorly (Step 10 returns to this).

Check it worked

Both blocks end with decoded equals input: True. If you see False or an exception, recheck the three parts of gorilla.py against the listings.

Step 6: Prove it is lossless, then plant bugs to see whether the tests notice

A compressor has two obligations. It must be lossless: decoding must return exactly what went in. And it must be small: it must stay within the bit budget you designed. Tests for the first obligation are easy to imagine; tests for the second are easy to forget. Floating point adds a third trap: you cannot compare results with ==. Under ==, NaN is not equal to NaN, and -0.0 is equal to 0.0, so a test could miss a lost sign bit. The tests below compare 64-bit patterns instead.

Write the tests

The first group checks the bit tool, the bucket edges, and the range limits. Save this as test_gorilla.py:

"""test_gorilla.py: round trip, boundary and size tests for bits.py and gorilla.py."""
import random

import pytest

import gorilla
from bits import BitReader, BitWriter

PAPER, PROM = gorilla.PAPER_BUCKETS, gorilla.PROMETHEUS_BUCKETS
T0 = 1_760_000_000_000


def as_bits(samples):
    """Compare by bit pattern: under ==, NaN != NaN and -0.0 == 0.0."""
    return [(t, gorilla.float_to_bits(v)) for t, v in samples]


def round_trips(samples, buckets=PAPER):
    return as_bits(gorilla.decode(gorilla.encode(samples, buckets), buckets)) == as_bits(samples)


def test_bit_fields_round_trip():
    rng = random.Random(1)
    fields = [(width, rng.getrandbits(width)) for width in (rng.randint(1, 64) for _ in range(500))]
    w = BitWriter()
    for width, value in fields:
        w.write_bits(value, width)
    r = BitReader(w.getvalue())
    assert [(width, r.read_bits(width)) for width, _ in fields] == fields


def test_reader_stops_at_the_end():
    r = BitReader(b"\xff")
    r.read_bits(8)
    with pytest.raises(EOFError):
        r.read_bits(1)


@pytest.mark.parametrize("dod, bits", [
    (0, 1), (1, 9), (-63, 9), (64, 9), (-64, 12), (65, 12), (-255, 12), (256, 12),
    (-256, 16), (257, 16), (-2047, 16), (2048, 16), (-2048, 36), (2049, 36),
    (2**31, 36), (-(2**31 - 1), 36),
])
def test_paper_bucket_edges(dod, bits):
    w = BitWriter()
    gorilla.write_dod(w, dod, PAPER)
    assert w.nbits == bits
    assert gorilla.read_dod(BitReader(w.getvalue()), PAPER) == dod


@pytest.mark.parametrize("buckets", [PAPER, PROM])
def test_every_small_dod_round_trips(buckets):
    for dod in range(-5000, 5001):
        w = BitWriter()
        gorilla.write_dod(w, dod, buckets)
        assert gorilla.read_dod(BitReader(w.getvalue()), buckets) == dod


def test_paper_table_refuses_what_does_not_fit():
    for dod in (2**31 + 1, -(2**31)):
        with pytest.raises(ValueError):
            gorilla.write_dod(BitWriter(), dod, PAPER)
    w = BitWriter()
    gorilla.write_dod(w, 2**40, PROM)         # the 64-bit bucket takes it
    assert w.nbits == 68                      # Prometheus's format notes say its timestamp codes use 1 to 68 bits
    assert gorilla.read_dod(BitReader(w.getvalue()), PROM) == 2**40

The second group covers the value cases and the awkward inputs. It pins the cost of the paper’s Figure 2 example, the cap of 31 on leading zeros (neighboring doubles such as 1.0 and 1.0000000000000002 XOR to a single bit, 63 leading zeros), and the 64-bit window stored as zero (1 + 1 + 5 + 6 + 64 = 77 bits, which is the longest value code). It round-trips special floating-point patterns, including positive and negative zero, infinities, subnormals and NaNs. One of the NaNs is 0x7ff0000000000002, the pattern Prometheus uses as its stale marker; its source describes it as “StaleNaN is a signaling NaN, due to the MSB of the mantissa being 0.” The codec turns every float into raw bits before it does anything else, and on this machine the marker survives the trip in both directions. The last test cuts the data at every possible length and demands an error each time, so truncated data can never turn into wrong numbers.

Append this to test_gorilla.py:

# test_gorilla.py (continued)
def test_value_cases_and_their_cost():
    values = [12.0, 24.0, 24.0, 12.0, 24.0]
    ledger = []
    gorilla.encode([(T0 + i, v) for i, v in enumerate(values)], ledger=ledger)
    assert [(row[2], row[3]) for row in ledger] == [("new", 14), ("same", 1), ("reuse", 3), ("reuse", 3)]


def test_the_papers_figure_2_costs():
    ledger = []
    gorilla.encode([(62, 12.0), (122, 12.0), (182, 24.0)], ledger=ledger)
    assert ledger == [("7", 9, "same", 1), ("0", 1, "new", 14)]   # the bit lengths printed in Figure 2


def test_leading_zero_field_is_capped_at_31():
    samples = [(T0, 1.0), (T0 + 1, 1.0000000000000002)]   # neighbouring doubles: the XOR is a single 1 bit
    assert round_trips(samples)
    ledger = []
    gorilla.encode(samples, ledger=ledger)
    assert ledger[0][2:] == ("new", 46)       # 1 + 1 + 5 + 6 + 33 meaningful bits


def test_a_64_bit_window_is_stored_as_zero():
    samples = [(T0, gorilla.bits_to_float(0)), (T0 + 1, gorilla.bits_to_float(0x8000000000000001))]
    assert round_trips(samples)
    ledger = []
    gorilla.encode(samples, ledger=ledger)
    assert ledger[0][2:] == ("new", 77)       # 1 + 1 + 5 + 6 + 64 meaningful bits


SPECIAL = [0x0, 0x8000000000000000, 0x7FF0000000000000, 0xFFF0000000000000, 0x7FF8000000000000,
           0x7FF0000000000002, 0xFFF8000000000001, 0x0000000000000001, 0x7FEFFFFFFFFFFFFF]


def test_special_floats_survive_bit_for_bit():
    values = [gorilla.bits_to_float(bits) for bits in SPECIAL]
    assert round_trips([(T0 + i, v) for i, v in enumerate(values + values[::-1])])


def test_zero_one_two_negative_and_unsorted_timestamps():
    assert gorilla.decode(gorilla.encode([])) == []
    assert round_trips([(T0, 1.5)])
    assert round_trips([(T0, 1.5), (T0 + 5, 2.5)])
    assert round_trips([(-5_000, 1.0), (-4_000, 2.0), (-9_000, 3.0), (7, 4.0)])


def test_cut_off_data_is_an_error_not_wrong_numbers():
    data = gorilla.encode([(T0 + 15_000 * i, float(i * i)) for i in range(40)])
    for cut in range(len(data)):
        with pytest.raises((EOFError, ValueError)):
            gorilla.decode(data[:cut])

The third group adds random series and two size checks. The random generator mixes steady clocks, jitter, gaps and even timestamps that go backwards, and value styles from integers to random 64-bit patterns. The size checks matter as much as the round trips: a steady series with a constant value must take exactly the bytes that 160 + 37 + 998 × 2 bits need, and a value that alternates between 1.0 and 2.0 must reuse one window for every sample after the first. Append this part too:

# test_gorilla.py (continued)
def random_series(rng):
    series, t, value = [], T0 + rng.randint(-10**6, 10**6), 0.0
    step = rng.choice([1, 15, 1000, 15_000, 300_000])
    kind = rng.choice(["ints", "decimals", "noise", "flat", "bits"])
    for _ in range(rng.randint(0, 60)):
        t += step + rng.choice([0, 0, 0, 1, -1, step, -step // 2, 10**7])   # jitter, gaps, even steps backwards
        if kind == "ints":
            value += rng.randint(-5, 50)
        elif kind == "decimals":
            value = round(value + rng.random() - 0.5, 2)
        elif kind == "noise":
            value = rng.random() * 10 ** rng.randint(-5, 9)
        elif kind == "bits":
            value = gorilla.bits_to_float(rng.getrandbits(64))               # NaNs, infinities, subnormals
        series.append((t, value))
    return series


@pytest.mark.parametrize("buckets", [PAPER, PROM])
def test_random_series_round_trip(buckets):
    rng = random.Random(2026)
    for _ in range(400):
        assert round_trips(random_series(rng), buckets)


def test_a_steady_series_costs_two_bits_a_sample():
    data = gorilla.encode([(T0 + 15_000 * i, 3.0) for i in range(1000)])
    expected_bits = 160 + (36 + 1) + 998 * 2   # header and first sample, second sample, then 1 + 1 bits each
    assert len(data) == -(-expected_bits // 8)


def test_alternating_values_reuse_one_window():
    ledger = []
    gorilla.encode([(T0 + 15_000 * i, 1.0 if i % 2 == 0 else 2.0) for i in range(50)], ledger=ledger)
    assert [row[2] for row in ledger] == ["new"] + ["reuse"] * 48
    assert {row[3] for row in ledger[1:]} == {13}      # 1 + 1 + 11 meaningful bits

python -m pytest -q

................................                                         [100%]
32 passed in 0.18s

Plant bugs on purpose

All 32 tests pass. But how good are the tests? A useful way to find out is mutation testing: plant one bug at a time in gorilla.py, run the tests, and see whether they fail. Save mutations.py. It copies the files to a temporary folder for each bug, so your real gorilla.py is never touched. It runs the full suite and then only the tests whose names contain round_trip:

"""mutations.py: plant one bug at a time in gorilla.py and see whether the tests notice."""
import re
import shutil
import subprocess
import sys
import tempfile
from pathlib import Path

HERE = Path(__file__).parent
FITS = "-((1 << (width - 1)) - 1) <= x <= (1 << (width - 1))"
MUTANTS = [
    ("two's complement range for deltas", FITS, "-(1 << (width - 1)) <= x < (1 << (width - 1))"),
    ("first bucket 8 bits wide", "(0b10, 2, 7)", "(0b10, 2, 8)"),
    ("no cap on the leading-zero field", "min(64 - x.bit_length(), 31)", "(64 - x.bit_length())"),
    ("reader forgets that 0 means 64", "meaningful = r.read_bits(6) or 64", "meaningful = r.read_bits(6)"),
    ("window reuse ignores trailing zeros",
     "leading >= window.leading and trailing >= window.trailing", "leading >= window.leading"),
    ("window reuse needs strictly more zeros",
     "leading >= window.leading and trailing >= window.trailing",
     "leading > window.leading and trailing > window.trailing"),
    ("no range check on the widest bucket", "if fits(dod, width):", "if fits(dod, width) or width == buckets[-1][2]:"),
    ("first timestamp read as unsigned", "t = to_signed64(r.read_bits(64))", "t = r.read_bits(64)"),
]


def run_tests(folder, *extra):
    run = subprocess.run([sys.executable, "-B", "-m", "pytest", "-q", "--tb=no", "-p", "no:cacheprovider", *extra],
                         cwd=folder, capture_output=True, text=True)
    failed = re.search(r"(\d+) failed", run.stdout)
    return f"caught ({failed.group(1)} failed)" if failed else "SURVIVED"


source = (HERE / "gorilla.py").read_text(encoding="utf-8")
print(f"{'planted bug':<40} {'all tests':<18} round-trip tests only")
for name, old, new in MUTANTS:
    assert source.count(old) == 1, name
    with tempfile.TemporaryDirectory() as tmp:
        for f in ("bits.py", "test_gorilla.py"):
            shutil.copy(HERE / f, tmp)
        Path(tmp, "gorilla.py").write_text(source.replace(old, new), encoding="utf-8")
        print(f"{name:<40} {run_tests(tmp):<18} {run_tests(tmp, '-k', 'round_trip')}")

python mutations.py

planted bug                              all tests          round-trip tests only
two's complement range for deltas        caught (9 failed)  caught (1 failed)
first bucket 8 bits wide                 caught (6 failed)  SURVIVED
no cap on the leading-zero field         caught (1 failed)  SURVIVED
reader forgets that 0 means 64           caught (4 failed)  caught (2 failed)
window reuse ignores trailing zeros      caught (4 failed)  caught (2 failed)
window reuse needs strictly more zeros   caught (2 failed)  SURVIVED
no range check on the widest bucket      caught (1 failed)  SURVIVED
first timestamp read as unsigned         caught (1 failed)  SURVIVED

The full suite catches all eight bugs. The right-hand column is the lesson: five of the eight bugs slip past the round-trip tests alone. Two of them (the 8-bit first bucket and the strict window reuse) leave the data intact and just waste bits, so only the exact bit-cost tests notice. The other three (no cap on the leading-zero field, an unsigned first timestamp, and no range check on the widest bucket) only matter for inputs that random testing rarely produces, so only the targeted edge tests notice. Random fuzzing is a good net, but it does not replace the specific tests.

Common mistake

Treating decode(encode(x)) == x as the whole test. For a compressor, a version that stores everything in 77 bits per value passes that test and defeats the point.

Step 7: Measure it on real and synthetic data

Build the comparison

Now compare Gorilla with general-purpose compressors on six series: the two real USGS series and four synthetic ones. The synthetic series use a perfectly regular 15 second clock and fixed seeds: a cpu_gauge (a random walk between 0 and 100 with one decimal digit), a request_counter (whole numbers that only grow), a constant, and noise (random floats, the worst case). They are built only from random(), because the Python documentation guarantees that “The generator’s random() method will continue to produce the same sequence when the compatible seeder is given the same seed”, while other methods may change between Python versions. Append this to datasets.py:

# datasets.py (continued)
import random

START_MS = 1_760_000_000_000   # a fixed start, so no result depends on the clock


def regular_times(n, step_ms=15_000):
    return [START_MS + i * step_ms for i in range(n)]


def jittered_times(n, step_ms=15_000, jitter_ms=3, seed=1):
    rng = random.Random(seed)
    return [START_MS + i * step_ms + int(rng.random() * (2 * jitter_ms + 1)) - jitter_ms for i in range(n)]


def drop_every(samples, k):
    """Remove every k-th sample, like a scrape that failed."""
    return [s for i, s in enumerate(samples) if (i + 1) % k != 0]


def cpu_gauge(n, seed=11):
    """A random walk between 0 and 100 with one decimal, like a CPU percentage."""
    rng, x, out = random.Random(seed), 40.0, []
    for _ in range(n):
        x = min(100.0, max(0.0, x + (rng.random() - 0.5) * 3))
        out.append(round(x, 1))
    return out


def request_counter(n, seed=7):
    """A counter that only goes up, stored as a float the way metrics systems store it."""
    rng, total, out = random.Random(seed), 0, []
    for _ in range(n):
        total += int(rng.random() * 41)
        out.append(float(total))
    return out


def constant(n, value=3.0):
    return [value] * n


def noise(n, seed=5):
    rng = random.Random(seed)
    return [rng.random() * 100 for _ in range(n)]


def synthetic(n=2017, step_ms=15_000):
    """Name -> samples, all on a perfectly regular 15 second clock."""
    times = regular_times(n, step_ms)
    return {
        "cpu_gauge": list(zip(times, cpu_gauge(n))),
        "request_counter": list(zip(times, request_counter(n))),
        "constant": list(zip(times, constant(n))),
        "noise": list(zip(times, noise(n))),
    }

The baselines are zlib at level 9 and xz (the lzma module at its strongest preset), both on the plain 16-byte records, plus xz after splitting the data into two columns (all timestamps first, then all values), which is how you would help a general-purpose compressor. Save step07_measure.py:

"""step07_measure.py: Gorilla against raw, zlib and lzma on real and synthetic series."""
import lzma
import struct
import zlib

import datasets
import gorilla


def interleaved(samples):
    return b"".join(struct.pack("<qd", t, v) for t, v in samples)


def columnar(samples):
    n = len(samples)
    return struct.pack(f"<{n}q", *(t for t, _ in samples)) + struct.pack(f"<{n}d", *(v for _, v in samples))


def xz(data):
    return len(lzma.compress(data, preset=9 | lzma.PRESET_EXTREME))


usgs = datasets.usgs_series()
series = {"usgs flow": usgs["flow"], "usgs stage": usgs["stage"], **datasets.synthetic()}

print(f"{'series':<16} {'n':>5} {'raw':>6} {'gorilla':>7} {'ratio':>6} {'ts bits':>7} {'val bits':>8} {'same':>6} {'zlib':>6} {'xz':>6} {'xz cols':>7}")
for name, samples in series.items():
    ledger = []
    data = gorilla.encode(samples, ledger=ledger)
    assert gorilla.decode(data) == samples
    raw = interleaved(samples)
    ts_bits = sum(row[1] for row in ledger) / len(ledger)
    val_bits = sum(row[3] for row in ledger) / len(ledger)
    same = 100 * sum(row[2] == "same" for row in ledger) / len(ledger)
    print(f"{name:<16} {len(samples):>5} {len(raw):>6} {len(data):>7} {len(raw) / len(data):>5.1f}x "
          f"{ts_bits:>7.2f} {val_bits:>8.2f} {same:>5.1f}% {len(zlib.compress(raw, 9)):>6} {xz(raw):>6} {xz(columnar(samples)):>7}")

python step07_measure.py

series               n    raw gorilla  ratio ts bits val bits   same   zlib     xz xz cols
usgs flow         2017  32272    1678  19.2x    1.02     5.56  66.9%   8254   6340    5676
usgs stage        2017  32272    4688   6.9x    1.02    17.51  66.8%   8540   6572    5784
cpu_gauge         2017  32272   14191   2.3x    1.02    55.22   2.9%  10976   7344    7140
request_counter   2017  32272    3956   8.2x    1.02    14.60   3.1%  11967   7916    7392
constant          2017  32272     529  61.0x    1.02     1.00 100.0%   5947   4380    4476
noise             2017  32272   16645   1.9x    1.02    64.95   0.0%  22064  19672   19308

Read the results

The columns ts bits and val bits are the average bits per sample spent on timestamps and values, and same is the share of values identical to the previous one. Reading the table row by row:

Timestamps are nearly free everywhere. The clock is regular, so each sample after the second costs one bit and the second pays 36. The paper reports that “about 96% of all time stamps can be compressed to a single bit” on Facebook’s data; here it is all but one.

River flow: 19.2 times smaller. The 1,678 bytes (about 6.7 bits per sample including the header) beat the best general-purpose result, xz on columns at 5,676 bytes, by 3.4 times. Two thirds of the readings repeat the previous one (the paper’s equivalent figure is “Roughly 51% of all values are compressed to a single bit since the current and previous values are identical”), and the rest are whole numbers that move a little.

Gage height: 6.9 times. The same share of repeats, but when the level moves, a value such as 4.51 changes noisy low bits, so values cost 17.51 bits each on average.

CPU gauge: only 2.3 times, and xz wins. This random walk changes almost every sample (only 2.9 percent repeat) and its one-decimal values have no short binary form, so each XOR is wide. Gorilla needs 14,191 bytes; xz on columns needs 7,140, about half as much. This is the weak spot, and Step 8 fixes it.

Counter and constant: 8.2 and 61.0 times. Whole numbers that only grow compress well, and a constant costs about 2 bits per sample.

Noise: 1.9 times, from the timestamps alone. Random floats cost 64.95 bits each, more than the 64 bits a plain float takes, because most XORs open a new window and pay its 13 bits of overhead. Gorilla is not magic: it exploits similarity, and noise has none.

Check it worked

The gorilla column must match exactly (1678, 4688, 14191, 3956, 529 and 16645 bytes), because the data and the seeds are fixed. The zlib and xz columns come from the compression libraries bundled with your Python, so they may differ a little.

Step 8: Fix the weak spot by storing whole numbers

Why do decimals compress badly? A value like 0.1 has no exact binary form (it repeats forever, like one third in decimal), so 37.4 and 37.5 differ in many scattered fraction bits. A whole number like 374 or 375 is stored exactly, with a long run of trailing zeros, and that is the case the paper praised. The trick is to store tenths (or hundredths) as whole numbers and divide when you read. This is the same idea as keeping money in integer minor units; our tutorial on handling money in Python covers why.

The division gets back the identical float. It works because 451 / 100 and float("4.51") are both the double nearest to 4.51: dividing two exactly representable numbers is correctly rounded, and so is parsing the text. That holds while the whole number is below 2⁵³ and the scale is a power of ten up to 10²². The script verifies it on all 2,017 readings of both series. Save step08_scale.py:

"""step08_scale.py: store 4.51 feet as the whole number 451 hundredths of a foot."""
import datasets
import gorilla

stage = datasets.usgs_series()["stage"]
scaled = [(t, float(round(v * 100))) for t, v in stage]
print("dividing by 100 restores every reading exactly:", [(t, v / 100) for t, v in scaled] == stage)

gauge = list(zip(datasets.regular_times(2017), datasets.cpu_gauge(2017)))
gauge_scaled = [(t, float(round(v * 10))) for t, v in gauge]
print("same for the cpu gauge in tenths of a percent:", [(t, v / 10) for t, v in gauge_scaled] == gauge)
print()
print(f"{'stored as':<34} {'bytes':>6} {'bits/sample':>12}")
for label, samples in (("stage in feet (4.51)", stage),
                       ("stage in hundredths of a foot", scaled),
                       ("cpu gauge in percent (37.4)", gauge),
                       ("cpu gauge in tenths of a percent", gauge_scaled)):
    data = gorilla.encode(samples)
    print(f"{label:<34} {len(data):>6} {len(data) * 8 / len(samples):>12.2f}")

python step08_scale.py

dividing by 100 restores every reading exactly: True
same for the cpu gauge in tenths of a percent: True

stored as                           bytes  bits/sample
stage in feet (4.51)                 4688        18.59
stage in hundredths of a foot        1521         6.03
cpu gauge in percent (37.4)         14191        56.29
cpu gauge in tenths of a percent     3766        14.94

Storing the gage height in hundredths of a foot cuts it from 4,688 to 1,521 bytes, and the CPU gauge in tenths of a percent falls from 14,191 to 3,766 bytes. Against xz on columns (5,784 and 7,140 bytes), Gorilla now wins on both: 3.8 and 1.9 times smaller. The gauge goes from 2.3 times to 8.6 times smaller than the raw 32,272 bytes.

Check it worked

Both verification lines must print True. If one prints False, a reading was not an exact multiple of the scale, so check the scale factor first.

Caveats

You must keep the scale factor somewhere (in your metric’s metadata or schema), because the codec does not know about it. Whole numbers above 2⁵³ are not exact as floats. And any sum or rate computed from scaled values must undo the scale.

Step 9: Choose the timestamp table that matches your clock

The paper tuned its bucket widths for timestamps in seconds. Prometheus stores milliseconds, and its source comments on the difference: “Gorilla has a max resolution of seconds, Prometheus milliseconds. Thus we use higher value range steps with larger bit size.” Its table uses payloads of 14, 17, 20 and 64 bits. Our encode takes the table as an argument, so you can compare them directly. The script below tries three clocks on 2,000 samples: an exact 15 second clock, the same clock with up to 3 milliseconds of jitter in either direction, and the exact clock with every 20th scrape lost. It also tries a 30 day outage and, last, decodes bytes with the wrong table. Save step09_buckets.py:

"""step09_buckets.py: the paper's timestamp table against Prometheus's, for millisecond timestamps."""
import datasets
import gorilla

N = 2000
values = datasets.cpu_gauge(N)
exact = list(zip(datasets.regular_times(N), values))
scenarios = {
    "exact 15 s clock": exact,
    "15 s clock, jitter of +-3 ms": list(zip(datasets.jittered_times(N), values)),
    "exact clock, every 20th scrape lost": datasets.drop_every(exact, 20),
}


def timestamp_bits(samples, buckets):
    ledger = []
    gorilla.encode(samples, buckets, ledger)
    return sum(row[1] for row in ledger) / len(ledger)


print(f"{'timestamps':<38} {'paper 7/9/12/32':>16} {'prometheus 14/17/20/64':>23}")
for name, samples in scenarios.items():
    print(f"{name:<38} {timestamp_bits(samples, gorilla.PAPER_BUCKETS):>16.2f} "
          f"{timestamp_bits(samples, gorilla.PROMETHEUS_BUCKETS):>23.2f}")

print()
thirty_days_ms = 30 * 24 * 3600 * 1000
outage = [(datasets.START_MS, 1.0), (datasets.START_MS + 15_000, 1.0),
          (datasets.START_MS + 15_000 + thirty_days_ms, 1.0)]
for name, buckets in (("paper", gorilla.PAPER_BUCKETS), ("prometheus", gorilla.PROMETHEUS_BUCKETS)):
    try:
        data = gorilla.encode(outage, buckets)
        print(f"30 day outage, {name}: {len(data)} bytes, round trip ok: {gorilla.decode(data, buckets) == outage}")
    except ValueError as error:
        print(f"30 day outage, {name}: ValueError: {error}")

print()
data = gorilla.encode(exact[:50], gorilla.PAPER_BUCKETS)
try:
    wrong = gorilla.decode(data, gorilla.PROMETHEUS_BUCKETS)
    print("paper bytes read with the Prometheus table, same samples back:", wrong == exact[:50])
except (EOFError, ValueError) as error:
    print("paper bytes read with the Prometheus table:", type(error).__name__ + ":", error)

python step09_buckets.py

timestamps                              paper 7/9/12/32  prometheus 14/17/20/64
exact 15 s clock                                   1.02                    1.01
15 s clock, jitter of +-3 ms                       8.33                   14.73
exact clock, every 20th scrape lost                4.67                    2.99

30 day outage, paper: ValueError: delta of delta 2591985000 does not fit in 32 bits
30 day outage, prometheus: 32 bytes, round trip ok: True

paper bytes read with the Prometheus table: ValueError: corrupt stream: window is wider than 64 bits

With an exact clock the two tables tie at about 1 bit. With 3 ms of jitter the paper’s table wins clearly, 8.33 bits against 14.73, because Prometheus’s smallest bucket always spends 16 bits even for tiny deviations. Its own source says so in a comment: “This seems to needlessly jump to large bit sizes even for very small deviations from zero.” When scrapes are lost, the deltas of deltas become +15000 and -15000, which overflow the paper’s 12-bit bucket and cost 36 bits each, while Prometheus’s 17-bit bucket holds them: 2.99 bits against 4.67. Neither table is best everywhere; the paper says its limits “were selected by sampling a set of real time series from the production system”, and you should do the same with yours.

The outage line shows a failure the paper’s table cannot avoid with milliseconds. Its 32-bit last bucket tops out at 2³¹ milliseconds, about 24.8 days, so a 30 day gap raises ValueError instead of being stored. Prometheus’s 64-bit last bucket takes it. A loud error is the right behavior, and the planted bug in Step 6 that removed the range check would have wrapped the value silently.

The final line is a reminder that the table is part of the file format. Here the bytes written with one table and read with another were caught as a corrupt stream, but nothing in this format guarantees that, because it has no checksum. Store an identifier of the encoding next to every chunk. Prometheus does: its chunk format has an encoding field whose values include XOR and XOR2.

Check it worked

On the exact clock the two tables should tie, on the jittery clock the paper’s table should win, and with lost scrapes Prometheus’s should win. Only the paper’s table should raise ValueError on the 30 day outage.

Step 10: Cut the stream into chunks so a short query does not decode a week

Gorilla is sequential. Each sample depends on the previous timestamp delta, the previous value and the previous window, so to read sample number 1,500 you must decode samples 0 to 1,499 first. The paper names the trade-off: “queries that wish to read data over a short time range might need to expend additional computational resources on decoding data.” It settles on two-hour blocks, because “blocks that extend longer than two hours provide diminishing returns for compressed size.” Prometheus aims for a fixed number of samples per chunk; its source defines DefaultSamplesPerChunk = 120 with the comment “DefaultSamplesPerChunk provides a default target number of samples per chunk.”

The remedy is to cut a long series into chunks and keep an index of each chunk’s first and last timestamp outside the compressed bytes, so a query can skip chunks that cannot contain its range. Save chunks.py:

"""chunks.py: cut a long series into fixed-size Gorilla chunks and read a time range back."""
from dataclasses import dataclass

import gorilla


@dataclass
class Chunk:
    first_t: int
    last_t: int
    data: bytes


def build(samples, per_chunk):
    parts = (samples[i:i + per_chunk] for i in range(0, len(samples), per_chunk))
    return [Chunk(part[0][0], part[-1][0], gorilla.encode(part)) for part in parts]


def query(chunks, t_from, t_to):
    """Return the samples with t_from <= t <= t_to and how many samples had to be decoded."""
    rows, decoded = [], 0
    for chunk in chunks:
        if chunk.last_t < t_from or chunk.first_t > t_to:
            continue                      # the index says this chunk cannot hold the range
        for t, v in gorilla.decode(chunk.data):
            decoded += 1
            if t_from <= t <= t_to:
                rows.append((t, v))
    return rows, decoded

Then run step10_chunks.py, which measures the size cost of different chunk lengths and the decoding cost of reading the most recent hour (12 samples):

"""step10_chunks.py: how long a chunk should be, and what a short query costs."""
import chunks
import datasets

flow = datasets.usgs_series()["flow"]
print(f"{'samples per chunk':>17} {'chunks':>6} {'total bytes':>11} {'bytes/sample':>12}")
for per_chunk in (2, 10, 60, 120, 288, 2017):
    built = chunks.build(flow, per_chunk)
    total = sum(len(c.data) for c in built)
    print(f"{per_chunk:>17} {len(built):>6} {total:>11} {total / len(flow):>12.2f}")

print()
t_from, t_to = flow[-12][0], flow[-1][0]   # the most recent hour: 12 samples, 5 minutes apart
for per_chunk in (120, 288, 2017):
    rows, decoded = chunks.query(chunks.build(flow, per_chunk), t_from, t_to)
    print(f"{per_chunk:>5} per chunk: {len(rows)} samples wanted, {decoded} decoded")

python step10_chunks.py

samples per chunk chunks total bytes bytes/sample
                2   1009       25983        12.88
               10    202        6406         3.18
               60     34        2214         1.10
              120     17        1800         0.89
              288      8        1669         0.83
             2017      1        1678         0.83

  120 per chunk: 12 samples wanted, 97 decoded
  288 per chunk: 12 samples wanted, 289 decoded
 2017 per chunk: 12 samples wanted, 2017 decoded

Every chunk pays our 20-byte header and restarts its windows, so tiny chunks are expensive (12.88 bytes per sample at 2 samples per chunk, close to raw) and the cost flattens near 0.83 bytes per sample from 288 samples upward. At 120 samples per chunk the week costs 1,800 bytes, about 7 percent more than one big chunk (1,678), but reading the last hour decodes 97 samples instead of 2,017. With 288 samples (one day), the last hour happens to straddle a boundary: 2,017 samples are 7 full chunks plus one lonely sample, so the query decodes 289. Choosing the chunk length is a balance between compression and read cost.

Check it worked

Bytes per sample should fall as chunks grow, and the three query lines should report 97, 289 and 2,017 samples decoded.

How this compares with Prometheus today

What Prometheus does

The Prometheus file tsdb/chunkenc/xor.go says that its code “was largely written by Damian Gryski” as part of go-tsz, a Go implementation of the Gorilla paper. Its chunk layout, as documented in the repository’s chunk format notes, is a 16-bit sample count, the first timestamp as a varint, the first value as a 64-bit float, the second timestamp delta as an unsigned varint, and then delta-of-delta and XOR codes for every later sample. Compared with the code you wrote:

The value code is the same one you built: the notes say it “has between 1 bit and 77 bits”, which matches the 1-bit same-value code and the 77-bit full window your tests pinned. The timestamp code uses the 14, 17, 20 and 64 bit table, and the notes give its range as 1 to 68 bits, which is the 68 bits one of your tests asserts. Prometheus special-cases the first two samples with varints where we use one uniform path, and its 16-bit count means a chunk cannot exceed 65,535 samples (the source panics with “chunk capacity exceeded”).

The format is also still moving. The same notes document a newer XOR2 encoding in which, from the third sample on, “a joint control prefix encodes both the timestamp delta-of-delta (dod) and whether the value changed”. The configuration documentation says “XOR2 gives better disk compression than XOR for typical Prometheus workloads and can store start timestamps”, and it warns that chunks written as XOR2 cannot be read by older Prometheus versions or by tools that do not support the encoding yet, such as a Thanos sidecar. The changelog lists “Stabilize the XOR2 float chunk encoding” under 3.15.0 (September 24, 2026). According to the main-branch configuration documentation I read on October 28, 2026, the encoding setting is storage.tsdb.chunk_encoding.floats, and when it is absent the encoding stays xor unless a feature flag for XOR2 or start-timestamp storage is set. Everything you built here is the foundation of both encodings.

What I did not verify

I did not run Prometheus, and I did not decode a real Prometheus chunk with this code, so I make no claim that our bytes match its bytes. I compared our bit costs with the paper’s worked example, our cost ranges with Prometheus’s format notes, and our constants with its source. I also did not benchmark speed: the code is pure Python, written for clarity.

Confirm it works end to end

From the gorilla_lab folder, run these three commands in order:

python fetch_usgs.py
python -m pytest -q
python step07_measure.py

You should see the row count and a sha256 line, then 32 passed, then the table from Step 7 with the same byte counts on the usgs flow row (1,678 bytes for Gorilla, if your USGS data matches mine). Every assert gorilla.decode(data) == samples inside the measuring script passed before it printed a row, so every number in that table comes from a lossless round trip.

Where to go next

• Make it streaming. Wrap the state (prev_t, prev_delta, prev_bits, window and the BitWriter) in an Appender class with an append(t, v) method, so a chunk can grow one sample at a time as Prometheus does.
• Replace the hand-made fuzzing. Our Hypothesis tutorial shows how property-based testing finds the edge cases you did not think of.
• Compare with a general-purpose codec. Our Zstandard dictionary tutorial measures a different way of shrinking small records, and our LSM-tree tutorial shows where compressed chunks like these could live on disk.
• Read the real thing. Open Prometheus’s xor.go next to your gorilla.py; you will recognize every function. Our Prometheus histograms tutorial covers another side of the same ecosystem.
• Try your own metrics. Export a few hours of a series you care about, store its values as scaled integers where that fits, and see which row of Step 7 it resembles.

Sources

• Pelkonen et al., Gorilla: A Fast, Scalable, In-Memory Time Series Database, Proceedings of the VLDB Endowment, 2015. Short quotations only; the paper is published under a Creative Commons Attribution-NonCommercial-NoDerivs license.
• Prometheus: tsdb/chunkenc/xor.go, chunk format notes, tsdb/head.go, model/value/value.go, configuration reference and CHANGELOG, all read on the main branch on October 28, 2026, plus the storage documentation.
• go-tsz by Damian Gryski, named in the Prometheus source as the origin of its code.
• U.S. Geological Survey: copyrights and credits, and the Instantaneous Values service queried in Step 1 (site 01646500).
• Python documentation: struct and random; Double-precision floating-point format on Wikipedia.