Blame for hdc/hypervec.py

"""Binary hypervector algebra for HDCSTORAGE.

Everything in this machine is a *bit-packed hypervector*: a Python ``int``
interpreted as ``d`` bits, bit 0 being the least significant.  Representation
choices, and why:

* **No numpy.**  Dense HDC is normally an array of +/-1 floats or of packed
  uint64 words.  A Python ``int`` of 10 000 bits is one object whose bitwise
  ops run in C over whole machine words, so ``bind`` is a single ``^`` and
  ``popcount`` is ``int.bit_count()``.  Memory per vector is d/8 bytes plus
  object overhead, no dependency, and no per-element Python loop.
* **Binary spatter codes (BSC).**  bind = XOR (self-inverse, so unbind ==
  bind), bundle = per-bit majority, permute = rotate within the d-bit mask.

Algebraic facts this module relies on, used everywhere downstream:

    bind(a, b) ^ b == a                       (self-inverse)
    sim(bind(a, x), bind(b, x)) ~= sim(a, b)  (binding preserves similarity)
    sim(permute(a, k), permute(b, k)) == sim(a, b)
    sim(bundle(S), x) ~ sqrt(1/|S|)	for x in S  -> cleanup recovers members
    E[sim(a, b)] == 0 for independent random a, b (sd ~ 1/sqrt(d))

The last line is the reason d is large: it sets the noise floor that cleanup
memory and the instruction decoder have to beat.
"""

from __future__ import annotations

import hashlib
import math
from typing import Iterable, Sequence

__all__ = [
    "DEFAULT_DIM",
    "mask",
    "popcount",
    "hamming",
    "random_hv",
    "hvc",
    "bind",
    "unbind",
    "bundle",
    "permute",
    "permute_inv",
    "sim",
    "sim_matrix",
    "int_to_hv",
    "hv_to_int",
    "hv_bits_of_int",
    "hv_from_bits",
    "hv_to_bits",
    "bit_probe",
    "shift",
    "CODEC_THRESHOLD",
    "MAX_INT_BITS",
    "density",
    "bits",
    "bitplanes",
    "role",
    "ROLES",
    "clear_cache",
]

#: Default dimensionality.  10 000 bits is the classic Kanerva/VSA size: the
#: noise floor 1/sqrt(d) ~ 0.01 is far below the ~0.4 signal of a bundle
#: member, so 4-6 field records decode reliably while 1 250 bytes per vector
#: still fits a program image in a few hundred kilobytes.
DEFAULT_DIM = 10_000


# --------------------------------------------------------------------------
# primitive bit operations
# --------------------------------------------------------------------------
def mask(d: int = DEFAULT_DIM) -> int:
    """The all-ones mask covering ``d`` bits."""
    return (1 << d) - 1


def popcount(x: int) -> int:
    """Number of set bits."""
    return x.bit_count()


def hamming(a: int, b: int) -> int:
    """Hamming distance between two hypervectors."""
    return (a ^ b).bit_count()


def density(v: int, d: int = DEFAULT_DIM) -> float:
    """Fraction of set bits (should hover near 0.5 for random vectors)."""
    return popcount(v & mask(d)) / d


def bits(v: int, d: int = DEFAULT_DIM) -> str:
    """Bit string, most-significant bit first (the usual left-to-right view)."""
    return format(v & mask(d), "0%db" % d)


def bitplanes(v: int, d: int = DEFAULT_DIM, width: int = 64) -> list[str]:
    """Split a vector into fixed-width bit planes, WSB first."""
    s = bits(v, d)
    return [s[i : i + width] for i in range(0, len(s), width)]


# --------------------------------------------------------------------------
# deterministic symbol -> hypervector
# --------------------------------------------------------------------------
_CACHE: dict[tuple[str, int], int] = {}


def clear_cache() -> None:
    _CACHE.clear()


def hvc(name: str, d: int = DEFAULT_DIM) -> int:
    """Deterministic hypervector for a symbol name.

    The same name always maps to the same vector, in-process and across runs,
    on any machine.  That is what makes a stored program image portable: the
    symbolic layer (opcodes, register names, literals, tool names) is
    *regenerated* from names rather than serialised, and only the composed
    instruction vectors live in the storage file.

    Uses BLAKE2b expanded over a 4-byte block counter, then masked to d bits.
    """
    key = (name, d)
    cached = _CACHE.get(key)
    if cached is not None:
        return cached
    out = 0
    block = 0
    h = hashlib.blake2b(name.encode("utf-8"), digest_size=64)
    seed = h.digest()
    while out.bit_length() < d:
        chunk = hashlib.blake2b(
            seed + block.to_bytes(4, "little"), digest_size=64
        ).digest()
        out = (out << 512) | int.from_bytes(chunk, "big")
        block += 1
    out &= mask(d)
    _CACHE[key] = out
    return out


def random_hv(d: int = DEFAULT_DIM, rng=None) -> int:
    """Random balanced hypervector (about d/2 set bits)."""
    if rng is None:
        import random

        rng = random
    return rng.getrandbits(d) & mask(d)


# Roles: named slots used to build records.  BSC's canonical record is
#   record = bundle([bind(role_i, value_i) for i in fields])
# and a field is read back with bind(record, role_i) -> noisy value_i.
ROLES: dict[str, str] = {
    "op": "role/operation",
    "ra": "role/dest-register",
    "rb": "role/source-register-b",
    "rc": "role/source-register-c",
    "imm": "role/immediate",
    "tgt": "role/target-pc",
    "addr": "role/address",
    "val": "role/value",
    "name": "role/name",
    "cap": "role/capability",
    "doc": "role/documentation",
    "arity": "role/arity",
    "tool": "role/tool",
    "arg": "role/argument",
    "tag": "role/tag",
    "kind": "role/kind",
}


def role(name: str, d: int = DEFAULT_DIM) -> int:
    """Hypervector for a role slot (``"op"``, ``"ra"``, ...)."""
    return hvc(ROLES.get(name, name), d)


# --------------------------------------------------------------------------
# the algebra
# --------------------------------------------------------------------------
def bind(a: int, b: int, d: int | None = None) -> int:
    """Binding = XOR.  Involutive: ``bind(bind(a, b), b) == a``."""
    v = a ^ b
    return v & mask(d) if d is not None else v


def unbind(record: int, key: int, d: int | None = None) -> int:
    """Unbinding is binding (XOR is self-inverse); named for readability."""
    return bind(record, key, d)


def permute(v: int, k: int = 1, d: int = DEFAULT_DIM) -> int:
    """Cyclic bit rotation left by ``k`` (a rotation, so no bits are lost).

    Permutation is how this machine encodes *order*: position p of a sequence
    is ``permute(item, p)``, which keeps items near-orthogonal while making
    "same item, different position" a different address.
    """
    m = mask(d)
    v &= m
    k %= d
    if k == 0:
        return v
    return ((v << k) | (v >> (d - k))) & m


def permute_inv(v: int, k: int = 1, d: int = DEFAULT_DIM) -> int:
    """Inverse rotation (rotate right by ``k``)."""
    return permute(v, -k, d)


def _add_bitplane(plates: list[int], v: int) -> None:
    """Ripple-carry add of one bit-plane into a binary counter."""
    i = 0
    carry = v
    while carry:
        if i == len(plates):
            plates.append(carry)
            return
        nxt = plates[i] & carry
        plates[i] ^= carry
        carry = nxt
        i += 1


def _plane_ge(plates: Sequence[int], threshold: int, d: int) -> int:
    """Bit-sliced ``count >= threshold`` over the counter planes.

    Walks planes from most to least significant holding one predicate, ``eq``
    ("equal to the threshold so far"), and one accumulator, ``gt`` ("already
    greater, therefore certainly >= threshold").  Lanes that fall below the
    threshold at some position are dropped from ``eq`` and die there, which is
    exactly a ripple comparator built from AND/NOT — the same gate logic a
    hardware majority element uses, widened to d lanes at once.
    """
    m = mask(d)
    nbits = max(len(plates), threshold.bit_length()) or 1
    eq = m
    gt = 0
    for i in range(nbits - 1, -1, -1):
        p = plates[i] if i < len(plates) else 0
        if (threshold >> i) & 1:
            eq &= p  # count bit 0 where threshold bit is 1  ->  too small
        else:
            gt |= eq & p  # count bit 1 where threshold bit is 0  ->  greater
            eq &= ~p & m
        if eq == 0:
            break
    return (gt | eq) & m


def bundle(vecs: Iterable[int], d: int = DEFAULT_DIM) -> int:
    """Majority bundling: per-bit vote, ties resolve to 0.

    Implemented bit-sliced: every vector is added into a d-lane binary
    counter (one Python int per bit position of the count), then a single
    vectorised comparator thresholds it at ``floor(n/2)+1``.  Cost is
    O(n * log n) big-int operations rather than O(n * d) element visits,
    which is what keeps a 10 000-bit bundle at tens of microseconds.
    """
    items = list(vecs)
    if not items:
        return 0
    if len(items) == 1:
        return items[0] & mask(d)
    plates: list[int] = []
    for v in items:
        _add_bitplane(plates, v & mask(d))
    return _plane_ge(plates, len(items) // 2 + 1, d)


def sim(a: int, b: int, d: int = DEFAULT_DIM) -> float:
    """Cosine-style similarity in [-1, 1]: agreement of bits, rescaled.

    ``sim = (d - 2*hamming) / d``; +1 identical, 0 orthogonal, -1 complement.
    """
    if d <= 0:
        raise ValueError("d must be positive")
    return (d - 2 * (a ^ b).bit_count()) / d


def sim_matrix(vecs: Sequence[int], d: int = DEFAULT_DIM) -> list[list[float]]:
    """All-pairs similarity matrix."""
    return [[sim(a, b, d) for b in vecs] for a in vecs]


# --------------------------------------------------------------------------
# integers stored *as* hypervectors
# --------------------------------------------------------------------------
# A scalar is represented by superposition of bit-position atoms:
#   int_to_hv(n) = bundle([bind(NUM, permute(BITBASE, i)) for i where n has bit i])
# so an integer occupies a vector like every other value and can be recovered
# by cleanup (one similarity probe per bit position).  This is the pixel that
# makes "all state is hyperdimensional" true rather than decorative.
#
# The layout is **two's complement**, 32 bits: bit MAX_INT_BITS-1 is the sign,
# and a negative value is dense where a sign-magnitude encoding would be
# sparse.  That choice is what lets arithmetic happen as a ripple-carry
# circuit over the bit atoms (see hdc/adder.py): two's complement makes
# addition a pure function of the bits, with no sign case analysis anywhere.
_NUM = hvc("codec/integer")
_BITBASE = hvc("codec/bit-position")

#: Width of the integer payload, in bits, two's complement.
MAX_INT_BITS = 32

#: Similarity above which a bit atom counts as present in a bundle.
#:
#: Measured on this implementation (40 random atom subsets per size, both
#: dimensions): a member of a k-atom bundle recovers at 1/sqrt(k) — 0.35 at
#: k=4 down to 0.136 at k=32 — while a non-member stays below 0.044 at
#: d=10 000 and below 0.093 at d=2 000.  0.10 is the value that separates all
#: of those cases and round-trips 415/415 sampled 32-bit values at both
#: dimensions; d=2 000 is therefore the supported floor for the codec.
CODEC_THRESHOLD = 0.10

#: Smallest dimension at which the integer codec round-trips reliably.
#:
#: Measured over all 256 byte values plus 200 random bytes: 0 errors at
#: d=1500, 2000 and 3000; **60 errors at d=1000** and 401 at d=500, with a
#: worst delta of 512.  Below this floor the bit atoms crosstalk enough that a
#: decode silently returns a neighbouring value, which would corrupt every
#: threshold and shift downstream, so the codec refuses instead.
MIN_CODEC_DIM = 1500


def _check_dim(d: int) -> None:
    if d < MIN_CODEC_DIM:
        raise ValueError(
            "integer codec needs d >= %d (got %d): bit atoms crosstalk below that "
            "and values decode wrong" % (MIN_CODEC_DIM, d)
        )


def _bit_probe(i: int, d: int) -> int:
    return bind(_NUM, permute(_BITBASE, i, d), d)


def bit_probe(i: int, d: int = DEFAULT_DIM) -> int:
    """The atom for bit position ``i``, for probing a coded value directly.

    Exposed because other layers (the shader) need single-bit tests — sign,
    ordering — without decoding the whole integer.
    """
    return _bit_probe(i, d)


def int_to_hv(n: int, d: int = DEFAULT_DIM) -> int:
    """Encode a signed 32-bit integer as a hypervector (two's complement).

    Zero encodes to the all-zero vector (the empty bundle), so a cleared
    register is literally a vector with no atoms in it.

    This codec is *symbolic, not similarity-preserving*: two integers that
    share set bits share bit-position atoms, so ``sim(int_to_hv(5),
    int_to_hv(37))`` is high even though 5 and 37 have nothing to do with each
    other.  That is deliberate — decode needs the atoms to be shared and
    probeable, exactly like a bit field.  Do not read a distance between two
    encoded integers as a semantic relationship, and do not expect binding two
    of them to decorrelate the way binding two random vectors does.
    """
    _check_dim(d)
    lo, hi = -(1 << (MAX_INT_BITS - 1)), (1 << (MAX_INT_BITS - 1)) - 1
    if not lo <= n <= hi:
        raise ValueError("integer outside the %d-bit codec range: %r" % (MAX_INT_BITS, n))
    return hv_from_bits(hv_bits_of_int(n), d)


def hv_bits_of_int(n: int) -> list[int]:
    """The two's complement bit pattern of an integer, least significant first."""
    pattern = n & ((1 << MAX_INT_BITS) - 1)
    return [(pattern >> i) & 1 for i in range(MAX_INT_BITS)]


def hv_from_bits(bits: Sequence[int], d: int = DEFAULT_DIM) -> int:
    """Assemble a hypervector from bit values, purely by bundling atoms.

    No integer arithmetic is involved: each set position contributes its atom
    to a majority bundle.  Two's complement sign extension is the caller's
    concern, exactly as it is in hardware.
    """
    parts = [_bit_probe(i, d) for i, bit in enumerate(bits) if bit]
    return bundle(parts, d) if parts else 0


def hv_to_bits(v: int, d: int = DEFAULT_DIM, threshold: float = CODEC_THRESHOLD) -> list[int]:
    """Read the bit pattern back with one similarity probe per position."""
    return [1 if sim(v, _bit_probe(i, d), d) >= threshold else 0 for i in range(MAX_INT_BITS)]


def shift(v: int, k: int, d: int = DEFAULT_DIM) -> int:
    """Multiply or divide a coded integer by 2**k, without decoding it.

    This works only because the bit atoms are a *permutation orbit*:
    ``atom(i) = NUM ^ perm(BITBASE, i)``, so rotating the whole vector by k
    rotates every atom's index by k.  Rotation commutes with bundling (it is a
    bit permutation, and majority is per-bit), so

        rotate(hv, k) = bundle_i( perm(NUM,k) ^ atom(i+k) )

    and binding away the ``perm(NUM,k) ^ NUM`` offset leaves exactly the
    encoding of the shifted value.  Two XORs and two rotations, no similarity
    probe, no threshold, no decode.

    Positive ``k`` shifts left (multiply), negative shifts right (divide).
    Exact for non-negative values; for negative values the shift is logical,
    so the sign bit is not extended.
    """
    if k == 0:
        return v
    return bind(permute(v, k, d), permute(_NUM, k, d) ^ _NUM, d)


def hv_to_int(v: int, d: int = DEFAULT_DIM, threshold: float = CODEC_THRESHOLD) -> int:
    """Decode an integer hypervector by per-bit cleanup (two's complement)."""
    _check_dim(d)
    m = 0
    for i, bit in enumerate(hv_to_bits(v, d, threshold)):
        if bit:
            m |= 1 << i
    if m & (1 << (MAX_INT_BITS - 1)):
        m -= 1 << MAX_INT_BITS
    return m


# --------------------------------------------------------------------------
# diagnostics
# --------------------------------------------------------------------------
def capacity_bits(d: int = DEFAULT_DIM) -> float:
    """Information content of one vector, in bits (log2 of sign patterns)."""
    return float(d)


def orthogonal_capacity(d: int = DEFAULT_DIM, noise: float = 0.35) -> float:
    """Rough count of distinguishable near-orthogonal items (~ d / (2 log d))."""
    return d / (2.0 * math.log2(d))