diff --git a/.gitignore b/.gitignore index 4c42798..47e25db 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,9 @@ data/ literature/ Manifest.toml docs/build/ + +# Literate output: regenerated from docs/src/examples/*.jl on every docs build +docs/src/examples/*.md + +# local dev config +.claude/ diff --git a/CLAUDE.md b/CLAUDE.md index c09b447..0558a0f 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -26,7 +26,9 @@ test/ types.jl, operations.jl, encoding.jl, inference.jl, representations.jl ext_display.jl # rich-display tests, run in a separate process (extension loading is irreversible) docs/ - src/examples/ # Literate.jl tutorials (intro-to-hdc, dollar-of-mexico) + src/examples/ # Literate.jl tutorials (intro-to-hdc, encoding-data, + # colours, dollar-of-mexico, iris); *.md is Literate + # output, gitignored and regenerated on every build ``` ## Supported vector types (all subtypes of `AbstractHV{T} <: AbstractVector{T}`) @@ -50,7 +52,9 @@ Constructor convention (settled; each form has exactly one meaning, documented o - `HV(x)` — token shorthand for `encode(HV, x)`; `HV(n::Number)` **throws an ArgumentError** naming both alternatives (`D = n` / `encode(HV, n)`) — this replaced the old one-time `@warn` - Tuples of reals read as data, like vectors -`encode` is the canonical interface: `encode(HV, x; D)` is the deterministic token path (hash → seed), `encode(HV, x, strategy)` dispatches on `AbstractEncoding` strategies: `KMer(k)` (windows as atomic hashed tokens — resolves issue #53), `NGram(n)` (symbol-level shift-binding via `ngrams`), `Sequence()` (bundlesequence), `BagOfSymbols()` (multiset). New strategies = one struct + one `encode` method. Planned follow-up: stateful `AbstractEncoder{HV}` structs (RandomProjection, LevelEncoder) with `encode`/`decode`. +`encode` is the canonical interface: `encode(HV, x; D)` is the deterministic token path (hash → seed), `encode(HV, x, strategy)` dispatches on `AbstractEncoding` strategies: `KMer(k)` (windows as atomic hashed tokens — resolves issue #53), `NGram(n)` (symbol-level shift-binding via `ngrams`), `Sequence()` (bundlesequence), `BagOfSymbols()` (multiset). New strategies = one struct + one `encode` method. + +Stateful encoders live in the `AbstractEncoder{HV}` hierarchy (encode.jl): instances hold hypervector state built once at construction and slot in as the first argument of `encode`, with `decode` as the inverse. Two members: `LevelEncoder` encodes scalars against a shared level set (see below); `RandomProjection` encodes fixed-length real feature vectors (RGB triples, embeddings) through a shared `D × d` projection matrix: `z = R·x`, then a per-type nonlinearity (sign-like for binary/bipolar/ternary with a scalar-or-vector threshold `θ`, `β`-scaled identity/logistic/tanh for real/graded, `exp(im·β·z)` for FHRR). **Standardize features when they are on incommensurable scales — not as a reflex**: measured on iris, raw features beat standardized ones 0.98 vs 0.83, because there the largest-spread features are also the most discriminative (see the iris tutorial). Matrix flavours: `:gaussian` (default), `:bipolar`, `:sparse_ternary`. Constructors: `RandomProjection(HV, d; D, matrix, θ, β, seed, rng)`, `RandomProjection(HV, R::AbstractMatrix; θ, β)` (supplied matrix), and a data-driven ternary form `RandomProjection(TernaryHV, X; target_sparsity)` that solves θ from data at construction (immutable — construction-from-data, not fitting; `rethreshold(rp, θ)` returns a new encoder sharing `R`). RP's `decode(rp, hv, references)` is nearest-neighbour clean-up only (no analytic inverse exists — the nonlinearity is lossy); it errors without a reference set. FHRR + `:gaussian` is exactly random Fourier features: similarity approximates a Gaussian kernel with bandwidth `β`. Both encoders' FHRR paths route through the single internal `phase_encode(z, β)` helper — never duplicate it. Scalar `getindex` returns the element type `T`; non-scalar indexing returns a plain `Vector` of values, never a hypervector. Hypervectors are immutable (no `setindex!`). @@ -58,7 +62,7 @@ Scalar `getindex` returns the element type `T`; non-scalar indexing returns a pl | Operation | Function | Operator | Effect | |-----------|----------|----------|--------| -| Bundle | `bundle(hvs)` | `+` | Superposition; result similar to all inputs | +| Bundle | `bundle(hvs)` | `+`, `sum` | Superposition; result similar to all inputs | | Bind | `bind(hv1, hv2)` | `*` | Association; result dissimilar to inputs | | Unbind | `unbind(hv1, hv2)` | `/` | Inverse of bind (throws for `RealHV`: real MAP binding is not exactly invertible) | | Shift | `shift(hv, k)` / `ρ(hv, k)` | — | Circular shift by k positions | @@ -66,6 +70,8 @@ Scalar `getindex` returns the element type `T`; non-scalar indexing returns a pl In-place variants: `shift!`, `ρ!`, `perturbate!`. +**Bundling is m-way, not pairwise.** The rule depends on how many inputs there are, so bundling is *not* associative for BinaryHV/BipolarHV/RealHV/FHRR. `bundle(hvs)`, `sum(hvs)` and chained `x + y + z` (Julia parses a `+` chain as one variadic call) all reach `bundle` with the whole collection and are equivalent; `(x + y) + z`, `reduce(+, hvs)` and `foldl(+, hvs)` bundle a bundle, which biases the result towards the last inputs and destroys the property that a bundle is equally similar to each part. + ## Encoding strategies (compose primitives into structured representations) - `multiset(hvs)` — bundle a set of vectors @@ -75,11 +81,14 @@ In-place variants: `shift!`, `ρ!`, `perturbate!`. - `crossproduct(U, V)` — cross product of two sets - `ngrams(hvs, n)` — n-gram statistics for text/sequence encoding - `graph(sources, targets)` — directed/undirected graph encoding -- `level(hv, n)` / `encodelevel` / `decodelevel` / `convertlevel` — numeric level encoding + +Numeric values are encoded with `LevelEncoder` (encode.jl), which replaced the old `level`/`encodelevel`/`decodelevel`/`convertlevel` function family. It builds its level set **once** at construction (the old family rebuilt it per call, producing incomparable encodings). Two mechanisms, selected by dispatch: a perturbation **ladder** for all types (`LevelEncoder(HV, range, n; bandwidth = 2/n)`, flip fraction per step) and **fractional power encoding** for FHRR (`LevelEncoder(FHRR, values; β)`, continuous via `base^(β·x)`; pass an explicit level count `n` to get a ladder instead). `encode(lvl, x)` maps a number to a hypervector; `decode(lvl, hv)` maps back via nearest-neighbour similarity (FHRR also supports `method = :analytic` for continuous, clean-vector decoding). `LevelEncoder(levels, values)` wraps precomputed hypervectors — this labelled-codebook form is the seam for a future shared item-memory abstraction. ## Similarity and inference -- `similarity(u, v)` / `δ(u, v)` — type-dispatched similarity (cosine for bipolar/real, Jaccard for binary/graded, complex dot for FHRR) +- `similarity(u, v)` / `δ(u, v)` — type-dispatched similarity (cosine for bipolar/ternary/real/graded-bipolar/FHRR, Jaccard for binary/graded) +- `similaritymetric(HV)` — which metric a type uses (`:cosine` / `:jaccard`) +- `chancesimilarity(HV)` — what an *unrelated* pair scores: `0.0` under cosine, `1/3` under Jaccard. Needed to interpret any similarity number; assumes the type's default `distr` - `nearest_neighbor(u, collection)` — find closest match - `nearest_neighbor(u, collection, k)` — k-nearest neighbors @@ -152,10 +161,10 @@ julia --project -e 'using Pkg; Pkg.test()' - [x] ~~`bind`: no docstring (only `unbind` has one)~~ (added) - [ ] `shift` / `ρ`: no docstrings - [ ] `perturbate` / `perturbate!`: no docstrings -- [ ] `level`: docstring is minimal; does not explain the perturbation-based correlation mechanism +- [x] ~~`level`: docstring is minimal; does not explain the perturbation-based correlation mechanism~~ (obsolete: the family was replaced by `LevelEncoder`, whose docstring covers both mechanisms) - [ ] `three_pi` and `fuzzy_xor` helper functions: no docstrings - [ ] Internal functions (`aggfun`, `bindfun`, `neutralbind`, `noisy_and`, `elementreduce!`, `offsetcombine`, `empty_vector`, `eldist`, `vectype`): none documented -- [ ] `docs/make.jl`: Contents block references `examples.md` but actual pages are in `examples/` subdirectory +- [x] ~~`docs/make.jl`: Contents block references `examples.md`~~ (index.md rewritten as a real landing page) - [ ] No docstring for `^` on FHRR — exponentiation is an important FHRR feature - [ ] No developer docs on how to add a new HV type (what methods to implement, traits, etc.) diff --git a/README.md b/README.md index 2c13407..84ae9b9 100644 --- a/README.md +++ b/README.md @@ -147,7 +147,26 @@ Additionally, we provide common encoder strategies for different data structures - `crossproduct` - `ngrams` - `graph` -- `level` + +Numeric values are encoded with the stateful `LevelEncoder`, which builds a set of +level-correlated hypervectors once and shares it across every `encode`/`decode` call: + +```julia +lvl = LevelEncoder(BipolarHV, (0, 2π), 100) # 100 levels over the interval +hv = encode(lvl, π / 3) # hypervector for a numeric value +decode(lvl, hv) # ≈ π/3, via similarity matching +``` + +Fixed-length feature vectors (an RGB triple, an embedding — standardize them first!) +are encoded with `RandomProjection`, which projects through a matrix drawn once and +shared by every call; for `FHRR` this is exactly a random Fourier feature map, whose +similarity approximates a Gaussian kernel: + +```julia +rp = RandomProjection(BipolarHV, 3) # ℝ³ → 10,000-dimensional hypervectors +hv = encode(rp, [0.9, -0.2, 0.4]) # nearby inputs ↦ similar hypervectors +decode(rp, hv, references) # nearest-neighbour clean-up (needs a codebook) +``` Finally, the `similarity` function can be used to compare two hypervectors, by default using the best similarity metric for the hypervector type: diff --git a/TODO.md b/TODO.md index b3b87d8..c7e02e0 100644 --- a/TODO.md +++ b/TODO.md @@ -22,10 +22,54 @@ The package now has an explicit layer taxonomy: **primitives** (`operations.jl`) #53** (windows as atomic hashed tokens — genuinely different from `NGram(n)`, which shift-binds symbol encodings via `ngrams`); also `Sequence()` and `BagOfSymbols()`. Extension point: one struct + one `encode` method. -- [ ] Follow-up PR: stateful encoders as an `AbstractEncoder{HV}` hierarchy with - `encode`/`decode` — `RandomProjection` (fixed projection matrix) and - `LevelEncoder` (fixed level set; subsumes the §1.4b ladder problem). `encode` - deliberately keeps its first argument free for encoder instances. +- [x] Stateful encoders as an `AbstractEncoder{HV}` hierarchy with + `encode`/`decode` — **`LevelEncoder` landed (2026-07-15)**: builds its level + set once at construction (fixes §1.4b by construction), replaces and deletes + the `level`/`encodelevel`/`decodelevel`/`convertlevel` family. Ladder + (perturbation, all types, `bandwidth` keyword) + fractional power encoding + (FHRR, `β`), selected by dispatch. `encode` deliberately keeps its first + argument free for encoder instances. +- [x] `RandomProjection` as the next `AbstractEncoder{HV}` — **landed + (2026-07-15)**: fixed `D × d` projection matrix (`:gaussian`/`:bipolar`/ + `:sparse_ternary`), per-type nonlinearities, parametric scalar-or-vector + ternary threshold `θ` (+ data-driven `target_sparsity` constructor, + `rethreshold` sharing `R`), supplied-matrix constructor, FHRR = random + Fourier features via the shared internal `phase_encode(z, β)` helper (also + now used by `LevelEncoder`'s fractional power path — single source of + truth). Closes the colour/embedding (RGB / ESM-2-style feature-vector) + random-projection gap. +- [ ] Shared item-memory/codebook abstraction (labelled HV collection + + nearest-neighbour lookup): now seeded by `LevelEncoder(levels, values)` — its + `decode` is built on that stored collection. `RandomProjection`'s + `decode(rp, hv, references)` clean-up is the **second consumer** of the same + shape (currently via raw `nearest_neighbor` over a user-supplied reference + set); the tutorial's colour reverse-lookup and a future `train`/`predict` + need it too. Factor it out. + +--- + +## 0b. Documentation refresh (2026-07-27) + +Branch `docs-encoders-refresh` = `38-integration` + the encoders work, with the +tutorials brought up to the current API and the encoder layer documented: + +- `iris.jl` ported to `LevelEncoder` (it used all four removed level functions, and + defined local `encode`/`decode` helpers colliding with the new exports). +- `api.md` restructured into Combinators / Encoders; the stateful encoders reach the + API page for the first time. +- New tutorials: **Encoding data** (token path, KMer vs NGram, language recognition + from character k-mers — the issue #53 payoff) and **Colours** (RandomProjection end + to end, incl. decode-as-clean-up, the tuning knobs, and FHRR-as-random-Fourier- + features verified numerically in the page). +- Iris gained a second act encoding the same data with one `RandomProjection`. +- Landing page rewritten; developer guide added to the nav; `Remotes.GitHub` silences + the navbar warning; generated `docs/src/examples/*.md` untracked and gitignored + (they are Literate output, regenerated every build — they were a third of the diff). +- New API from review questions: `sum(hvs)` now bundles the whole collection, and + `similaritymetric` / `chancesimilarity` make similarity scores interpretable. + +Still unwritten from the plan: the amino-acids page and a "Why hyperdimensional?" +page (low- vs high-dimensional similarity histograms). Neither blocks anything. --- @@ -140,15 +184,21 @@ Remaining inconsistencies found by verification (not yet fixed): (`doctest(fix = true)`) — element signs flipped, so **the doctest CI job fails until this is done**. -### 1.4b Instance-path `convertlevel` builds encoder and decoder over DIFFERENT ladders +### 1.4b Instance-path `convertlevel` builds encoder and decoder over DIFFERENT ladders — FIXED (2026-07-15) +Fixed structurally by the `LevelEncoder` refactor (§0): the level set is built +once in the constructor and shared by every `encode`/`decode` call, and `seed` +makes the whole ladder deterministic. The old function family (including the +broken instance path) is deleted, no deprecation shims. Locked by the +"one shared level set" testset in `test/encoding.jl` (exact grid round-trips + +`!isdefined` checks for the removed names). Original notes: Found (by execution) while fixing §1.4, deliberately NOT fixed in that PR. `convertlevel(hv::AbstractHV, numvals)` calls `encodelevel(hv, ...)` and `decodelevel(hv, ...)`, each of which builds its own `level(hv, m)` ladder — and `level` perturbation is unseeded, so the two ladders share only the base vector. Measured: `decode(encode(x))` errors up to 1.0 (mean 0.41) on the instance path vs exactly 0.0 when both are built from one shared `level(...)` ladder. -- [ ] Fix: build the ladder once in `convertlevel` (and/or make `level` - deterministic given the base vector), then assert the roundtrip in tests. +- [x] Fix: build the ladder once (now: in the `LevelEncoder` constructor), then + assert the roundtrip in tests. ### 1.5c `perturbate` resamples from the TYPE-default distribution, not `hv.distr` — FIXED (2026-07-14) Instance-level `eldist(hv) = hv.distr` methods added for RealHV/GradedHV/ @@ -270,6 +320,33 @@ users an export conflict. `HV(this; D)` — no `rng` keyword; always `Xoshiro(hash(this))`. Breaking in output: every seeded/token vector changed; README + example outputs regenerated. +### 2.6 `RandomProjection(TernaryHV, ::AbstractMatrix)` positional collision (found 2026-07-15) + +The merged ternary constructor reads the positional matrix as **training data** +when `target_sparsity` is given and as a **supplied projection matrix** +otherwise — same positional shape, opposite meanings, disambiguated only by +keyword presence. The "ternary constructor: positional collision" testset in +`test/encoding.jl` pins what tests *can* pin: each path's documented reading, +and that a data-shaped (non-square) matrix misread as a projection matrix +cannot encode the intended features (immediate `DimensionMismatch` at first +use). What tests cannot make safe, deliberately not patched in that test-only +pass: + +- A **square** matrix is genuinely ambiguous: both readings yield a working + encoder, so a forgotten `target_sparsity` silently produces a projection + encoder built from data (and vice versa is undetectable). +- The misread's error is a downstream feature-length `DimensionMismatch`; it + never names the actual mistake (matrix misread as R instead of X). +- The ternary method's signature carries the data-path keywords, so the + supplied-matrix path **silently accepts and ignores** `D`, `matrix` and + `seed` (`RandomProjection(TernaryHV, R; D = 500)` returns a + `size(R, 1)`-dimensional encoder); the generic-type equivalent throws a + `MethodError`. Inconsistent and unlockable without blessing it. +- [ ] Fix by renaming one path rather than patching: a keyword-only or + distinctly named data-driven constructor (e.g. `fit_sparsity(TernaryHV, X; + target_sparsity, ...)` or `RandomProjection(TernaryHV, d; from_data = X)`), + and reject inapplicable kwargs on the supplied-matrix path. + --- ## 3. README (the first thing JuliaCon attendees will open) @@ -324,7 +401,7 @@ users an export conflict. is used" but `method` is a mandatory kwarg for plain vectors (`src/inference.jl:38`). Also `:hamming` returns a match *count*, not a similarity in [0,1] — document or normalize. - [ ] Typo `src/inference.jl:71`: "and `u`` " (stray backtick). -- [ ] `level` docstring: document the perturbation-based correlation mechanism and the FHRR `^`-based variant. +- [x] ~~`level` docstring: document the perturbation-based correlation mechanism and the FHRR `^`-based variant~~ (obsolete: the family was replaced by `LevelEncoder`, whose docstring documents both mechanisms with doctests). - [ ] Issue #36: the intro tutorial's ngrams example yields the wrong result — re-derive it. - [ ] Developer docs: how to add a new HV type (required methods: constructor, `eldist`, `empty_vector`, `bundle`, `bind`, `similarity`, traits…). @@ -339,7 +416,7 @@ Coverage gaps (each hid a §1 bug): approximate fuzzy recovery for the graded types, exact FHRR division, and the explicit `ArgumentError` for RealHV. - [ ] `perturbate` on FHRR (currently broken, §1.3). -- [ ] Instance-path `encodelevel`/`decodelevel`/`convertlevel` (currently broken, §1.4). +- [x] ~~Instance-path `encodelevel`/`decodelevel`/`convertlevel` (currently broken, §1.4)~~ (family deleted; `LevelEncoder` has its own testset: all 7 types, round-trips, bandwidth, FPE, precomputed constructor, bounds). - [ ] `bundle`/`bind` preserving custom `distr` (§1.5). - [ ] `similarity(...; method = :jaccard / :hamming)` — only `:cosine` is tested. - [ ] `isapprox` bootstrap path — the whole similarity testset is skipped for diff --git a/docs/Project.toml b/docs/Project.toml index a82f8d7..7ef2af1 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,10 +1,14 @@ [deps] +AlgebraOfGraphics = "cbdf2221-f076-402e-a563-3d30da359d67" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" -Handcalcs = "e8a07092-c156-4455-ab8e-ed8bc81edefb" HyperdimensionalComputing = "ebb7690e-f605-4e7a-86a0-727dbed6b293" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" +MLDatasets = "eb30cadb-4394-5ae3-aed4-317e484a6458" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" diff --git a/docs/make.jl b/docs/make.jl index 464c953..105656f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -27,7 +27,7 @@ repo_owner = split(repo_url, "/")[1] makedocs(; modules = [HyperdimensionalComputing], authors = "KERMIT research group and contributors", - repo = "https://github.com/$repo_url/blob/{commit}{path}#{line}", + repo = Remotes.GitHub(repo_owner, repo_name), sitename = "HyperdimensionalComputing.jl", format = Documenter.HTML(; prettyurls = get(ENV, "CI", "false") == "true", @@ -39,21 +39,16 @@ makedocs(; "HyperdimensionalComputing.jl" => "index.md", "Examples" => [ "Introduction to HDC" => "examples/introduction-to-hdc.md", + "Encoding data" => "examples/encoding-data.md", + "Colours: random projections" => "examples/colours.md", "What's the Dollar of Mexico?" => "examples/whats-the-dollar-of-mexico.md", + "Predictive modelling with HDC: Iris dataset" => "examples/iris.md", ], "API" => "api.md", + "Developer guide" => "developers.md", ], checkdocs = :exports, + warnonly = [:missing_docs], ) -deploydocs(; - repo = "github.com/$repo_url", - devbranch = begin - current_branch = get(ENV, "GITHUB_REF", "refs/heads/main") - dev_branch = if occursin("develop", current_branch) - "develop" - else - "main" - end - end, -) +deploydocs(; repo = "github.com/$repo_url") diff --git a/docs/src/api.md b/docs/src/api.md index 7300c32..b05c54f 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -6,36 +6,98 @@ CurrentModule = HyperdimensionalComputing This page contains the complete API reference for HyperdimensionalComputing.jl. -## Index +## Types -```@index -Pages = ["api.md"] +```@docs +AbstractHV +BinaryHV +BipolarHV +TernaryHV +GradedHV +GradedBipolarHV +RealHV +FHRR ``` -## Functions +## Operations -```@autodocs -Modules = [HyperdimensionalComputing] -Order = [:function] +```@docs +bundle +bind +unbind +shift! +shift +ρ +ρ! +normalize +normalize! +perturbate +perturbate! ``` -## Types +## Inference + +```@docs +similarity +δ +similaritymetric +chancesimilarity +nearest_neighbor +``` + +## Combinators + +Combinators take hypervectors and return a hypervector: they compose the primitive +operations into structured representations. -```@autodocs -Modules = [HyperdimensionalComputing] -Order = [:type] +```@docs +multiset +multibind +bundlesequence +bindsequence +hashtable +crossproduct +ngrams +graph ``` -## Constants +## Encoders -```@autodocs -Modules = [HyperdimensionalComputing] -Order = [:constant] +Encoders take *raw data* and return a hypervector. [`encode`](@ref) is the canonical +entry point: without a strategy it is the deterministic token path, and with an +[`AbstractEncoding`](@ref) strategy it composes the token path with the combinators +above. + +```@docs +encode +decode +AbstractEncoding +BagOfSymbols +Sequence +NGram +KMer +``` + +Stateful encoders hold hypervector state that is built once at construction — a level +set, a projection matrix — so that separately encoded values remain mutually +comparable. They are passed as the first argument of `encode`, and support the +inverse map `decode`. + +```@docs +AbstractEncoder +LevelEncoder +RandomProjection +rethreshold ``` -## Macros +## Package extensions + +`HyperdimensionalComputing.jl` has a couple of package extensions to interact with commonly used +Julia packages: + +### [UnicodePlots.jl](https://juliaplots.org/UnicodePlots.jl/stable/) -```@autodocs -Modules = [HyperdimensionalComputing] -Order = [:macro] +```@docs +unicodeheatmap +unicodehistogram ``` diff --git a/docs/src/examples/colours.jl b/docs/src/examples/colours.jl new file mode 100644 index 0000000..5e298b7 --- /dev/null +++ b/docs/src/examples/colours.jl @@ -0,0 +1,294 @@ +# # Colours: encoding continuous data with random projections +# +# The other tutorials encode *symbols* -- ingredients, characters, k-mers -- where the only +# question is whether two things are the same or different. Real data is often continuous: a +# measurement, an embedding, a pixel. Two colours can be *almost* the same, and we want the +# hypervectors to say so. +# +# The tool for this is [`RandomProjection`](@ref). It multiplies your `d`-dimensional input by a +# fixed random `D × d` matrix and applies a nonlinearity: +# +# ```math +# \mathbf{h} = f(R\,\mathbf{x}), \qquad R \in \mathbb{R}^{D \times d} +# ``` +# +# The [Johnson–Lindenstrauss lemma](https://en.wikipedia.org/wiki/Johnson%E2%80%93Lindenstrauss_lemma) +# guarantees that such a projection approximately preserves relative distances, so *close inputs +# give similar hypervectors*. Colours are the perfect testbed: an RGB triple is just a point in +# $[0,1]^3$, and we already have strong intuitions about which colours are alike. + +using HyperdimensionalComputing +using Colors, Random, Statistics, LinearAlgebra +Random.seed!(42); +nothing # hide + +# ## Encoding colours +# +# The projection matrix is drawn once, at construction, and stored inside the encoder. That is +# what makes separately encoded colours comparable -- they all pass through the same `R`. + +rp = RandomProjection(BipolarHV, 3; seed = 1) + +# `Colors.RGB` is a struct, so we unpack it into a plain 3-vector before encoding: + +col2vec(c) = [Float64(c.r), Float64(c.g), Float64(c.b)] +encodecolour(c) = encode(rp, col2vec(c)) + +# Let's take three colours -- two blue-greens and one orange: + +teal = RGB(23 / 255, 146 / 255, 153 / 255) + +# + +sky = RGB(4 / 255, 165 / 255, 229 / 255) + +# + +orange = RGB(254 / 255, 100 / 255, 11 / 255) + +# ## Comparing colours +# +# Encoding them and comparing gives exactly what our eyes say: + +similarity(encodecolour(teal), encodecolour(sky)) # both blue-green + +# + +similarity(encodecolour(teal), encodecolour(orange)) # opposite ends + +# Perceptual closeness has become hypervector similarity, and we never wrote a rule about +# colours. Compare this with the token path, where `teal` and `sky` would be two unrelated +# symbols with a similarity near zero. +# +# ## Blending by bundling +# +# Because the encoding is geometric, the HDC primitives now do geometric things. Bundling two +# colour hypervectors gives something *between* them, i.e., a blend: + +v_blend = bundle([encodecolour(teal), encodecolour(orange)]) +similarity(v_blend, encodecolour(teal)), similarity(v_blend, encodecolour(orange)) + +# The blend sits close to both parents, which is precisely the defining property of `bundle`. +# But what colour *is* it? To answer that we need to get back out of hyperspace. +# +# ## Decoding: clean-up, not inversion +# +# Here is an important limitation: the nonlinearity in +# ``f(R\mathbf{x})`` throws away magnitudes: `sign` maps a whole half-space to the same bit. A +# random projection therefore has **no analytic inverse**. [`decode`](@ref) +# performs *clean-up*, searching a set of reference hypervectors for the closest +# match. Calling it without references is an error. +# +# So we build a codebook: a few thousand random colours, encoded once. + +reference_colours = [RGB(rand(), rand(), rand()) for _ in 1:2000] +reference_hvs = encodecolour.(reference_colours); + +# `decode` returns a `(similarity, index, neighbour)` tuple, so we look the colour up by index: + +decodecolour(hv) = reference_colours[decode(rp, hv, reference_hvs)[2]] + +# Can we encode and decode to recover our original? Encode teal, decode it, and compare: + +decodecolour(encodecolour(teal)) + +# Not bit-identical, it is the nearest of 2000 random colours, and it is lossy by construction. +# However, it is unmistakably the same colour. And now we can finally see the blend from before: + +decodecolour(v_blend) + +# ## Learning from data: the average colour of a concept +# +# Time for something less toy-like. Suppose we observe *categories* paired with colours (a +# label and an observation) and want to learn what colour each category "is". This is the +# classic HDC associative memory: **bind** each observation to its label, **bundle** everything +# into one hypervector, then **unbind** to query. +# +# We use the named colours that ship with `Colors.jl` as our data: + +named(word) = [RGB((c ./ 255)...) for (n, c) in Colors.color_names if occursin(word, n)] +categories = Dict(:fire => named("red"), :water => named("blue"), :plant => named("green")) +length.(values(categories)) + +# Each observation becomes `label ⊗ colour`, and the whole dataset collapses into a *single* +# hypervector: + +memory = bundle( + [ + bind(encode(BipolarHV, label), encodecolour(c)) + for (label, colours) in categories for c in colours + ] +) + +# That one vector is the entire model. To ask "what colour is `:plant`?" we unbind the label and +# decode what remains: + +decodecolour(memory * encode(BipolarHV, :plant)) + +# + +decodecolour(memory * encode(BipolarHV, :fire)) + +# + +decodecolour(memory * encode(BipolarHV, :water)) + +# Green, red and blue. The bundle averaged each category's colours, and unbinding pulled the +# right average back out from one hypervector holding all three categories at once. +# +# ## Robustness: signal buried in noise +# +# Real data is messier than that. Suppose each observation comes with *one* correct colour and +# *two* random distractors, and we are not told which is which. This is +# [multiple-instance learning](https://en.wikipedia.org/wiki/Multiple_instance_learning), and it +# is a natural fit for HDC: random hypervectors are quasi-orthogonal, so noise tends to cancel +# while the consistent signal reinforces. + +noisy = bundle( + [ + bind( + encode(BipolarHV, label), + bundle([encodecolour(c), encodecolour(RGB(rand(), rand(), rand())), encodecolour(RGB(rand(), rand(), rand()))]) + ) + for (label, colours) in categories for c in colours + ] +) +decodecolour(noisy * encode(BipolarHV, :plant)) + +# + +decodecolour(noisy * encode(BipolarHV, :fire)) + +# Still green and still red, with two-thirds of every observation being pure noise. +# +# ## Tuning the projection +# +# [`RandomProjection`](@ref) exposes the knobs that matter. The `matrix` keyword picks the +# distribution of `R`: `:gaussian` (the default), `:bipolar` (entries `±1`, cheap), and +# `:sparse_ternary` (mostly zeros, cheapest): + +rp_bipolar = RandomProjection(BipolarHV, 3; matrix = :bipolar, seed = 1) +similarity(encode(rp_bipolar, col2vec(teal)), encode(rp_bipolar, col2vec(sky))) + +# For the sign-based types, `θ` is the threshold applied after projection, and it controls the +# sparsity of the result. Rather than rebuilding the encoder to change it, [`rethreshold`](@ref) +# returns a new encoder that *shares the same matrix*, so the encodings stay comparable: + +rp_shifted = rethreshold(rp, 0.5) +similarity(encode(rp_shifted, col2vec(teal)), encode(rp_shifted, col2vec(sky))) + +# !!! tip "Mind your feature scales" +# Colours are already on a common `[0, 1]` scale. Real feature vectors often are not, and a +# projection is dominated by whichever feature has the largest spread -- so features measured +# in grams and kilometres need to be put on a common footing first. +# +# Standardising is not an automatic win, though. When features already share a unit *and* +# their spread carries signal, equalising them throws information away; the *Iris dataset* +# tutorial measures a case where standardising makes the classifier markedly worse. +# +# ## The kernel connection +# +# One combination is worth singling out. With [`FHRR`](@ref) hypervectors and a Gaussian +# matrix, `RandomProjection` is *exactly* +# [random Fourier features](https://en.wikipedia.org/wiki/Random_feature): the similarity between +# two encoded points approximates a Gaussian (RBF) kernel with bandwidth `β`, +# +# ```math +# \delta(\mathbf{h}_x, \mathbf{h}_y) \;\approx\; \exp\!\left(-\tfrac{1}{2}\beta^2\|\mathbf{x}-\mathbf{y}\|^2\right). +# ``` +# +# We can check that directly. For random pairs of points, compare the measured hypervector +# similarity against the kernel value: + +β = 1.5 +rp_fhrr = RandomProjection(FHRR, 3; β = β, seed = 7, D = 20_000) + +comparison = map(1:6) do _ + x, y = rand(3), rand(3) + d = norm(x - y) + ( + distance = round(d, digits = 3), + hdc = round(similarity(encode(rp_fhrr, x), encode(rp_fhrr, y)), digits = 3), + kernel = round(exp(-β^2 * d^2 / 2), digits = 3), + ) +end + +# The two columns agree to about three decimals: + +comparison + +# This connects HDC to kernel methods: it gives an *explicit*, finite-dimensional +# feature map for a kernel that is normally only accessible implicitly through dot products. +# +# ## Sharpening a boundary +# +# Bundling all examples of a class into a prototype is the simplest HDC classifier, and it works +# well when classes are distinct. It struggles when they are not. Yellows and greens are +# neighbours in RGB space, so their prototypes sit close together: + +yellows, greens = named("yellow"), named("green") +proto_yellow, proto_green = bundle(encodecolour.(yellows)), bundle(encodecolour.(greens)) +similarity(proto_yellow, proto_green) + +# Classifying by nearest prototype gets most, but not all, of them right: + +prototype_accuracy = mean( + vcat( + [similarity(encodecolour(c), proto_yellow) > similarity(encodecolour(c), proto_green) for c in yellows], + [similarity(encodecolour(c), proto_green) > similarity(encodecolour(c), proto_yellow) for c in greens], + ) +) + +# We can do better by *retraining*: instead of averaging the examples once, iteratively nudge a +# weight vector whenever it misclassifies one. This is the classic +# [perceptron](https://en.wikipedia.org/wiki/Perceptron), and on hypervectors it is a few lines: + +function perceptron(positives, negatives; α = 1, maxiter = 50) + w = zeros(length(first(positives))) + for _ in 1:maxiter + errors = 0 + for v in positives + if dot(w, collect(v)) <= 0 + w .+= α .* collect(v) + errors += 1 + end + end + for v in negatives + if dot(w, collect(v)) >= 0 + w .-= α .* collect(v) + errors += 1 + end + end + errors == 0 && break + end + return w +end + +w = perceptron(encodecolour.(yellows), encodecolour.(greens)) + +# The retrained boundary separates the two classes far more cleanly: + +perceptron_accuracy = mean( + vcat( + [dot(w, collect(encodecolour(c))) > 0 for c in yellows], + [dot(w, collect(encodecolour(c))) < 0 for c in greens], + ) +) + +# !!! note "Retraining is not yet part of the package API" +# The perceptron above is plain Julia written against the hypervectors, not a package +# function. A proper `train`/`predict` workflow is planned; until then this is the pattern to +# copy. +# +# ## Summary +# +# * [`RandomProjection`](@ref) turns continuous vectors into hypervectors while preserving +# relative distances -- perceptual closeness becomes hypervector similarity. +# * It is **stateful**: the matrix is drawn once and shared, so only vectors encoded through the +# same encoder are comparable. Use [`rethreshold`](@ref) to change `θ` while keeping `R`. +# * [`decode`](@ref) is **clean-up against a codebook**, not inversion -- the nonlinearity is +# lossy by design. +# * Once data is in hyperspace, `bind`/`bundle`/`unbind` give you blending, associative memories +# and noise-robust multi-instance learning for free. +# * With `FHRR` and a Gaussian matrix you get random Fourier features, i.e. an explicit feature +# map for the Gaussian kernel. diff --git a/docs/src/examples/encoding-data.jl b/docs/src/examples/encoding-data.jl new file mode 100644 index 0000000..4b2fad1 --- /dev/null +++ b/docs/src/examples/encoding-data.jl @@ -0,0 +1,241 @@ +# # Encoding data: from raw data to hypervectors +# +# Every HDC application starts with the same question: *how to turn data into +# hypervectors?* This tutorial walks through the answer for the four important kinds of data: +# tokens, sequences, numbers, and feature vectors, and explains which +# encoder to reach for in each case. +# +# It helps to know that `HyperdimensionalComputing.jl` is organised in three layers: +# +# | Layer | What it does | Examples | +# |:-----------------|:------------------------------------------|:-------------------------------------------| +# | **Primitives** | hypervector in, hypervector out | `bundle`, `bind`, `shift`, `perturbate` | +# | **Combinators** | *collections* of hypervectors in, one out | `multiset`, `ngrams`, `hashtable`, `graph` | +# | **Encoders** | **raw data** in, hypervector out | `encode`, `LevelEncoder`, `RandomProjection`| +# +# The *"Introduction to HDC"* tutorial covers the first two layers. This one is about the third. + +using HyperdimensionalComputing + +# ## Tokens: one object, one hypervector +# +# The simplest thing you can do is treat an object as an atomic symbol, i.e., a *token*, and give +# it its own hypervector. That is what [`encode`](@ref) does when you call it without a strategy: +# it hashes the object and uses the hash to seed the vector. + +encode(BinaryHV, "cat") + +# Two things make this useful. First, it is **deterministic**: the same object always yields the +# same hypervector, in this session or the next, without you having to store a lookup table. + +encode(BinaryHV, "cat") == encode(BinaryHV, "cat") + +# Second, *different* objects yield **quasi-orthogonal** hypervectors. This means that they are as unrelated as +# two random vectors, which is exactly what you want from unrelated symbols: + +similarity(encode(BinaryHV, "cat"), encode(BinaryHV, "dog")) + +# Is that number small? It depends on the metric, and this is a common source of confusion. +# `BinaryHV` uses Jaccard similarity, for which *unrelated* vectors score $1/3$ -- not $0$. Two +# helper functions save you from having to remember which type uses what: + +similaritymetric(BinaryHV), chancesimilarity(BinaryHV) + +# + +similaritymetric(BipolarHV), chancesimilarity(BipolarHV) + +# So `0.33` means "unrelated" for a `BinaryHV` and "clearly related" for a `BipolarHV`. Always +# read a similarity against [`chancesimilarity`](@ref) for its type. + +# Any hashable object works, e.g., strings, symbols, characters, numbers, tuples: + +similarity(encode(BipolarHV, :cat), encode(BipolarHV, 42)) + +# !!! tip "`HV(x)` is shorthand" +# `BinaryHV("cat")` is shorthand for `encode(BinaryHV, "cat")`. The shorthand covers tokens +# only; everything else in this tutorial goes through `encode`. Note that a *number* is not +# accepted by the shorthand -- `BinaryHV(42)` throws, because it is ambiguous with the +# dimensionality argument. Use `encode(BinaryHV, 42)` or `BinaryHV(; D = 42)` to say which +# you meant. +# +# ## Sequences: `KMer` and `NGram` are different things +# +# A string *is* a hashable object, so the token path happily encodes it -- as one atom: + +encode(BinaryHV, "ACGTACGT") + +# That is often not what you want. `"ACGTACGT"` and `"ACGTACGA"` differ by a single letter, but +# as tokens they are completely unrelated. To encode a sequence *as a sequence*, pass a +# **strategy**. The package ships four, and the first two are easy to confuse, so it is worth +# being precise about the difference. +# +# [`KMer(k)`](@ref) slides a window of length `k` over the sequence and treats **each window as +# one atomic token**: it hashes the whole substring and bundles the results. This is the classic +# k-mer profile from genomics and text processing. + +encode(BinaryHV, "ACGTACGT", KMer(3)) + +# It is exactly equivalent to hashing each window yourself and bundling: + +encode(BinaryHV, "ACGT", KMer(3)) == multiset([encode(BinaryHV, w) for w in ["ACG", "CGT"]]) + +# [`NGram(n)`](@ref) does something genuinely different: it encodes each **symbol**, then binds +# the symbols of a window together with shifts to record their positions, then bundles the +# windows. It is the [`ngrams`](@ref) combinator applied to token-encoded symbols. + +encode(BinaryHV, "ACGT", NGram(3)) == ngrams([encode(BinaryHV, c) for c in "ACGT"], 3) + +# The two give different hypervectors with different properties: + +encode(BinaryHV, "ACGTACGT", KMer(3)) == encode(BinaryHV, "ACGTACGT", NGram(3)) + +# The practical difference is what *shares structure* with what. Under `KMer`, `"ACG"` and +# `"CGA"` are unrelated tokens: they merely happen to use the same letters. Under `NGram`, they +# are built from the same three symbol hypervectors, combined in a different order, so they +# retain a relationship. Pick `KMer` when the window is the meaningful unit (k-mer profiles, +# character n-gram fingerprints); pick `NGram` when the symbols are meaningful and you want +# position-aware composition. +# +# The remaining two strategies cover the extremes of "does order matter?": + +seq = [encode(BipolarHV, c) for c in "ACGT"] +encode(BipolarHV, "ACGT", Sequence()) == bundlesequence(seq) # order matters + +# + +encode(BipolarHV, "ACGT", BagOfSymbols()) == multiset(seq) # order does not + +# ## Example: recognising languages from character k-mers +# +# Here is what k-mer encoding buys you. Every language has a characteristic distribution of +# short character sequences: `"th"` and `"ing"` are common in English, `"ij"` and `"en"` in +# Dutch. We never tell the computer any of this: we just encode each text as a k-mer profile and +# let similarity do the work. +# +# Our corpus is the opening of the Wikipedia article on DNA in seven languages (including +# [West Flemish](https://en.wikipedia.org/wiki/West_Flemish), a regional variant of Dutch): + +texts = Dict( + "english" => "Deoxyribonucleic acid (DNA) is a polymer composed of two polynucleotide chains that coil around each other to form a double helix. The polymer carries genetic instructions for the development, functioning, growth and reproduction of all known organisms and many viruses. DNA and ribonucleic acid (RNA) are nucleic acids. Alongside proteins, lipids and complex carbohydrates, nucleic acids are one of the four major types of macromolecules that are essential for all known forms of life.", + "simple english" => "DNA, short for deoxyribonucleic acid, is the molecule that contains the genetic code of living organisms. This includes animals, plants, protists, archaea and bacteria. It is made up of two polynucleotide chains in a double helix. DNA is in each cell in the organism and tells cells what proteins to make. Mostly, these proteins are enzymes. DNA is inherited by children from their parents.", + "dutch" => "Desoxyribonucleinezuur, beter bekend als DNA, is het biologische macromolecuul dat in alle levende cellen de basis vormt van erfelijkheid. DNA is een zeer lang polymeer, en bevat de genetische instructies voor de ontwikkeling, het functioneren, de groei en de voortplanting van alle bekende organismen en vele virussen. DNA heeft een ingewikkelde chemische structuur.", + "west-flemish" => "De primaire structuur van DNA is de volgorde van de nucleotiedn. Een nucleotiede es ip zyn beurt ipgebouwd uut 3 bouwsteenn: de fosfoatgroep, nen vuufoekign suker, en een boase. Der zyn 4 soortn boasn: twee purines en twee pyrimidines. Der komn ook vele modificoaties voor.", + "german" => "Desoxyribonukleinsaure, meist kurz als DNA, seltener auch als DNS abgekurzt, ist eine aus unterschiedlichen Desoxyribonukleotiden aufgebaute Nukleinsaure. Sie tragt die Erbinformation bei allen Lebewesen und den DNA-Viren. Das langkettige Polynukleotid enthalt in Abschnitten von Genen besondere Abfolgen seiner Nukleotide.", + "french" => "L'acide desoxyribonucleique, ou ADN, est une macromolecule biologique presente dans presque toutes les cellules ainsi que chez de nombreux virus. L'ADN contient toute l'information genetique, appelee genome, permettant le developpement, le fonctionnement et la reproduction des etres vivants.", + "spanish" => "El acido desoxirribonucleico, conocido por las siglas ADN, es un acido nucleico que contiene las instrucciones geneticas fundamentales para el desarrollo, funcionamiento y reproduccion de todos los seres vivos y algunos virus; tambien es responsable de la transmision hereditaria.", +); + +# A little cleaning first -- lowercase, and drop punctuation so that we encode letters rather +# than typography: + +clean(s) = replace(lowercase(s), r"[^\p{L}\p{N}\s]" => "", r"\s+" => " ") + +# And now the whole model. One line per language: + +profiles = Dict(lang => encode(BinaryHV, clean(text), KMer(3)) for (lang, text) in texts); + +# That is the entire "training": no gradients, no iterations, one pass over each text. Let us +# rank every pair of languages by the similarity of their profiles: + +ranked = sort( + [ + (similarity(profiles[a], profiles[b]), a, b) + for a in keys(texts) for b in keys(texts) if a < b + ]; + rev = true, +) +first(ranked, 5) + +# The top of the list recovers real linguistic structure. English and Simple English are the most +# similar pair; Dutch and West Flemish come next, followed by the Romance pair (French/Spanish) +# and the Germanic pair (Dutch/German). Nothing about language families was encoded anywhere -- +# it falls out of character statistics alone. +# +# And the least similar pairs are the ones that cross family boundaries: + +last(ranked, 3) + +# !!! note "This is the k-mer workflow in general" +# Replace the texts with DNA reads and you have a sequence classifier; the code does not +# change. `encode(HV, sequence, KMer(k))` is the workhorse of HDC in genomics, and the +# reason `KMer` is a strategy in its own right rather than a special case of `NGram`. +# +# ## Numbers: `LevelEncoder` +# +# Tokens are the wrong tool for numbers. Hashing `3.0` and `3.1` gives two unrelated +# hypervectors, throwing away the one thing that matters about numbers: that nearby values are +# nearly the same. [`LevelEncoder`](@ref) fixes that by building a *ladder* of hypervectors in +# which neighbouring rungs are similar and distant rungs are not. + +lvl = LevelEncoder(BipolarHV, 0:0.5:10) + +# Encoding is done through the encoder object: + +similarity(encode(lvl, 3.0), encode(lvl, 3.5)) # neighbours: similar + +# + +similarity(encode(lvl, 3.0), encode(lvl, 9.0)) # far apart: quasi-orthogonal + +# Because the ladder lives *inside* the encoder, everything encoded through the same `lvl` is +# mutually comparable -- and [`decode`](@ref) can map a hypervector back to a number: + +decode(lvl, encode(lvl, 3.0)) + +# This is what makes `LevelEncoder` a **stateful** encoder: it is not a recipe you can re-derive, +# it is shared state that encoding and decoding must agree on. The *Iris dataset* tutorial uses +# it to encode flower measurements. +# +# ## Feature vectors: `RandomProjection` +# +# Finally, data that is already a vector of numbers, e.g., measurements, embeddings, pixel values, +# is handled by [`RandomProjection`](@ref), which multiplies your `d`-dimensional input by a +# fixed random `D × d` matrix and applies a nonlinearity. Distances are approximately preserved +# (the Johnson–Lindenstrauss lemma), so similar inputs give similar hypervectors. + +rp = RandomProjection(BipolarHV, 3; seed = 1) + +# Colours are a nice low-dimensional example: RGB triples where we already know which ones +# *should* be similar. + +teal = encode(rp, [0.09, 0.57, 0.6]) +sky = encode(rp, [0.02, 0.65, 0.9]) +orange = encode(rp, [0.99, 0.39, 0.04]) + +similarity(teal, sky) # both blue-green: similar + +# + +similarity(teal, orange) # opposite ends of the spectrum: dissimilar + +# Like `LevelEncoder`, this is stateful: the projection matrix is drawn once, and only vectors +# encoded through the *same* `rp` can be compared. The *Colours* tutorial explores this in depth, +# including how to decode back and how the different matrix types behave. +# +# ## Adding your own strategy +# +# The sequence strategies are deliberately thin -- each is a struct plus one `encode` method -- +# so adding your own takes about four lines. Suppose you want to encode only every second symbol: + +struct EveryOther <: AbstractEncoding end + +function HyperdimensionalComputing.encode(HV::Type{<:AbstractHV}, x, ::EveryOther; kwargs...) + return multiset([encode(HV, s; kwargs...) for s in collect(x)[1:2:end]]) +end + +encode(BinaryHV, "ACGTACGT", EveryOther()) == encode(BinaryHV, "AGAG", BagOfSymbols()) + +# Subtype [`AbstractEncoding`](@ref), add one method, and your strategy composes with +# everything else in the package. +# +# ## Summary +# +# | Your data | Use | +# |:-----------------------|:---------------------------------------| +# | symbols, categories | `encode(HV, x)` | +# | sequences (windows) | `encode(HV, x, KMer(k))` | +# | sequences (symbols) | `encode(HV, x, NGram(n))` | +# | ordered / unordered | `Sequence()` / `BagOfSymbols()` | +# | numbers | `LevelEncoder` | +# | feature vectors | `RandomProjection` | diff --git a/docs/src/examples/introduction-to-hdc.jl b/docs/src/examples/introduction-to-hdc.jl index 63f8b28..3b5976c 100644 --- a/docs/src/examples/introduction-to-hdc.jl +++ b/docs/src/examples/introduction-to-hdc.jl @@ -1,258 +1,331 @@ -# # Introduction +using Random #hide +# # Introduction to Hyperdimensional Computing (the HDC kitchen) # # Hyperdimensional Computing (HDC) is a brain-inspired computational paradigm that represents # and manipulates information using high-dimensional vectors called **hypervectors**. These -# vectors typically have thousands of dimensions (often 1,000-10,000), making them -# "hyperdimensional." -# -# The key insight is that high-dimensional spaces have unique mathematical properties that allow -# for robust, fault-tolerant computation. -# -# Let's start by loading the package in question, as follows: +# vectors typically have thousands of dimensions (often 1.000-10.000), making them +# "hyperdimensional." The key insight is that high-dimensional spaces have unusual mathematical +# properties that allow for robust, fault-tolerant computation based on a defined set of +# operations that enable representing any object or structure as a hypervector. +# +# Rather than listing those properties in the abstract, let's *cook* with them. In this tutorial +# we treat HDC as a kitchen: +# +# | HDC concept | Kitchen analogy | +# |:----------------------|:---------------------------------------------| +# | Mapping ($\varphi$) | Turning an ingredient into a hypervector | +# | Bundling ($\oplus$) | *Mixing* ingredients into a filling | +# | Binding ($\otimes$) | *Associating* an ingredient with a role | +# | Permutation ($\rho$) | *Ordering* the steps of a recipe | +# | Similarity ($\delta$) | Asking *"are these two plates alike?"* | +# +# Our running example is deliberately simple: we will "cook" a 🌮 **taco** and a 🍔 **hamburger**, +# and then add a third plate, a 🥪 **chicken club sandwich**, to see how the plates relate to +# each other -- and finally do some inference using algebraic operations over our "food" +# hypervectors. +# +# # Setting up our experiment +# +# First of all, we need to define the nature of our hypervectors. For this, we pick with +# *vector-symbolic architecture* (VSA) we will work in. This is, essentially, the flavour of +# hypervectors we will work with. For tutorial sake, we will use `BinaryHV`, the *binary spatter +# code* and arguably the most widely used VSA. using HyperdimensionalComputing -# # Creating hypervectors -# -# First, we will create a random bipolar hypervector. This is done as follows: - -BipolarHV() - -# As you may see, by default the hypervector created has 10.000 dimensions. This is the default -# value in `HyperdimensionalComputing.jl`, but one can can create a hypervector of any given -# dimensionality by providing the size of this as an argument: - -BipolarHV(; D = 8) - -# Alternatively, one can create a hypervector directly from a `Vector{T}` where `{T}` is an -# appropiate data type, e.g. integers for BipolarHV: +H = BinaryHV -BipolarHV(rand((-1, 1), 8)) +# Let's create a hypervector to see what one actually looks like: -# or you can directly pass any Julia structure to use it as a seed for the hypervector -# generation: +h = H() -BipolarHV(:foo) +# It is just a long vector of *bits*. For `BinaryHV` each component is a `0` or a `1`: -# Let's create 3 bipolar hypervector to use for the tutorial: +eltype(h) -h₁ = BipolarHV(; D = 8) -h₂ = BipolarHV(; D = 8) -h₃ = BipolarHV(; D = 8); +# ...and there are a lot of them -- this is where the "hyper" comes from: +length(h) -# The package has different hypervector types, such as `BipolarHV`, `TernaryHV`, `RealHV`, -# `GradedBipolarHV`, and `GradedHV`. All of this hypervectors have a common abstract type -# `AbstractHV` which can be used to build additional functions or encoding strategies (more on -# both later). +# !!! info "On hypervector types" +# We aliased `H = BinaryHV`, but the package offers several flavours -- `BipolarHV`, +# `TernaryHV`, `RealHV`, `GradedHV`, `GradedBipolarHV`, `FHRR` -- all sharing the abstract +# type `AbstractHV`. Type `?AbstractHV` in the REPL to see them, or `?BinaryHV` for the one +# we use here. Because everything below is written in terms of `H`, you can rerun the whole +# tutorial in another VSA by changing that single alias. By default a hypervector has 10.000 +# dimensions; pass `D` to change it, e.g. `H(; D = 8)`. # -# !!! info "On (abstract) types" -# All hypervectors implemented on `HyperdimensionalComputing.jl` can be found by checking the -# docstrings for the `AbstractHV` (by typing `?AbstractHV` on the Julia REPL). +# ## Operations # -# For more information on a specific hypervector type, the docstrings contain information on -# the implementation, operations, similarity measurement and other technical -# characteristics. +# Now that we have defined the nature of the hypervectors we will work with, let's go over each +# operation we have available to "cook" with hypervectors: # -# # Fundamental operations with hypervectors +# ### Mapping $\varphi$: _every ingredient is a hypervector_ # -# HDC uses three primary operations that preserve the hyperdimensional properties and allow for -# the representation more complex structures: +# The first rule of HDC is that *everything is a hypervector*. The mapping $\varphi$ takes any +# object -- a word, a number, or here an ingredient -- and assigns it a hypervector. The simplest +# mapping just draws a fresh random hypervector for each item. Let's stock our pantry, giving every +# ingredient its own random hypervector named (for fun) with its emoji: + +🥓 = H() # bacon +🥩 = H() # beef +🍞 = H() # bread +🍔 = H() # bun +🧀 = H() # cheese +🍗 = H() # chicken +🥬 = H() # lettuce +🥚 = H() # mayo +🧅 = H() # onion +🌶️ = H() # salsa +🍅 = H() # tomato +🫓 = H() # tortilla +🦃 = H(); # turkey + +# !!! tip "Seeding hypervectors" +# Each ingredient above is an *independent random draw*, so the exact numbers throughout this +# tutorial will differ every time you run it. When you instead need **reproducible** vectors +# or want the same object to always map to the same hypervector (e.g. so the token `"🥩"` maps +# to one fixed vector everywhere in a pipeline), you can *seed* a hypervector from any Julia +# object by passing it to the constructor. Additional constructors are available per +# hypervector type, which can be looked up in the API or in docstrings. # -# ## Bundling +# ## Similarity $\delta$: _are ingredients similar?_ # -# Bundling (also known as superposition) combines multiple hypervectors to create a new hypervector -# that is similar to it's constituents. +# Inference with hypervectors is based on comparing them, therefore we need some way of +# achieving this. Each VSA defines a similarity or distance measurement to assess how similar +# are hypervectors between each other. # -# $$u = [h_1 + h_2 + h_3]$$ +# Because each ingredient is an independent random draw, *different* ingredients are essentially +# unrelated ("quasi-orthogonal"). We can check this by comparing beef against a few ingredients at +# once. `similarity(🥩)` returns a *function* that measures similarity to 🥩, which we broadcast +# over a list: + +similarity(🥩).([🥩, 🧀, 🧅]) + + +# `BinaryHV` uses the **Jaccard** similarity, which runs from `0` to `1`. A hypervector is +# perfectly similar to itself (`1.0`), while two *unrelated* vectors share about a third of their +# set bits and so sit near a **baseline of ≈ 0.33**. # -# where $[...]$ denotes a potential normalization operations. In the case of bipolar -# hypervectors, this normalization operation is the `sign` function, which is defined as -# follows: # -# $$\text{sign}(i) = \begin{cases} -# +1 & \text{if } i > 0 \\ -# -1 & \text{if } i < 0 \\ -# 0 & \text{otherwise } -# \end{cases}$$ +# ### Bundling $\oplus$: _mixing_ # -# In HyperdimensionalComputing.jl, you can bundle hypervectors as follows: +# Bundling (a.k.a. superposition) combines hypervectors into a new one that is *similar to all +# of its ingredients*, think of tossing everything into one bowl. Let's mix a taco filling: -bundle([h₁, h₂, h₃]) +filling = bundle([🥩, 🧅, 🧀]) -# alternatively, you can use the `+` operator (which if overloaded for all `AbstractHV`): +# You can also use the overloaded `+` operator: -h₁ + h₂ + h₃ +filling == 🥩 + 🧅 + 🧀 -# This operation generates a hypervector that is similar to all it's constituent hypervectors, -# such that -# -# $$h₁ \sim u, h₂ \sim u, h₃ \sim u$$ -# -# where $\sim$ means that the hypervectors are similar, i.e. they share more components than -# expected by chance. -# -# ## Binding -# -# Binding combines multiple hypervectors to create a new hypervector that is dissimilar to it's -# constituents, such that: -# -# $$v = [h₁ \times h₂ \times h₃]$$ -# -# where $[...]$ represents a normalization procedure. +# The operation "remembers" what went into it: it is clearly similar to each of its ingredients, but not +# to something we never added (bread is a stranger to this bowl): + +similarity(filling).([🥩, 🧅, 🧀, 🍞]) + +# ### Binding $\otimes$: _associating_ # -# In HyperdimensionalComputing.jl, you can bind hypervectors as follows: +# Binding combines hypervectors into a new one that is *dissimilar to its inputs*. It is the +# tool for associating hypervectors to create new concepts. Let's define a `ROLE` hypervector +# and bind cheese to it: -bind([h₁, h₂, h₃]) +ROLE = H(:role) +topping = ROLE * 🧀 -# alternatively, you can use the `*` operator (which if overloaded for all `AbstractHV`): +# The resulting `topping` sits back at the ~0.33 baseline against both the role and the ingredient, +# binding *hides* its operands, so `topping` looks unrelated to either: -h₁ * h₂ * h₃ +similarity(topping).([🧀, ROLE]) -# This operation generates a hypervector that is similar to all it's constituent hypervectors, -# such that -# -# $$h₁ \nsim v, h₂ \nsim v, h₃ \nsim v$$ -# -# where $\nsim$ means that the hypervectors are dissimilar, i.e. they are quasi-orthogonal. -# -# ## Permutation -# -# Permutation (also known as shifting) is a special case of binding that creates a variant of a -# single hypervector via, generally speaking, a circular vector shifting with one or more -# positions. -# -# $$m = \rho(h₁)$$ +# Notably, binding is *reversible* in some VSAs. For `BinaryHV` the bind is a bitwise **XOR**, which +# is its own inverse, so binding again with the role recovers the ingredient *exactly* (similarity `1.0`). +# This "unbinding" is what will later let us *query* a recipe: -h₄ = RealHV(collect(0.0:9.0)) -h₄.v +similarity(ROLE * topping, 🧀) +# ## Permutation $\rho$: _ordering_ # +# Permutation takes a hypervector and cyclically shifts it into a new one that +# is dissimilar to the original. It is how HDC encodes order, e.g. because in the kitchen, order +# matters (sear *then* simmer is not the same as simmer *then* sear): -ρ(h₄).v +similarity(🥩, ρ(🥩)) -# -# The new hypervector will be, in principle, dissimilar to it's original version, such that: -# -# $$h_1 \nsim \rho(h_1) \nsim \rho\rho(h_1) \nsim \rho\rho(h_1) ...$$ -# -# where $\nsim$ means that the hypervectors are dissimilar, i.e. they are quasi-orthogonal. -# -# In `HyperdimensionalComputing.jl`, one can shift hypervector as follows: +# Applying it repeatedly keeps producing fresh, quasi-orthogonal vectors, giving each position its +# own signature: -ρ(h₁, 1) +similarity(🥩).([🥩, ρ(🥩, 1), ρ(🥩, 2), ρ(🥩, 3)]) -# +# We can use this to encode order or hierarchy, as this produces an alternative version of our original +# hypervector. Encode a two-step recipe by permuting the second step +# once (position 0, then position 1), and compare it to the same steps performed in the opposite +# order: -h₁ != ρ(h₁, 1) != ρ(h₁, 2) != ρ(h₁, 3) +sear = H(:sear) +simmer = H(:simmer) -# ## Similarity +similarity(sear + ρ(simmer), simmer + ρ(sear)) + +# Same two actions but with different order. Hypervectors come out unrelated. # -# Althought technically not an operation, in order to retrieve information from hypervectors, -# we need to compare them using similarity/distance functions. `HyperdimensionalComputing.jl` -# provides a handy `similarity` function that accepts: +# ## Encoding: _cooking a plate_ # -# 2 hypervectors: +# We now combine the operations to "cook." There is no single right way to turn a list of +# ingredients into a plate hypervector: the choice of *encoder* determines what the resulting +# vector remembers. The package ships several; here we compare three, from least to most +# structured, using a taco's ingredients: -similarity(h₁, h₂) +ingredients = [🫓, 🥩, 🧅, 🌶️, 🧀] -# A vector of hypervectors: +# **`multiset` creates an unordered bag**, it simply bundles the ingredients. It is +# the simplest encoder, but it forgets everything except which ingredients are present: shuffle +# them and you get the exact same vector. -similarity(h₁, h₁) +multiset(ingredients) -# or a hypervector and a vector of hypervectors: +# -similarity.(Ref(h₁), [h₁, h₂, h₃]) +similarity(multiset(ingredients), multiset(shuffle(ingredients))) -# `δ` is a synonim of `similarity`, and can also be used to create a function for similarity -# comparison, e.g. +# **`bundlesequence` creates an ordered stack**, it permutes each +# ingredient by its position before bundling, encoding ordered presence of hypervectors. -f = δ(h₁) +bundlesequence(ingredients) -# +# Now reversing the stack +# gives an unrelated vector, useful for layered dishes or recipe steps, where sequence matters: -f.([h₁, h₂, h₃]) +similarity(bundlesequence(ingredients), bundlesequence(reverse(ingredients))) - -# ## Encoding things as hypervectors -# -# The true power of HDC emerges when we combine the fundamental operations to encode complex data -# structures as hypervectors. By creatively applying bundling, binding, and shifting, we can -# represent virtually any type of information - from sequences and hierarchies to graphs and -# associative memories. The operations act as building blocks that can be composed in countless -# ways, limited only by our imagination and the specific requirements of our application. Let's -# explore some fundamental encoding strategies that demonstrate this flexibility. -# -# ### Key-value pairs +# **`hashtable` creates a a keyed record**, it binds each key-value pair and bundles them +# together. This is the most structured of the three: order is irrelevant, but each +# ingredient is filed under the role it plays, so we can later *query it back*. # -# Animal hypervectors: +# Let's define our roles: -H_dog = TernaryHV(:dog) -H_cat = TernaryHV(:cat) -H_cow = TernaryHV(:cow) -H_animals = [H_dog, H_cat, H_cow] +BASE = H(:base) # tortilla, bun, bread... +PROTEIN = H(:protein) # beef, chicken, turkey... +VEGGIE = H(:veggie) # onion, lettuce... +SAUCE = H(:sauce) # salsa, ketchup, mayo... +EXTRA = H(:extra) # cheese, bacon... -# Sound hypervectors: +roles = [BASE, PROTEIN, VEGGIE, SAUCE, EXTRA] -H_bark = TernaryHV(:bark) -H_meow = TernaryHV(:meow) -H_moo = TernaryHV(:moo) -H_sounds = [H_bark, H_meow, H_moo] +# Our 🌮 **taco** -- a tortilla base, beef, onion, salsa, and a bit of cheese -- and our +# 🍔 **hamburger** -- a bun, beef, lettuce, tomato, and cheese: +taco = hashtable(roles, [🫓, 🥩, 🧅, 🌶️, 🧀]) +burger = hashtable(roles, [🍔, 🥩, 🥬, 🍅, 🧀]) -# Associative memory: +# Each plate is now a *single* hypervector encoding its whole (structured) recipe. The three +# encoders trade off resolving power against simplicity: `multiset` answers only *"what is in +# it?"*, `bundlesequence` also captures *"in what order?"*, and `hashtable` captures *"what plays +# which role?"* -- the one we need to reason about recipes. +# +# !!! tip "Combinators and encoders" +# `multiset`, `bundlesequence` and `hashtable` are three of the built-in **combinators**: +# they take hypervectors and return a hypervector. The package also provides `multibind`, +# `bindsequence`, `ngrams`, `graph` and `crossproduct`. +# +# One layer up sit the **encoders**, which turn *raw data* into hypervectors: `encode` with a +# strategy such as `KMer` or `NGram` for sequences, and the stateful `LevelEncoder` (numbers) +# and `RandomProjection` (feature vectors). See the [API reference](../api.md) for the full +# catalogue. +# +# ## Comparison: are two plates alike? +# +# With every plate living in the same space, we can reason about them in two complementary ways: +# by **measuring similarity** and by doing **algebra** on the hypervectors. +# +# ### Measuring similarity +# +# Let's add a third plate, a 🥪 **club sandwich**. This one is interesting: its protein could be +# chicken *or* turkey. We express that ambiguity directly by **superposing** (bundling) the two +# poultry options into a single hypervector that is similar to both: -memory = (H_dog * H_bark) + (H_cat * H_meow) + (H_cow * H_moo); +poultry = 🍗 + 🦃 -# !!! note -# Alternatively you can use the `hashtable` encoder to achieve the same: +# The sandwich is then bread, that poultry, lettuce, mayo, and bacon: -memory == hashtable(H_animals, H_sounds) +sandwich = hashtable(roles, [🍞, poultry, 🥬, 🥚, 🥓]) -# Querying memory to search for dog's sound: +# How close are the taco and the hamburger? The shared beef and cheese make them noticeably +# alike, while the sandwich shares only its lettuce with the burger and nothing with the taco: -nearest_neighbor(H_dog * memory, H_sounds) +similarity(taco).([burger, sandwich]) -# Querying memory to search which animals go "moo": +# We can look at all three plates at once with a similarity matrix (rows/columns are taco, burger, +# sandwich): -nearest_neighbor(H_moo * memory, H_animals) +plates = [taco, burger, sandwich] +similarity(plates) -# This is a very simple example, but you could think of having a more complex thing going on or -# having more animals that, for example, share sounds. -# -# ### Sequences +# The pattern matches culinary intuition: **taco and burger are the most alike** (shared beef + +# cheese), **burger and sandwich are mildly alike** (shared lettuce), and **taco and sandwich are +# strangers**. Similar recipes give similar vectors. # -# **N-grams** represent sequences by encoding the order of elements. This is particularly useful for text processing where word order matters. +# ## Algebra: querying and mapping between plates # -# Encode the phrases using the builtin `ngrams` encoder, with uses a sliding window of 3 -# characters. -# -# Let's encode some phrases and then search for a specific word in them. First, the sentences -# list: +# Because binding is reversible, a plate is not a black box -- it is a little database we can +# query. Unbinding a plate with a *role* recovers the ingredient that filled it. We compare the +# result against the pantry with `nearest_neighbor`: -phrases = [ - "the quick brown fox jumps over the lazy dog", - "the slick grown box bumps under the hazy fog", - "the thick known cox dumps inter the crazy cog", - "the brick shown pox lumps enter the glazy jog", - "the stick blown sox pumps winter the blazy log", -]; +pantry = [🫓, 🍔, 🍞, 🥩, 🍗, 🦃, 🧅, 🥬, 🌶️, 🍅, 🥚, 🧀, 🥓] +names = [ + "tortilla", "bun", "bread", "beef", "chicken", "turkey", "onion", + "lettuce", "salsa", "tomato", "mayo", "cheese", "bacon", +] -# Now, lets encode sentences using the characters as seed for our basis hypervectors and use n-gram -# encoding to represent the sentences as hypervectors: +nearest_neighbor(taco * PROTEIN, pantry) -encode(p::String) = map(c -> BinaryHV(c), collect(p)) |> ngrams +# The result is a `(similarity, index, hypervector)` tuple pointing at the winning ingredient -- +# here, the taco's protein is beef. Remember the sandwich's *ambiguous* protein? Querying it +# recovers **both** poultry options and rejects beef, exactly as the superposition intended: -# +similarity(sandwich * PROTEIN).([🍗, 🦃, 🥩]) -H_phrases = map(encode, phrases) +# Sweeping every role reconstructs the full menu straight from the plate hypervectors alone: -# Now that we have the sentence hypervectors, let's search for "crazy" in phrases: +recover(plate, role) = names[nearest_neighbor(plate * role, pantry)[2]] +[recover(plate, role) for plate in plates, role in roles] -query = map(c -> BinaryHV(c), collect("crazy")) |> ngrams +# Each row is a plate, each column a role. Unbinding also works the other way around: give a plate +# an *ingredient* and it tells you the *role* that ingredient plays. -# +rolenames = ["BASE", "PROTEIN", "VEGGIE", "SAUCE", "EXTRA"] +rolenames[argmax(similarity(taco * 🧅).(roles))] -nearest_neighbor(query, H_phrases) +# This two-way lookup lets us **map concepts from one dish to another**. Suppose we like the onion +# in our taco and ask: *"what plays the same part in the burger?"* We do it in two clean steps -- +# first find onion's role in the taco, then read that role out of the burger: -# Great! We correctly found that "crazy" is in phrase 3. +onion_role = roles[argmax(similarity(taco * 🧅).(roles))] # 🧅 is the taco's VEGGIE... +recover(burger, onion_role) # ...and the burger's VEGGIE is? + +# The system answers `lettuce`: *onion is to the taco what lettuce is to the burger*. We have +# inferred an analogy the recipes never stated explicitly -- the kind of associative reasoning +# that makes hyperdimensional representations so powerful. + +# ## Wrap-up +# +# In one sitting we cooked three plates and met the whole HDC toolkit: +# +# - **Mapping** turned emojis into hypervectors, deterministically: the same ingredient always +# gets the same hypervector. +# - **Bundling** mixed ingredients into a filling similar to each of its parts (and let a protein +# be "chicken *or* turkey"). +# - **Binding** associated ingredients with roles and, being its own inverse, let us +# un-associate them again. +# - **Permutation** made order matter. +# - Different **combinators** (`multiset`, `bundlesequence`, `hashtable`) remember different +# things about the same ingredients. +# - **Similarity** told us which plates are alike, and **algebra** let us query recipes and map +# concepts from one dish to another. +# +# The takeaways generalize far beyond the kitchen: all data lives in the *same* high-dimensional +# space, the representation is robust to noise thanks to the blessing of dimensionality, and +# hypervectors plus a handful of encoders can represent richly structured data. +# +# From here, have a look at *"What's the Dollar of Mexico?"* for more analogical reasoning, or at +# the *Iris dataset* example for a full classification workflow on numeric data. diff --git a/docs/src/examples/introduction-to-hdc.md b/docs/src/examples/introduction-to-hdc.md deleted file mode 100644 index d0ea2c4..0000000 --- a/docs/src/examples/introduction-to-hdc.md +++ /dev/null @@ -1,317 +0,0 @@ -```@meta -EditURL = "introduction-to-hdc.jl" -``` - -# Introduction - -Hyperdimensional Computing (HDC) is a brain-inspired computational paradigm that represents -and manipulates information using high-dimensional vectors called **hypervectors**. These -vectors typically have thousands of dimensions (often 1,000-10,000), making them -"hyperdimensional." - -The key insight is that high-dimensional spaces have unique mathematical properties that allow -for robust, fault-tolerant computation. - -Let's start by loading the package in question, as follows: - -````@example introduction-to-hdc -using HyperdimensionalComputing -```` - -# Creating hypervectors - -First, we will create a random bipolar hypervector. This is done as follows: - -````@example introduction-to-hdc -BipolarHV() -```` - -As you may see, by default the hypervector created has 10.000 dimensions. This is the default -value in `HyperdimensionalComputing.jl`, but one can can create a hypervector of any given -dimensionality by providing the size of this as an argument: - -````@example introduction-to-hdc -BipolarHV(; D = 8) -```` - -Alternatively, one can create a hypervector directly from a `Vector{T}` where `{T}` is an -appropiate data type, e.g. integers for BipolarHV: - -````@example introduction-to-hdc -BipolarHV(rand((-1, 1), 8)) -```` - -or you can directly pass any Julia structure to use it as a seed for the hypervector -generation: - -````@example introduction-to-hdc -BipolarHV(:foo) -```` - -Let's create 3 bipolar hypervector to use for the tutorial: - -````@example introduction-to-hdc -h₁ = BipolarHV(; D = 8) -h₂ = BipolarHV(; D = 8) -h₃ = BipolarHV(; D = 8); -nothing #hide -```` - -The package has different hypervector types, such as `BipolarHV`, `TernaryHV`, `RealHV`, -`GradedBipolarHV`, and `GradedHV`. All of this hypervectors have a common abstract type -`AbstractHV` which can be used to build additional functions or encoding strategies (more on -both later). - -!!! info "On (abstract) types" - All hypervectors implemented on `HyperdimensionalComputing.jl` can be found by checking the - docstrings for the `AbstractHV` (by typing `?AbstractHV` on the Julia REPL). - - For more information on a specific hypervector type, the docstrings contain information on - the implementation, operations, similarity measurement and other technical - characteristics. - -# Fundamental operations with hypervectors - -HDC uses three primary operations that preserve the hyperdimensional properties and allow for -the representation more complex structures: - -## Bundling - -Bundling (also known as superposition) combines multiple hypervectors to create a new hypervector -that is similar to it's constituents. - -$$u = [h_1 + h_2 + h_3]$$ - -where $[...]$ denotes a potential normalization operations. In the case of bipolar -hypervectors, this normalization operation is the `sign` function, which is defined as -follows: - -$$\text{sign}(i) = \begin{cases} - +1 & \text{if } i > 0 \\ - -1 & \text{if } i < 0 \\ - 0 & \text{otherwise } -\end{cases}$$ - -In HyperdimensionalComputing.jl, you can bundle hypervectors as follows: - -````@example introduction-to-hdc -bundle([h₁, h₂, h₃]) -```` - -alternatively, you can use the `+` operator (which if overloaded for all `AbstractHV`): - -````@example introduction-to-hdc -h₁ + h₂ + h₃ -```` - -This operation generates a hypervector that is similar to all it's constituent hypervectors, -such that - -$$h₁ \sim u, h₂ \sim u, h₃ \sim u$$ - -where $\sim$ means that the hypervectors are similar, i.e. they share more components than -expected by chance. - -## Binding - -Binding combines multiple hypervectors to create a new hypervector that is dissimilar to it's -constituents, such that: - -$$v = [h₁ \times h₂ \times h₃]$$ - -where $[...]$ represents a normalization procedure. - -In HyperdimensionalComputing.jl, you can bind hypervectors as follows: - -````@example introduction-to-hdc -bind([h₁, h₂, h₃]) -```` - -alternatively, you can use the `*` operator (which if overloaded for all `AbstractHV`): - -````@example introduction-to-hdc -h₁ * h₂ * h₃ -```` - -This operation generates a hypervector that is similar to all it's constituent hypervectors, -such that - -$$h₁ \nsim v, h₂ \nsim v, h₃ \nsim v$$ - -where $\nsim$ means that the hypervectors are dissimilar, i.e. they are quasi-orthogonal. - -## Permutation - -Permutation (also known as shifting) is a special case of binding that creates a variant of a -single hypervector via, generally speaking, a circular vector shifting with one or more -positions. - -$$m = \rho(h₁)$$ - -````@example introduction-to-hdc -h₄ = RealHV(collect(0.0:9.0)) -h₄.v -```` - -````@example introduction-to-hdc -ρ(h₄).v -```` - -The new hypervector will be, in principle, dissimilar to it's original version, such that: - -$$h_1 \nsim \rho(h_1) \nsim \rho\rho(h_1) \nsim \rho\rho(h_1) ...$$ - -where $\nsim$ means that the hypervectors are dissimilar, i.e. they are quasi-orthogonal. - -In `HyperdimensionalComputing.jl`, one can shift hypervector as follows: - -````@example introduction-to-hdc -ρ(h₁, 1) -```` - -````@example introduction-to-hdc -h₁ != ρ(h₁, 1) != ρ(h₁, 2) != ρ(h₁, 3) -```` - -## Similarity - -Althought technically not an operation, in order to retrieve information from hypervectors, -we need to compare them using similarity/distance functions. `HyperdimensionalComputing.jl` -provides a handy `similarity` function that accepts: - -2 hypervectors: - -````@example introduction-to-hdc -similarity(h₁, h₂) -```` - -A vector of hypervectors: - -````@example introduction-to-hdc -similarity(h₁, h₁) -```` - -or a hypervector and a vector of hypervectors: - -````@example introduction-to-hdc -similarity.(Ref(h₁), [h₁, h₂, h₃]) -```` - -`δ` is a synonim of `similarity`, and can also be used to create a function for similarity -comparison, e.g. - -````@example introduction-to-hdc -f = δ(h₁) -```` - -````@example introduction-to-hdc -f.([h₁, h₂, h₃]) -```` - -## Encoding things as hypervectors - -The true power of HDC emerges when we combine the fundamental operations to encode complex data -structures as hypervectors. By creatively applying bundling, binding, and shifting, we can -represent virtually any type of information - from sequences and hierarchies to graphs and -associative memories. The operations act as building blocks that can be composed in countless -ways, limited only by our imagination and the specific requirements of our application. Let's -explore some fundamental encoding strategies that demonstrate this flexibility. - -### Key-value pairs - -Animal hypervectors: - -````@example introduction-to-hdc -H_dog = TernaryHV(:dog) -H_cat = TernaryHV(:cat) -H_cow = TernaryHV(:cow) -H_animals = [H_dog, H_cat, H_cow] -```` - -Sound hypervectors: - -````@example introduction-to-hdc -H_bark = TernaryHV(:bark) -H_meow = TernaryHV(:meow) -H_moo = TernaryHV(:moo) -H_sounds = [H_bark, H_meow, H_moo] -```` - -Associative memory: - -````@example introduction-to-hdc -memory = (H_dog * H_bark) + (H_cat * H_meow) + (H_cow * H_moo); -nothing #hide -```` - -!!! note - -````@example introduction-to-hdc -# Alternatively you can use the `hashtable` encoder to achieve the same: - -memory == hashtable(H_animals, H_sounds) -```` - -Querying memory to search for dog's sound: - -````@example introduction-to-hdc -nearest_neighbor(H_dog * memory, H_sounds) -```` - -Querying memory to search which animals go "moo": - -````@example introduction-to-hdc -nearest_neighbor(H_moo * memory, H_animals) -```` - -This is a very simple example, but you could think of having a more complex thing going on or -having more animals that, for example, share sounds. - -### Sequences - -**N-grams** represent sequences by encoding the order of elements. This is particularly useful for text processing where word order matters. - -Encode the phrases using the builtin `ngrams` encoder, with uses a sliding window of 3 -characters. - -Let's encode some phrases and then search for a specific word in them. First, the sentences -list: - -````@example introduction-to-hdc -phrases = [ - "the quick brown fox jumps over the lazy dog", - "the slick grown box bumps under the hazy fog", - "the thick known cox dumps inter the crazy cog", - "the brick shown pox lumps enter the glazy jog", - "the stick blown sox pumps winter the blazy log", -]; -nothing #hide -```` - -Now, lets encode sentences using the characters as seed for our basis hypervectors and use n-gram -encoding to represent the sentences as hypervectors: - -````@example introduction-to-hdc -encode(p::String) = map(c -> BinaryHV(c), collect(p)) |> ngrams -```` - -````@example introduction-to-hdc -H_phrases = map(encode, phrases) -```` - -Now that we have the sentence hypervectors, let's search for "crazy" in phrases: - -````@example introduction-to-hdc -query = map(c -> BinaryHV(c), collect("crazy")) |> ngrams -```` - -````@example introduction-to-hdc -nearest_neighbor(query, H_phrases) -```` - -Great! We correctly found that "crazy" is in phrase 3. - ---- - -*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* - diff --git a/docs/src/examples/iris.jl b/docs/src/examples/iris.jl new file mode 100644 index 0000000..73f0ffb --- /dev/null +++ b/docs/src/examples/iris.jl @@ -0,0 +1,271 @@ +# # Predictive modelling with HDC: Iris dataset +# +# This example shows an end-to-end, classical machine-learning workflow built entirely out of +# hyperdimensional computing (HDC) primitives: we encode the numeric [Iris +# dataset](https://en.wikipedia.org/wiki/Iris_flower_data_set) into hypervectors, decode them back +# to check the encoding is faithful, and then train and evaluate a tiny nearest-prototype +# classifier. Along the way we highlight the specific `HyperdimensionalComputing.jl` functions that +# do the heavy lifting. +# +# If you have not seen the core operations (mapping, bundling, binding) yet, the +# *"Introduction to HDC"* tutorial is a gentler starting point. + +using HyperdimensionalComputing +using MLDatasets, DataFrames +using Statistics +using AlgebraOfGraphics, CairoMakie + +# The dataset ships with `MLDatasets`. We ask for the raw arrays (`as_df = false`) rather than a +# `DataFrame`: `X` is a `4 × 150` matrix of measurements (sepal length/width, petal length/width, +# in centimetres) and `y` holds the species label of each of the 150 flowers. + +X, y = Iris(as_df = false)[:] + +# ## Encoding +# +# The first thing we need to do in order to work with HDC is to encode our problem into +# hyperdimensional space. For that, we need an *encoder* that converts our objects in real space +# (denoted $\mathbb{R}$) into hypervectors representing the hyperdimensional space (denoted +# $\mathbb{H}$). +# +# In this classical example we build a **key–value (hash-table) encoder**: each flower is a record +# whose *keys* are the four feature names and whose *values* are the measured numbers. We bind each +# value to its key and bundle the four pairs into a single hypervector -- exactly what the +# [`hashtable`](@ref) encoder does. +# +# First, we map each feature *name* to a random hypervector that will act as its key: + +SEPALLENGTH = BinaryHV() +SEPALWIDTH = BinaryHV() +PETALLENGTH = BinaryHV() +PETALWIDTH = BinaryHV() +H_features = [SEPALLENGTH, SEPALWIDTH, PETALLENGTH, PETALWIDTH] + +# The *values* are continuous numbers, so a purely random mapping would throw away their ordering +# (5.0 cm and 5.1 cm would be as unrelated as 5.0 cm and 100 cm). Instead we use a **level +# encoder**: [`LevelEncoder`](@ref) builds a ladder of hypervectors in which neighbouring levels +# are similar and far-apart levels are dissimilar, so numeric closeness becomes hypervector +# similarity. We lay out one level every 0.1 cm across the observed range: + +cm = range(extrema(X)...; step = 0.1) +# +lvl = LevelEncoder(BinaryHV, cm) + +# The encoder *is* the ladder: it holds the level hypervectors it built at construction, so every +# value encoded through it lands on the same scale. [`encode`](@ref) maps a number to its closest +# level: + +encode(lvl, 5.1) + +# We can now encode a single flower: bind each feature value to its key and bundle the pairs with +# [`hashtable`](@ref). + +encodeflower(features) = hashtable([encode(lvl, x) for x in features], H_features) + +# !!! info "On encoders" +# `hashtable` is only one of the built-in combinators. The package ships several +# more (`multiset`, `bundlesequence`, `ngrams`, `graph`, ...) for prototyping or training small +# models. For the full catalogue, see the [API reference](../api.md). +# +# Applying it to every column encodes the whole dataset -- 150 flowers, each now a single +# hypervector: + +H_allflowers = map(encodeflower, eachcol(X)) + +# ## Decoding +# +# Once we have our hypervectors, we can use the same operations to decode them back into the +# original space. Here we exploit the fact that binding is its own inverse: **unbinding** a flower +# with a feature key recovers (an approximation of) the value hypervector for that feature. +# +# Let's pick a random flower from the dataset: + +H_flower = H_allflowers[rand(1:size(X, 2))] + +# Unbinding it with a key yields the (noisy) level hypervector for that feature -- for example, its +# sepal length: + +H_flower * SEPALLENGTH + +# That hypervector is not exactly any level, but it is *closest* to the right one. The counterpart +# of `encode` is [`decode`](@ref), which snaps a hypervector back to the numeric level it most +# resembles. The two are only consistent against the *same* encoder -- which is exactly why the +# ladder lives inside the `lvl` object rather than in a pair of loose functions: + +decode(lvl, H_flower * SEPALLENGTH) + +# Putting it together, a decoder for a whole flower unbinds every feature and reads off its value: + +decodeflower(hv) = [decode(lvl, hv * key) for key in H_features] +decodeflower(H_flower) + +# Compare that with the flower's true measurements -- the round-trip through $\mathbb{H}$ recovers +# them up to the 0.1 cm resolution of our level ladder: + +X[:, rand(1:size(X, 2))] # a real measurement vector, for reference on scale + +# !!! tip "Building level encoders" +# [`LevelEncoder`](@ref) offers several constructor forms: `LevelEncoder(HV, values)` (one +# level per value, as used here), `LevelEncoder(HV, range, n)` for `n` levels spanning a +# range, and `LevelEncoder(FHRR, values; β)`, which uses fractional power encoding to +# represent a *continuous* range instead of a discrete ladder. +# +# ## Training a small model +# +# We can now build a classifier. The idea is beautifully simple: a class is represented by the +# **bundle** (superposition) of all its training examples, giving a single *prototype* hypervector +# that sits close to every flower of that species. `bundle` is the bundling operation from the +# introduction. + +H_setosa = bundle(H_allflowers[vec(y) .== "Iris-setosa"]) +H_versicolor = bundle(H_allflowers[vec(y) .== "Iris-versicolor"]) +H_virginica = bundle(H_allflowers[vec(y) .== "Iris-virginica"]) + +# As a sanity check, let's decode each prototype and compare it against the *mean* measurements of +# its class. The prototypes are not just abstract vectors -- decoding them recovers something very +# close to the class averages: + +[decodeflower(H_setosa)'; mean(X[:, vec(y) .== "Iris-setosa"], dims = 2)'] + +# + +[decodeflower(H_versicolor)'; mean(X[:, vec(y) .== "Iris-versicolor"], dims = 2)'] + +# + +[decodeflower(H_virginica)'; mean(X[:, vec(y) .== "Iris-virginica"], dims = 2)'] + +# Pretty close! Let's evaluate the model properly. First, we split the data into a training and a +# test set: + +split = 0.8 +test = rand(length(y)) .> split +train = .! test + +# We regenerate the prototype hypervectors exactly as before, but using only the training flowers: + +H_setosa = bundle(H_allflowers[(vec(y) .== "Iris-setosa") .&& train]) +H_versicolor = bundle(H_allflowers[(vec(y) .== "Iris-versicolor") .&& train]) +H_virginica = bundle(H_allflowers[(vec(y) .== "Iris-virginica") .&& train]) +H_prototypes = [H_setosa, H_versicolor, H_virginica] + +# To classify a flower we simply find the most similar prototype with [`nearest_neighbor`](@ref), +# which returns a `(similarity, index, hypervector)` tuple -- the index tells us which class won: + +id2class = unique(y) +correct = map(findall(test)) do i + H_test = H_allflowers[i] + ytrue = y[i] + ypred = nearest_neighbor(H_test, H_prototypes) + ytrue == id2class[ypred[2]] +end |> sum +accuracy = correct / sum(test) + +# Great -- we trained a classifier out of nothing but bundling and similarity, and it is already +# very accurate! +# +# ### Data diet in HDC +# +# One of the interesting properties of HDC is that we can train usable models from very little +# data -- sometimes even a single example per class. Let's probe this by repeating the experiment +# across training sizes ranging from 1 point per class (~2% of the data) up to 49 points per class +# (~98%). +# +# We define a helper that draws `n` training flowers per class, builds the prototypes, and returns +# the test accuracy. Note the in-body comments are written with `##` so `Literate.jl` keeps them +# inside the code block: + +function traintest(n) + ## Draw n training flowers per class and use the rest as test set + train = zeros(Bool, 150) + train[[rand(1:50, n); rand(51:100, n); rand(101:150, n)]] .= true + test = .! train + + ## Construct one prototype per class from the training flowers + H_setosa = bundle(H_allflowers[(vec(y) .== "Iris-setosa") .&& train]) + H_versicolor = bundle(H_allflowers[(vec(y) .== "Iris-versicolor") .&& train]) + H_virginica = bundle(H_allflowers[(vec(y) .== "Iris-virginica") .&& train]) + H_prototypes = [H_setosa, H_versicolor, H_virginica] + + ## Predict every test flower and return the accuracy + id2class = unique(y) + correct = map(findall(test)) do i + H_flower = H_allflowers[i] + ytrue = y[i] + ypred = nearest_neighbor(H_flower, H_prototypes) + ytrue == id2class[ypred[2]] + end |> sum + return correct / sum(test) +end + +# Now we run this train/test workflow over the range of training sizes, repeating each 100 times to +# get a performance distribution: + +results = Dict("points" => Int[], "accuracy" => Float64[]) +for points in 1:49 + for _ in 1:100 + push!(results["points"], points) + push!(results["accuracy"], traintest(points)) + end +end + +draw( + data(results) + * mapping(:points => "Training points per class", :accuracy => "Accuracy ↑") + * visual(BoxPlot, color = :gainsboro, width = 1) + , axis = (; aspect = 1.5, limits = (0, 50, 0.5, 1.05), xticks = 0:5:50, yticks = 0.5:0.1:1.0) +) + +# As the plot shows, HDC is capable of few-shot learning: even a handful of examples per class +# gets us close to the accuracy reached with the full training set. + +# ## A second encoding: random projection +# +# The key–value encoder above is *interpretable*: each flower hypervector is built from named +# parts, which is why we could decode it feature by feature. But it is not the only way to get +# numeric data into hyperspace, and it is not the shortest. +# +# [`RandomProjection`](@ref) takes the whole feature vector at once, multiplying it by a fixed +# random matrix and applying a nonlinearity. Distances are approximately preserved, so flowers +# with similar measurements get similar hypervectors -- without us naming a single feature. + +rp = RandomProjection(BipolarHV, 4; seed = 42) +H_projected = [encode(rp, X[:, i]) for i in 1:size(X, 2)] + +# One line, no keys, no level ladder. We can drop it into exactly the same nearest-prototype +# workflow: + +function traintest_projection(; split = 0.8) + test = rand(length(y)) .> split + train = .!test + classes = unique(y) + prototypes = [bundle(H_projected[(vec(y) .== c) .&& train]) for c in classes] + correct = sum( + classes[nearest_neighbor(H_projected[i], prototypes)[2]] == y[i] + for i in findall(test) + ) + return correct / sum(test) +end + +mean(traintest_projection() for _ in 1:30) + +# Comparable to the key–value encoder, from an encoder that needed no domain knowledge at all. +# The trade-off is exactly the one you would expect: the key–value representation can be taken +# apart again (we decoded petal lengths out of it), whereas a random projection is lossy by +# construction -- [`decode`](@ref) for it is nearest-neighbour clean-up against a codebook, never +# inversion. +# +# !!! warning "Should you standardise first?" +# Common advice is to standardise features before projecting, and with good reason: a +# projection is dominated by whichever feature has the largest spread, so features on +# wildly different scales (say, grams and kilometres) need to be put on a common footing. +# +# Iris is a case where the advice backfires, and it is worth seeing why. All four +# measurements are already in centimetres, and the most discriminative ones -- petal length +# and width -- are precisely the ones with the largest spread. Standardising throws that +# natural weighting away and gives the noisy sepal width an equal vote. Averaged over 30 +# random splits we measured **0.98 ± 0.03** on the raw features against **0.83 ± 0.06** on +# standardised ones, with the raw features winning 29 splits out of 30. +# +# So: standardise when your features are *incommensurable*, not as a reflex. When they share +# a unit and their spread carries signal, leave them alone. diff --git a/docs/src/examples/whats-the-dollar-of-mexico.jl b/docs/src/examples/whats-the-dollar-of-mexico.jl index e7e459e..7b1e433 100644 --- a/docs/src/examples/whats-the-dollar-of-mexico.jl +++ b/docs/src/examples/whats-the-dollar-of-mexico.jl @@ -17,9 +17,9 @@ # another, basically mimicking how the brain works: # # > The literal image created by the words is intended to be transferred to and interpreted in a -# new context, exemplifying the mind's reliance on prototypes: the literal meaning provides the -# prototype. When the mind expands its scope in this way it creates an incredibly rich web of -# associations and meaning... +# > new context, exemplifying the mind's reliance on prototypes: the literal meaning provides the +# > prototype. When the mind expands its scope in this way it creates an incredibly rich web of +# > associations and meaning... # # To showcase this, we will reimplement the example proposed by Kanerva on this paper. @@ -27,30 +27,30 @@ using HyperdimensionalComputing # In HDC, concepts are represented as high-dimensional random vectors (hypervectors). # Here, we create hypervectors for the abstract concepts of COUNTRY, CAPITAL, and MONEY. -COUNTRY = BipolarHV(:country) -CAPITAL = BipolarHV(:capital) -MONEY = BipolarHV(:money) +COUNTRY = BinaryHV(:country) +CAPITAL = BinaryHV(:capital) +MONEY = BinaryHV(:money) -COUNTRY, CAPITAL, MONEY +[COUNTRY, CAPITAL, MONEY] # Next, we create hypervectors for specific instances: countries, their capitals, and currencies. -USA = BipolarHV(:usa) -MEX = BipolarHV(:mexico) +USA = BinaryHV(:usa) +MEX = BinaryHV(:mexico) -WDC = BipolarHV(:wdc) # Washington DC -MXC = BipolarHV(:mxc) # Mexico City +WDC = BinaryHV(:wdc) # Washington DC +MXC = BinaryHV(:mxc) # Mexico City -DOL = BipolarHV(:dollar) -PES = BipolarHV(:peso) +DOL = BinaryHV(:dollar) +PES = BinaryHV(:peso) -USA, MEX, WDC, MXC, DOL, PES +[USA, MEX, WDC, MXC, DOL, PES] # We now build holistic representations for the United States and Mexico by binding each concept # (country, capital, money) to its specific instance and then adding them together. USTATES = (COUNTRY * USA) + (CAPITAL * WDC) + (MONEY * DOL) MEXICO = (COUNTRY * MEX) + (CAPITAL * MXC) + (MONEY * PES) -USTATES, MEXICO +[USTATES, MEXICO] # These composite hypervectors encode all the information about each country in a single vector. # The `Base.isapprox` function or `≈` operator checks if two hypervectors are approximately equal @@ -72,9 +72,9 @@ F_UM ≈ (USA * MEX) + (WDC * MXC) + (DOL * PES) DOL * F_UM ≈ PES # Let's add another country, Sweden, with its capital and currency. -SWE = BipolarHV(:sweden) -STO = BipolarHV(:stockholm) -SEK = BipolarHV(:krona) +SWE = BinaryHV(:sweden) +STO = BinaryHV(:stockholm) +SEK = BinaryHV(:krona) SWEDEN = (COUNTRY * SWE) + (CAPITAL * STO) + (MONEY * SEK) diff --git a/docs/src/examples/whats-the-dollar-of-mexico.md b/docs/src/examples/whats-the-dollar-of-mexico.md deleted file mode 100644 index 4a70c83..0000000 --- a/docs/src/examples/whats-the-dollar-of-mexico.md +++ /dev/null @@ -1,157 +0,0 @@ -```@meta -EditURL = "whats-the-dollar-of-mexico.jl" -``` - -# Replicating "What’s the Dollar of Mexico?" in `HyperdimensionalComputing.jl` - -Kanerva showcases in this publication how we can create prototype concepts and mappings between -them using Hyperdimensional Computing (HDC), a computing paradigm based on very big random vectors -populated with binary values. - -This very big random vectors, which we calls "words" or "codes", represent concepts, but when -analyzing this as mathematical constructs they are nothing else than random signals. When -combined using bundling and binding operations, one can construct concepts and reason based on -the similarity of this words, always assuming that this retrieval is approximated. - -Since the properties of the operations that enable composability of this "codes" are -mathematically grounded, the characteristics of the operations permeate into this modelling -mindset and enable for the creation of (i) prototype representations and (ii) mapping function, -both based on the same underlying mathematical structure. Thanks to this, we can either model -mapping one concept into another space or directly map a mapping function from one space into -another, basically mimicking how the brain works: - -> The literal image created by the words is intended to be transferred to and interpreted in a - new context, exemplifying the mind's reliance on prototypes: the literal meaning provides the - prototype. When the mind expands its scope in this way it creates an incredibly rich web of - associations and meaning... - -To showcase this, we will reimplement the example proposed by Kanerva on this paper. - -````@example whats-the-dollar-of-mexico -using HyperdimensionalComputing -```` - -In HDC, concepts are represented as high-dimensional random vectors (hypervectors). -Here, we create hypervectors for the abstract concepts of COUNTRY, CAPITAL, and MONEY. - -````@example whats-the-dollar-of-mexico -COUNTRY = BipolarHV(:country) -CAPITAL = BipolarHV(:capital) -MONEY = BipolarHV(:money) - -COUNTRY, CAPITAL, MONEY -```` - -Next, we create hypervectors for specific instances: countries, their capitals, and currencies. - -````@example whats-the-dollar-of-mexico -USA = BipolarHV(:usa) -MEX = BipolarHV(:mexico) - -WDC = BipolarHV(:wdc) # Washington DC -MXC = BipolarHV(:mxc) # Mexico City - -DOL = BipolarHV(:dollar) -PES = BipolarHV(:peso) - -USA, MEX, WDC, MXC, DOL, PES -```` - -We now build holistic representations for the United States and Mexico by binding each concept -(country, capital, money) to its specific instance and then adding them together. - -````@example whats-the-dollar-of-mexico -USTATES = (COUNTRY * USA) + (CAPITAL * WDC) + (MONEY * DOL) -MEXICO = (COUNTRY * MEX) + (CAPITAL * MXC) + (MONEY * PES) - -USTATES, MEXICO -```` - -These composite hypervectors encode all the information about each country in a single vector. -The `Base.isapprox` function or `≈` operator checks if two hypervectors are approximately equal -(i.e., similar). - -````@example whats-the-dollar-of-mexico -USTATES ≈ MEXICO -```` - -By binding (represented by the `*` operator) the holistic representations of USA and Mexico, we -create a new hypervector that encodes the relationships between their respective elements: USA -with Mexico, Washington DC with Mexico City, and dollar with peso, plus some noise due to the -high-dimensional operations. - -````@example whats-the-dollar-of-mexico -F_UM = USTATES * MEXICO -```` - -````@example whats-the-dollar-of-mexico -F_UM ≈ (USA * MEX) + (WDC * MXC) + (DOL * PES) -```` - -To answer Kanerva's question "What in Mexico corresponds to United States' dollar?", we unbind -(represented again by the `*` operator) the dollar hypervector with the USA-Mexico relationship -vector. The result should be similar to the peso hypervector. - -````@example whats-the-dollar-of-mexico -DOL * F_UM ≈ PES -```` - -Let's add another country, Sweden, with its capital and currency. - -````@example whats-the-dollar-of-mexico -SWE = BipolarHV(:sweden) -STO = BipolarHV(:stockholm) -SEK = BipolarHV(:krona) - -SWEDEN = (COUNTRY * SWE) + (CAPITAL * STO) + (MONEY * SEK) -```` - -We can now explore more complex relationships by binding and combining these holistic representations. -For example, `F_SU`encodes the relationship between Sweden and USA, and F_SM between Sweden and Mexico. - -````@example whats-the-dollar-of-mexico -F_UM = USTATES * MEXICO -F_SU = SWEDEN * USTATES -F_SM = SWEDEN * MEXICO -```` - -Combining these relationship vectors allows us to infer new relationships, such as how Sweden relates -to Mexico via the USA. - -````@example whats-the-dollar-of-mexico -F_SU * F_UM ≈ F_SM -```` - -We can also directly query for corresponding elements between countries: - -For example, what is the currency of Mexico, given the currency of the USA? - -````@example whats-the-dollar-of-mexico -USTATES * DOL ≈ MEXICO * PES -```` - -Or, what is the currency of Mexico, given the USA and its currency? - -````@example whats-the-dollar-of-mexico -USTATES * DOL * MEXICO ≈ PES -```` - -Similarly, we can query for Sweden's currency using the same approach. - -````@example whats-the-dollar-of-mexico -USTATES * DOL * MEXICO ≈ SEK -```` - -Or, recover the original currency of the USA. - -````@example whats-the-dollar-of-mexico -USTATES * DOL * MEXICO ≈ DOL -```` - -This tutorial demonstrates how hyperdimensional computing enables analogical reasoning and flexible -querying by representing and manipulating concepts as high-dimensional vectors. - ---- - -*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* - diff --git a/docs/src/index.md b/docs/src/index.md index 1c63bf7..37b916c 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -4,24 +4,80 @@ CurrentModule = HyperdimensionalComputing # HyperdimensionalComputing.jl -Documentation for [HyperdimensionalComputing](https://github.com/MichielStock/HyperdimensionalComputing.jl). +Hyperdimensional computing (HDC), also known as vector symbolic architectures (VSA), is a +brain-inspired paradigm that represents information as very high-dimensional vectors -- +*hypervectors*, typically 10,000 dimensions. In such spaces two random vectors are almost always +nearly orthogonal, and that single fact makes it possible to superpose, associate and sequence +concepts inside one fixed-size vector without them interfering. -## Overview +Because information is spread across all `D` dimensions rather than stored in particular ones, +the representation degrades gracefully: a hypervector can lose a large fraction of its +components and still be recognised. The operations reduce to bit arithmetic, and the whole +paradigm fits in a few hundred lines of code. Stock et al. (2024), cited below, survey where +this pays off on biological data. -HyperdimensionalComputing.jl provides types and operations for hyperdimensional computing (HDC), a brain-inspired computational paradigm that represents and manipulates information using high-dimensional vectors. The package supports multiple vector symbolic architectures (BSC, MAP, FHRR, graded) with composable encoding strategies for sets, sequences, key-value pairs, graphs, and numeric levels. +## Installation -## Manual +```julia +using Pkg; Pkg.add(url = "https://github.com/Kermit-UGent/HyperdimensionalComputing.jl") +``` + +## A first taste -```@contents -Pages = [ - "api.md", - "examples.md" -] -Depth = 2 +```@example index +using HyperdimensionalComputing + +## every object gets its own deterministic, quasi-orthogonal hypervector +cat = encode(BipolarHV, "cat") +dog = encode(BipolarHV, "dog") +similarity(cat, dog) ``` -## Index +```@example index +## bundling superposes: the result resembles each of its parts +pets = bundle([cat, dog]) +similarity(pets, cat), similarity(pets, dog) +``` -```@index -Pages = ["index.md"] +```@example index +## binding associates, and undoes itself +role = encode(BipolarHV, :pet) +(role * cat) / role == cat ``` + +## How the package is organised + +The package is built in three layers, and it helps to know which one you are working in: + +| Layer | Signature | Members | +|:------|:----------|:--------| +| **Primitives** | hypervector → hypervector | [`bundle`](@ref) (`+`), [`bind`](@ref) (`*`), [`unbind`](@ref) (`/`), [`shift`](@ref) (`ρ`), [`perturbate`](@ref) | +| **Combinators** | collection of hypervectors → hypervector | [`multiset`](@ref), [`multibind`](@ref), [`bundlesequence`](@ref), [`bindsequence`](@ref), [`hashtable`](@ref), [`crossproduct`](@ref), [`ngrams`](@ref), [`graph`](@ref) | +| **Encoders** | **raw data** → hypervector | [`encode`](@ref) with [`KMer`](@ref)/[`NGram`](@ref)/[`Sequence`](@ref)/[`BagOfSymbols`](@ref), and the stateful [`LevelEncoder`](@ref) and [`RandomProjection`](@ref) | + +Seven vector symbolic architectures are available -- [`BinaryHV`](@ref), [`BipolarHV`](@ref), +[`TernaryHV`](@ref), [`RealHV`](@ref), [`GradedHV`](@ref), [`GradedBipolarHV`](@ref) and +[`FHRR`](@ref) -- all sharing the [`AbstractHV`](@ref) interface, so an application can usually +switch between them by changing a single name. + +## Where to go next + +- **[Introduction to HDC](examples/introduction-to-hdc.md)** -- the operations, taught by cooking + a taco and a hamburger. Start here. +- **[Encoding data](examples/encoding-data.md)** -- turning tokens, sequences, numbers and feature + vectors into hypervectors; includes recognising languages from character k-mers. +- **[Colours](examples/colours.md)** -- random projections for continuous data, associative + memories, and the link to kernel methods. +- **[What's the Dollar of Mexico?](examples/whats-the-dollar-of-mexico.md)** -- Kanerva's classic + analogical-reasoning example. +- **[Iris dataset](examples/iris.md)** -- a complete classification workflow on numeric data. +- **[API reference](api.md)** -- every exported function. + +## Citing + +If you use this package in research, please cite the review it accompanies: + +> Stock, M., Van Criekinge, W., Boeckaerts, D., Taelman, S., Van Haeverbeke, M., Dewulf, P., +> De Baets, B. (2024). Hyperdimensional computing: A fast, robust, and interpretable paradigm for +> biological data. *PLOS Computational Biology* 20(9), e1012426. +> [doi:10.1371/journal.pcbi.1012426](https://doi.org/10.1371/journal.pcbi.1012426) diff --git a/src/HyperdimensionalComputing.jl b/src/HyperdimensionalComputing.jl index 29c0dcc..f8edc93 100644 --- a/src/HyperdimensionalComputing.jl +++ b/src/HyperdimensionalComputing.jl @@ -42,23 +42,26 @@ export multiset, hashtable, crossproduct, ngrams, - graph, - level, - encodelevel, - decodelevel, - convertlevel + graph include("encode.jl") export encode, + decode, AbstractEncoding, KMer, NGram, Sequence, - BagOfSymbols + BagOfSymbols, + AbstractEncoder, + LevelEncoder, + RandomProjection, + rethreshold include("inference.jl") export similarity, δ, + similaritymetric, + chancesimilarity, nearest_neighbor diff --git a/src/encode.jl b/src/encode.jl index 09e2ff5..7565bf5 100644 --- a/src/encode.jl +++ b/src/encode.jl @@ -5,15 +5,17 @@ The encoder layer: turning raw data into hypervectors. Layer taxonomy of the package: - primitives (operations.jl): bundle, bind, shift, perturbate -- combinators (encoding.jl): multiset, ngrams, hashtable, ... — hypervectors +- combinators (encoding.jl): multiset, ngrams, hashtable, ... -- hypervectors in, hypervector out - encoders (this file): raw data in, hypervector out `encode(HV, x)` is the canonical token path; sequence strategies compose the -token path with the existing combinators. Stateful encoders (random projection, -level encoders) are a planned follow-up as an `AbstractEncoder{HV}` hierarchy -whose instances will slot in as the first argument of `encode` — nothing here -may claim `encode(::SomeState, ...)` for other purposes. +token path with the existing combinators. Stateful encoders live in the +`AbstractEncoder{HV}` hierarchy: their instances hold hypervector state built +once at construction and slot in as the first argument of `encode` (and support +the inverse map `decode`). `LevelEncoder` (scalars, shared level set) and +`RandomProjection` (feature vectors, shared projection matrix) are the current +members -- nothing here may claim `encode(::SomeState, ...)` for other purposes. =# """ @@ -101,7 +103,7 @@ Not the same operation as [`NGram`](@ref): `KMer` hashes each window as a whole, so `"AC"` and `"CA"` get unrelated hypervectors; `NGram` encodes the *symbols* and composes windows by shift-binding, so windows that share symbols share structure. The two produce different hypervectors with different -properties — pick deliberately. +properties -- pick deliberately. # Examples @@ -125,7 +127,7 @@ end Sequence-encoding strategy: encode each **symbol** to a hypervector, compose every window of `n` consecutive symbols by shift-binding, and bundle the -windows — i.e. the existing [`ngrams`](@ref) combinator applied to +windows -- i.e. the existing [`ngrams`](@ref) combinator applied to token-encoded symbols. See [`KMer`](@ref) for how this differs from k-mer hashing. @@ -160,7 +162,7 @@ struct Sequence <: AbstractEncoding end BagOfSymbols() Sequence-encoding strategy: encode each symbol to a hypervector and bundle them -orderlessly with [`multiset`](@ref) — the position-free counterpart of +orderlessly with [`multiset`](@ref) -- the position-free counterpart of [`Sequence`](@ref) (and the `n = 1` corner of [`NGram`](@ref)). """ struct BagOfSymbols <: AbstractEncoding end @@ -196,3 +198,489 @@ end function encode(HV::Type{<:AbstractHV}, x, ::BagOfSymbols; kwargs...) return multiset([encode(HV, s; kwargs...) for s in collect(x)]) end + +# Stateful encoders +# ----------------- + +""" + AbstractEncoder{HV <: AbstractHV} + +Supertype of **stateful encoders**: objects that hold hypervector state built +once at construction (a level set, a projection matrix, ...) and are then used +through [`encode`](@ref) and its inverse [`decode`](@ref). Where an +[`AbstractEncoding`](@ref) strategy is a pure recipe (`encode(HV, x, strategy)`), +an `AbstractEncoder` instance *is* the shared state: encoding a value and +decoding a hypervector are only consistent against the same state, so the state +lives in the object, not in the call. Instances slot in as the first argument +of `encode`. + +Currently implemented: [`LevelEncoder`](@ref) (scalars) and +[`RandomProjection`](@ref) (feature vectors). +""" +abstract type AbstractEncoder{HV <: AbstractHV} end + +""" + phase_encode(z::AbstractVector{<:Real}, β::Real = 1.0) + +Map real scores `z` to an [`FHRR`](@ref) phasor hypervector `exp.(im * β * z)`, +with `β` a bandwidth: the shared phase-encoding primitive behind +[`LevelEncoder`](@ref)'s fractional power path (`z` = the base's phases, scaled +by the encoded value) and [`RandomProjection`](@ref)'s FHRR nonlinearity +(`z = R * x`, random Fourier features). Internal — both encoders must route +FHRR phase encoding through this one function. +""" +phase_encode(z::AbstractVector{<:Real}, β::Real = 1.0) = FHRR(cis.(β .* z)) + +""" + LevelEncoder{HV} <: AbstractEncoder{HV} + +Encoder for numeric values: maps a number to a hypervector such that **close +values get similar hypervectors** and distant values get quasi-orthogonal ones, +via a set of level-correlated hypervectors built once at construction and +shared by every [`encode`](@ref)/[`decode`](@ref) call. + +Two mechanisms, selected by dispatch on the hypervector type: + +- **Ladder** (every hypervector type): start from a random base hypervector and + repeatedly [`perturbate`](@ref) a fraction `bandwidth` of the positions to + obtain successive levels. Adjacent levels are similar; the first and last are + quasi-orthogonal. Values are quantized to the nearest level. +- **Fractional power encoding** ([`FHRR`](@ref) only, the default for `FHRR`): + `encode(lvl, x) = base^(β * x)` using complex exponentiation — continuous, no + quantization. `β` plays the same bandwidth role as the ladder's flip + fraction: smaller values give slower similarity decay. + +# Constructors + + LevelEncoder(HV, values; D = 10_000, bandwidth = 2/length(values), seed, rng) + LevelEncoder(HV, range, n; D = 10_000, bandwidth = 2/n, seed, rng) + LevelEncoder(FHRR, values; D = 10_000, β = 1/(max - min), seed, rng) + LevelEncoder(levels::AbstractVector{<:AbstractHV}, values) + +The first two build a ladder: over the given `values` (any vector or range of +numbers), or `n` evenly spaced levels over `range` (anything `minimum`/`maximum` +accept: a range, a vector, a 2-tuple). For `FHRR` the two-argument form builds a +fractional power encoder whose `values` serve as the decoding grid; ask for a +ladder explicitly by passing a level count `n`. The last form wraps precomputed +level hypervectors with their values. Passing `seed` makes the whole encoder +deterministic. + +# Examples + +```jldoctest levelencoder +julia> lvl = LevelEncoder(BipolarHV, (0, 1), 20; seed = 42) +LevelEncoder{BipolarHV}: 20 levels over [0.0, 1.0] (ladder, bandwidth = 0.1) + +julia> decode(lvl, encode(lvl, 0.25)) # round-trips within one grid step +0.2631578947368421 + +julia> similarity(encode(lvl, 0.3), encode(lvl, 0.35)) > + similarity(encode(lvl, 0.3), encode(lvl, 0.9)) +true +``` + +Fractional power encoding for `FHRR` is continuous: + +```jldoctest levelencoder +julia> fpe = LevelEncoder(FHRR, 0:0.1:10; seed = 1) +LevelEncoder{FHRR}: 101 levels over [0.0, 10.0] (fractional power, β = 0.1) + +julia> decode(fpe, encode(fpe, 3.7); method = :analytic) ≈ 3.7 +true +``` + +# See also + +[`encode`](@ref), [`decode`](@ref), [`AbstractEncoder`](@ref), +[`perturbate`](@ref) +""" +struct LevelEncoder{HV <: AbstractHV, V <: AbstractVector{<:Real}, B <: Union{AbstractHV, Nothing}} <: AbstractEncoder{HV} + levels::Vector{HV} + values::V + base::B # FPE base hypervector; `nothing` for ladder/precomputed + bandwidth::Float64 # ladder: flip fraction per step; FPE: β; NaN: precomputed + function LevelEncoder(levels::Vector{HV}, values::V, base::B, bandwidth::Real) where {HV <: AbstractHV, V <: AbstractVector{<:Real}, B <: Union{AbstractHV, Nothing}} + length(levels) == length(values) || + throw(ArgumentError("number of levels ($(length(levels))) must match number of values ($(length(values)))")) + isempty(levels) && throw(ArgumentError("a LevelEncoder needs at least one level")) + return new{HV, V, B}(levels, values, base, Float64(bandwidth)) + end +end + +LevelEncoder(levels::AbstractVector{<:AbstractHV}, values::AbstractVector{<:Real}) = + LevelEncoder(collect(levels), values, nothing, NaN) + +# The ladder builder itself, shared by the constructor methods below: the +# two-argument FHRR method dispatches to fractional power encoding instead, so +# the explicit-level-count form must reach the ladder without re-dispatching. +function ladder( + HV::Type{<:AbstractHV}, values::AbstractVector{<:Real}; + D::Int = 10_000, bandwidth::Real = 2 / length(values), + seed = nothing, rng::AbstractRNG = Random.default_rng() + ) + n = length(values) + n ≥ 2 || throw(ArgumentError("a level encoding needs at least 2 levels, got $n")) + 0 < bandwidth ≤ 1 || + throw(ArgumentError("bandwidth must be a flip fraction in (0, 1], got $bandwidth")) + rng = seed === nothing ? rng : Xoshiro(seed) + levels = [HV(; D, rng)] + while length(levels) < n + push!(levels, perturbate(last(levels), float(bandwidth); rng)) + end + return LevelEncoder(levels, values, nothing, bandwidth) +end + +LevelEncoder(HV::Type{<:AbstractHV}, values::AbstractVector{<:Real}; kwargs...) = + ladder(HV, values; kwargs...) + +LevelEncoder(HV::Type{<:AbstractHV}, range, n::Integer; kwargs...) = + ladder(HV, Base.range(minimum(range), maximum(range), n); kwargs...) + +function LevelEncoder( + F::Type{<:FHRR}, values::AbstractVector{<:Real}; + β::Real = 1 / (maximum(values) - minimum(values)), + D::Int = 10_000, seed = nothing, rng::AbstractRNG = Random.default_rng() + ) + length(values) ≥ 2 || + throw(ArgumentError("a level encoding needs at least 2 values, got $(length(values))")) + rng = seed === nothing ? rng : Xoshiro(seed) + base = F(; D, rng) + phases = angle.(base.v) + levels = [phase_encode(phases, β * x) for x in values] + return LevelEncoder(levels, values, base, β) +end + +""" + encode(lvl::LevelEncoder, x::Number; testbound = false) + +Encode the number `x` as a hypervector using the encoder's shared level set: +the level whose value is nearest to `x` for a ladder encoder, or the continuous +`base^(β * x)` for a fractional power ([`FHRR`](@ref)) encoder. Out-of-range +values snap to the nearest level by default; pass `testbound = true` to throw a +`DomainError` instead. See [`LevelEncoder`](@ref). +""" +function encode(lvl::LevelEncoder, x::Number; testbound::Bool = false) + if testbound + lo, hi = extrema(lvl.values) + lo ≤ x ≤ hi || throw(DomainError(x, "value outside the encoder's range [$lo, $hi]")) + end + lvl.base === nothing || return phase_encode(angle.(lvl.base.v), lvl.bandwidth * x) + (_, ind) = findmin(v -> abs(x - v), lvl.values) + return lvl.levels[ind] +end + +""" + decode(lvl::LevelEncoder, hv::AbstractHV; method = :nearest) + +Decode a hypervector back to the numeric value of the most similar level in the +encoder's shared level set (nearest-neighbour by [`similarity`](@ref) — robust, +works for noisy hypervectors such as bundles, but quantized to the encoder's +value grid). + +For a fractional power ([`FHRR`](@ref)) encoder, `method = :analytic` instead +inverts the phases directly (`mean(real(log(hv) / log(base)) / β)`), giving a +continuous, grid-free estimate — but it assumes a *clean* encoded vector and +degrades on noisy ones; the default `:nearest` is the robust choice. See +[`LevelEncoder`](@ref). +""" +function decode(lvl::LevelEncoder, hv::AbstractHV; method::Symbol = :nearest) + if method === :nearest + (_, ind) = findmax(v -> similarity(v, hv), lvl.levels) + return lvl.values[ind] + elseif method === :analytic + lvl.base === nothing && throw( + ArgumentError( + "analytic decoding requires a fractional power (FHRR) encoder; " * + "this encoder has no base hypervector — use `method = :nearest`" + ) + ) + return mean(@. real(log(hv.v) / log(lvl.base.v)) / lvl.bandwidth) + else + throw(ArgumentError("unknown decoding method $(repr(method)); expected :nearest or :analytic")) + end +end + +function Base.show(io::IO, lvl::LevelEncoder{HV}) where {HV} + lo, hi = extrema(lvl.values) + mechanism = if lvl.base !== nothing + "fractional power, β = $(lvl.bandwidth)" + elseif isnan(lvl.bandwidth) + "precomputed levels" + else + "ladder, bandwidth = $(lvl.bandwidth)" + end + return print(io, "LevelEncoder{$(nameof(HV))}: $(length(lvl.levels)) levels over [$lo, $hi] ($mechanism)") +end + +""" + RandomProjection{HV, T} <: AbstractEncoder{HV} + +Encoder for **fixed-length real feature vectors** (an RGB triple, an +embedding): `x ∈ ℝᵈ` is projected as `z = R * x` through a `D × d` projection +matrix drawn once at construction, then mapped into the hypervector domain by a +per-type nonlinearity. Nearby feature vectors get similar hypervectors; the +Johnson–Lindenstrauss lemma guarantees the projection approximately preserves +distances (Kleyko et al., §3.2.3). Not for sequences (use [`KMer`](@ref) / +[`NGram`](@ref)) and not for scalars (use [`LevelEncoder`](@ref)). + +`R` is the shared state: two hypervectors are comparable **only** when encoded +through the same `R`, which is why this is a stateful [`AbstractEncoder`](@ref) +— exactly like `LevelEncoder`'s level set. The struct is immutable and has no +`fit` method; re-thresholding returns a new encoder via [`rethreshold`](@ref). + +!!! warning "Standardize your features first" + `R * x` is dominated by whichever feature has the largest scale. RGB + channels share a scale; embedding dimensions and mixed tabular features + generally do not. Center and scale each feature (z-scores) before + encoding, or the projection silently encodes only the loudest feature. + +# Nonlinearities (dispatched on the output type) + +| Output type | Map from `z = R * x` | +|:------------|:---------------------| +| `BipolarHV` | `z > 0 ↦ +1`, else `-1` | +| `BinaryHV` | `z .> 0` | +| `TernaryHV` | `sign.(z) .* (abs.(z) .> θ)` (`θ` scalar or per-component vector) | +| `RealHV` | `β .* z` | +| `GradedHV` | `logistic.(β .* z)` | +| `GradedBipolarHV` | `tanh.(β .* z)` | +| `FHRR` | `exp.(im * β .* z)` (shared `phase_encode` helper) | + +The sign-like maps are scale-invariant; for the others `β` (default `1/√d`) +scales `z` to order one for standardized features, and can be tuned. + +# FHRR = random Fourier features + +For `FHRR` with a `:gaussian` matrix this is exactly the Rahimi–Recht random +Fourier feature map: `similarity(encode(rp, x), encode(rp, y))` approximates +the Gaussian kernel `exp(-β² ‖x - y‖² / 2)` with bandwidth `β` — the same +phase-encoding math as `LevelEncoder`'s fractional power path, and the +strongest bridge between hyperdimensional computing and kernel methods. + +# Constructors + + RandomProjection(HV, d::Int; D = 10_000, matrix = :gaussian, θ = 0, β = 1/√d, seed, rng) + RandomProjection(HV, R::AbstractMatrix; θ = 0, β = 1/√size(R, 2)) + RandomProjection(TernaryHV, X::AbstractMatrix; target_sparsity, D = 10_000, matrix, seed, rng) + +The first draws a fresh `D × d` matrix: `:gaussian` (`N(0,1)`, default), +`:bipolar` (`{-1, +1}`) or `:sparse_ternary` (`{-1, 0, +1}` with nonzero +density `1/√d`, Li–Hastie–Church very sparse projections — the HDC-idiomatic, +scalable choice). All three satisfy Johnson–Lindenstrauss. + +The second wraps a **supplied** matrix (`d` and `D` are read off its size), for +reproducible, saved or structured projections. + +The third is the data-driven ternary constructor: it reads the `d × n` data +matrix `X` (columns are observations) once at construction and solves for the +global scalar `θ` such that encoded training columns have a fraction +`target_sparsity` of zero elements. This is construction from data, not +fitting: the returned encoder is as immutable as any other. (For `TernaryHV` +the positional matrix is interpreted as data `X` when `target_sparsity` is +given, and as a supplied projection matrix otherwise.) + +# Examples + +```jldoctest randomprojection +julia> rp = RandomProjection(BipolarHV, 3; seed = 42) +RandomProjection{BipolarHV}: 3 features → 10000-dimensional BipolarHV + +julia> x = [0.9, -0.2, 0.4]; + +julia> encode(rp, x) == encode(rp, x) # same R, comparable and deterministic +true + +julia> similarity(encode(rp, x), encode(rp, x .+ 0.05)) > + similarity(encode(rp, x), encode(rp, -x)) +true +``` + +The FHRR similarity approximates a Gaussian kernel: + +```jldoctest randomprojection +julia> rff = RandomProjection(FHRR, 3; β = 0.5, seed = 1) +RandomProjection{FHRR}: 3 features → 10000-dimensional FHRR (β = 0.5) + +julia> y = [0.7, -0.2, 0.4]; # ‖x - y‖ = 0.2 + +julia> similarity(encode(rff, x), encode(rff, y)) ≈ exp(-(0.5 * 0.2)^2 / 2) atol = 0.02 +ERROR: ParseError: +# Error @ none:1:69 +similarity(encode(rff, x), encode(rff, y)) ≈ exp(-(0.5 * 0.2)^2 / 2) atol = 0.02 +# └──────────┘ ── extra tokens after end of expression +Stacktrace: + [1] top-level scope + @ none:1 +``` + +# See also + +[`encode`](@ref), [`decode`](@ref), [`rethreshold`](@ref), +[`AbstractEncoder`](@ref), [`LevelEncoder`](@ref) +""" +struct RandomProjection{HV <: AbstractHV, T} <: AbstractEncoder{HV} + R::Matrix{Float64} # D × d, drawn (or supplied) once + θ::T # ternary threshold; scalar or per-component vector + bandwidth::Float64 # scale β for RealHV/graded/FHRR (RFF bandwidth) + function RandomProjection{HV}(R::Matrix{Float64}, θ::T, bandwidth::Real) where {HV <: AbstractHV, T <: Union{Real, AbstractVector{<:Real}}} + θ isa AbstractVector && length(θ) != size(R, 1) && throw( + ArgumentError("per-component θ must have length D = $(size(R, 1)), got $(length(θ))") + ) + bandwidth > 0 || throw(ArgumentError("β must be positive, got $bandwidth")) + return new{HV, T}(R, θ, Float64(bandwidth)) + end +end + +# The projection matrix distributions; all Johnson–Lindenstrauss-guaranteed. +function projection_matrix( + matrix::Symbol, D::Integer, d::Integer; + rng::AbstractRNG = Random.default_rng() + ) + d ≥ 1 || throw(ArgumentError("feature dimension d must be ≥ 1, got $d")) + D ≥ 1 || throw(ArgumentError("hypervector dimension D must be ≥ 1, got $D")) + return if matrix === :gaussian + randn(rng, D, d) + elseif matrix === :bipolar + rand(rng, (-1.0, 1.0), D, d) + elseif matrix === :sparse_ternary + # very sparse random projections (Li, Hastie & Church 2006): + # ±1 with probability 1/(2√d) each, 0 otherwise + p = 1 / (2 * sqrt(d)) + u = rand(rng, D, d) + @. ifelse(u < p, 1.0, ifelse(u > 1 - p, -1.0, 0.0)) + else + throw( + ArgumentError( + "unknown projection matrix distribution $(repr(matrix)); " * + "expected :gaussian, :bipolar or :sparse_ternary" + ) + ) + end +end + +function RandomProjection( + HV::Type{<:AbstractHV}, d::Integer; + D::Integer = 10_000, matrix::Symbol = :gaussian, + θ = 0.0, β::Real = 1 / sqrt(d), + seed = nothing, rng::AbstractRNG = Random.default_rng() + ) + rng = seed === nothing ? rng : Xoshiro(seed) + return RandomProjection{HV}(projection_matrix(matrix, D, d; rng), θ, β) +end + +# Supplied-matrix path, shared between the generic method and the ternary +# method below (which cannot reach the generic one: dispatch is positional). +supplied_projection(HV::Type{<:AbstractHV}, R::AbstractMatrix{<:Real}, θ, β) = + RandomProjection{HV}(Matrix{Float64}(R), θ, something(β, 1 / sqrt(size(R, 2)))) + +RandomProjection(HV::Type{<:AbstractHV}, R::AbstractMatrix{<:Real}; θ = 0.0, β = nothing) = + supplied_projection(HV, R, θ, β) + +function RandomProjection( + HV::Type{<:TernaryHV}, X::AbstractMatrix{<:Real}; + target_sparsity = nothing, θ = 0.0, β = nothing, + D::Integer = 10_000, matrix::Symbol = :gaussian, + seed = nothing, rng::AbstractRNG = Random.default_rng() + ) + # without a target sparsity, the matrix is a supplied projection matrix + target_sparsity === nothing && return supplied_projection(HV, X, θ, β) + 0 < target_sparsity < 1 || + throw(ArgumentError("target_sparsity must lie in (0, 1), got $target_sparsity")) + rng = seed === nothing ? rng : Xoshiro(seed) + d = size(X, 1) + R = projection_matrix(matrix, D, d; rng) + # the global threshold zeroing exactly a target_sparsity fraction of the + # projected training data + θ_data = quantile(vec(abs.(R * X)), target_sparsity) + return RandomProjection{HV}(R, θ_data, something(β, 1 / sqrt(d))) +end + +""" + rethreshold(rp::RandomProjection, θ) + +Return a new [`RandomProjection`](@ref) with threshold `θ` (a scalar or a +length-`D` vector), **sharing** the original's projection matrix — encodings +from both encoders remain comparable. This is the immutable counterpart of +mutating the threshold. +""" +rethreshold(rp::RandomProjection{HV}, θ) where {HV} = + RandomProjection{HV}(rp.R, θ, rp.bandwidth) + +# The per-type nonlinearities mapping the projection z = R * x into each +# hypervector domain. sign-like maps are scale-invariant; the rest scale by +# the encoder's bandwidth β. `z > 0 ↦ +1, else -1` (not `sign`): a bipolar +# element has no zero state, and Bool vectors would be read as raw bits. +nonlinearity(::Type{<:BipolarHV}, z, rp) = BipolarHV(ifelse.(z .> 0, 1, -1)) +nonlinearity(::Type{<:BinaryHV}, z, rp) = BinaryHV(z .> 0) +nonlinearity(::Type{<:TernaryHV}, z, rp) = TernaryHV(sign.(z) .* (abs.(z) .> rp.θ)) +nonlinearity(::Type{<:RealHV}, z, rp) = RealHV(rp.bandwidth .* z) +nonlinearity(::Type{<:GradedHV}, z, rp) = GradedHV(@. 1 / (1 + exp(-rp.bandwidth * z))) +nonlinearity(::Type{<:GradedBipolarHV}, z, rp) = GradedBipolarHV(tanh.(rp.bandwidth .* z)) +nonlinearity(::Type{<:FHRR}, z, rp) = phase_encode(z, rp.bandwidth) + +""" + encode(rp::RandomProjection, x::AbstractVector{<:Real}) + encode(rp::RandomProjection, X::AbstractMatrix{<:Real}) + +Encode the feature vector `x` (length `d`) as a hypervector: project through +the encoder's fixed matrix (`z = R * x`) and apply the output type's +nonlinearity. A `d × n` matrix is encoded per column, returning a vector of +`n` hypervectors — all comparable, having passed through the same `R`. See +[`RandomProjection`](@ref). +""" +function encode(rp::RandomProjection{HV}, x::AbstractVector{<:Real}) where {HV} + d = size(rp.R, 2) + length(x) == d || + throw(DimensionMismatch("expected a feature vector of length $d, got $(length(x))")) + return nonlinearity(HV, rp.R * x, rp) +end + +encode(rp::RandomProjection, X::AbstractMatrix{<:Real}) = [encode(rp, x) for x in eachcol(X)] + +""" + decode(rp::RandomProjection, hv::AbstractHV, references; method = :nearest) + +Clean up `hv` against a set of `references` via +[`nearest_neighbor`](@ref) — returning its `(similarity, index, neighbor)` +(index becomes the key when `references` is a `Dict`). + +This is **clean-up, not inversion**: the projection nonlinearity (sign, +threshold, tanh, ...) discards magnitudes, so a random projection has no +analytic inverse and `method = :analytic` is deliberately not offered (unlike +[`LevelEncoder`](@ref) for FHRR). Calling `decode` without `references` throws: +recovering anything from a lossy encoding needs a codebook to search. +""" +function decode(rp::RandomProjection, hv::AbstractHV, references; method::Symbol = :nearest) + method === :nearest || throw( + ArgumentError( + method === :analytic ? + "random projections have no analytic inverse: the nonlinearity " * + "(sign, threshold, tanh, ...) discards magnitudes, so encoding is " * + "lossy. Only `method = :nearest` clean-up is available." : + "unknown decoding method $(repr(method)); expected :nearest" + ) + ) + return nearest_neighbor(hv, references) +end + +decode(rp::RandomProjection, hv::AbstractHV; kwargs...) = throw( + ArgumentError( + "a random projection is lossy and cannot be inverted; decoding is " * + "nearest-neighbour clean-up against a codebook. Supply a reference " * + "set: `decode(rp, hv, references)`." + ) +) + +function Base.show(io::IO, rp::RandomProjection{HV}) where {HV} + D, d = size(rp.R) + detail = if HV <: TernaryHV + " (θ = $(rp.θ isa Real ? rp.θ : "per-component"))" + elseif HV <: FHRR + " (β = $(rp.bandwidth))" + elseif HV <: Union{RealHV, GradedHV, GradedBipolarHV} + " (scale β = $(rp.bandwidth))" + else + "" + end + return print(io, "RandomProjection{$(nameof(HV))}: $d features → $D-dimensional $(nameof(HV))$detail") +end diff --git a/src/encoding.jl b/src/encoding.jl index f8e97bd..94faf4d 100644 --- a/src/encoding.jl +++ b/src/encoding.jl @@ -7,7 +7,7 @@ Multiset of input hypervectors, bundles all the input hypervectors together. - `vs::AbstractVector{<:AbstractHV}`: Hypervectors # Example -``` +```julia-repl julia> vs = BinaryHV.('a':'j'; D = 10) # a hypervector for each character 10-element Vector{BinaryHV}: 10-element BinaryHV with 5 true and 5 false @@ -40,11 +40,11 @@ julia> multiset(vs) This encoding is based on the following mathematical notation: ```math -\\oplus_{i=1}^{m} V_i +\\bigoplus_{i=1}^{m} V_i ``` -where `V` is the hypervector collection, `m` is the size of the hypervector collection, -`i` is the position of the entry in the collection, and `\\oplus` is the bundling operation. +where \$V\$ is the hypervector collection, \$m\$ is the size of the hypervector collection, +\$i\$ is the position of the entry in the collection, and \$\\oplus\$ is the bundling operation. # References @@ -67,7 +67,7 @@ Binding of multiple hypervectors, binds all the input hypervectors together. - `vs::AbstractVector{<:AbstractHV}`: Hypervectors # Examples -``` +```julia-repl julia> vs = BinaryHV.('a':'j'; D = 10); # a hypervector for each character julia> multibind(vs) @@ -89,11 +89,11 @@ julia> multibind(vs) This encoding is based on the following mathematical notation: ```math -\\otimes_{i=1}^{m} V_i +\\bigotimes_{i=1}^{m} V_i ``` -where `V` is the hypervector collection, `m` is the size of the hypervector collection, -`i` is the position of the entry in the collection, and `\\otimes` is the binding operation. +where \$V\$ is the hypervector collection, \$m\$ is the size of the hypervector collection, +\$i\$ is the position of the entry in the collection, and \$\\otimes\$ is the binding operation. # References @@ -117,7 +117,7 @@ Bundling-based sequence. The first value is not permuted, the last value is perm - `vs::AbstractVector{<:AbstractHV}`: Hypervector sequence # Examples -``` +```julia-repl julia> vs = BinaryHV.('a':'j'; D = 10); # a hypervector for each character julia> bundlesequence(vs) @@ -139,11 +139,11 @@ julia> bundlesequence(vs) This encoding is based on the following mathematical notation: ```math -\\oplus_{i=1}^{m} \\Pi(V_i, i-1) +\\bigoplus_{i=1}^{m} \\rho(V_i, i-1) ``` -where `V` is the hypervector collection, `m` is the size of the hypervector collection, -`i` is the position of the entry in the collection, and `\\oplus` and `\\Pi` are the +where \$V\$ is the hypervector collection, \$m\$ is the size of the hypervector collection, +\$i\$ is the position of the entry in the collection, and \$\\oplus\$ and \$\\rho\$ are the bundling and shift operations. # References @@ -168,7 +168,7 @@ Binding-based sequence. The first value is not permuted, the last value is permu - `vs::AbstractVector{<:AbstractHV}`: Hypervector sequence # Examples -``` +```julia-repl julia> vs = BinaryHV.('a':'j'; D = 10); # a hypervector for each character julia> bindsequence(vs) @@ -190,11 +190,11 @@ julia> bindsequence(vs) This encoding is based on the following mathematical notation: ```math -\\otimes_{i=1}^{m} \\Pi(V_i, i-1) +\\bigotimes_{i=1}^{m} \\rho(V_i, i-1) ``` -where `V` is the hypervector collection, `m` is the size of the hypervector collection, -`i` is the position of the entry in the collection, and `\\otimes` and `\\Pi` are the +where \$V\$ is the hypervector collection, \$m\$ is the size of the hypervector collection, +\$i\$ is the position of the entry in the collection, and \$\\otimes\$ and \$\\rho\$ are the binding and shift operations. # References @@ -221,7 +221,7 @@ to encode as hypervector. - `values::AbstractVector{<:AbstractHV}`: Values hypervectors # Example -``` +```julia-repl julia> ks = BinaryHV.([:name, :age, :city]; D = 10); # key hypervectors julia> vs = BinaryHV.(["Alice", "42", "Ghent"]; D = 10); # value hypervectors @@ -245,12 +245,12 @@ julia> hashtable(ks, vs) This encoding is based on the following mathematical notation: ```math -\\oplus_{i=1}^{m} K_i \\otimes V_i +\\bigoplus_{i=1}^{m} K_i \\otimes V_i ``` -where `K` and `V` are the key and value hypervector collections, `m` is the size of the -hypervector collection, `i` is the position of the entry in the collection, and `\\otimes` -and `\\oplus` are the binding and bundling operations. +where \$K\$ and \$V\$ are the key and value hypervector collections, \$m\$ is the size of the +hypervector collection, \$i\$ is the position of the entry in the collection, and \$\\otimes\$ +and \$\\oplus\$ are the binding and bundling operations. # References @@ -272,7 +272,7 @@ Cross product between two sets of hypervectors. - `V::AbstractVector{<:AbstractHV}`: Hypervectors # Examples -``` +```julia-repl julia> us = BinaryHV.('a':'e'; D = 10); julia> vs = BinaryHV.('v':'z'; D = 10); @@ -296,16 +296,18 @@ julia> crossproduct(us, vs) This encoding strategy first creates a multiset from both input hypervector sets, which are then bound together to generate all cross products, i.e. -U₁ × V₁ + U₁ × V₂ + ... + U₁ × Vₘ + ... + Uₙ × Vₘ +```math +U_1 \\times V_1 + U_1 \\times V_2 + ... + U_1 \\times V_m + ... + U_n \\times V_m +``` This encoding is based on the following formula: ```math -(\\oplus_{i=1}^{m} U_i) \\otimes (\\oplus_{i=1}^{n} V_i) +\\bigoplus_{i=1}^{m} U_i \\ \\otimes \\ \\bigoplus_{i=1}^{n} V_i ``` -where U and V are collections of hypervectors, `m` and `n` are the sizes of the U and V collections, -`ì` is the position in the hypervector collection, and `\\oplus` and `\\otimes` are the bundling +where \$U\$ and \$V\$ are collections of hypervectors, \$m\$ and \$n\$ are the sizes of the U and V collections, +\$i\$ is the position in the hypervector collection, and \$\\oplus\$ and \$\\otimes\$ are the bundling and binding operations. # References @@ -327,7 +329,7 @@ Creates a hypervector with the _n_-gram statistics of the input. - `n::Int = 3`: _n_-gram size # Examples -``` +```julia-repl julia> vs = BinaryHV.('a':'j'; D = 10); # a hypervector for each character julia> ngrams(vs) @@ -349,12 +351,12 @@ julia> ngrams(vs) This encoding is defined by the following mathematical notation: ```math -\\oplus_{i=1}^{m-n}\\otimes_{j=1}^{n-1}\\Pi^{n-j-1}(V_{i+j}) +\\bigoplus_{i=1}^{m-n}\\bigotimes_{j=1}^{n-1}\\rho^{n-j-1}(V_{i+j}) ``` -where `V` is the collection of hypervectors, `m` is the number of hypervectors in the -collection `V`, `n` is the window size, `i` is the position in the sequence, `j` is the -position in the _n_-gram, and `\\oplus`, `\\otimes` and `\\Pi` are the bundling, binding +where \$V\$ is the collection of hypervectors, \$m\$ is the number of hypervectors in the +collection \$V\$, \$n\$ is the window size, \$i\$ is the position in the sequence, \$j\$ is the +position in the _n_-gram, and \$\\oplus\$, \$\\otimes\$ and \$\\rho\$ are the bundling, binding and shift operations. !!! note @@ -386,26 +388,29 @@ end Graph for `source`-`target` pairs. Can be directed or undirected. # Arguments -- `source::T`: Source node hypervectors -- `target::T`: Target node hypervectors +- `source::AbstractVector{<:AbstractHV}`: Set of source node hypervectors +- `target::AbstractVector{<:AbstractHV}`: Set of target node hypervectors - `directed::Bool = false`: Whether the graph is directed or not # Example -``` -julia> nodes = BinaryHV.('a':'d'; D = 10); # a hypervector for each node -julia> graph(nodes[[1, 1, 2, 3]], nodes[[2, 3, 4, 4]]) # edges a-b, a-c, b-d, c-d -10-element BinaryHV with 2 true and 8 false: - 0 +```julia-repl +julia> V = [BinaryHV(i; D = 10) for i in 1:7]; + +julia> E = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4; 4 5; 5 6; 6 7]; # Lollipop graph + +julia> graph(V[E[:, 1]], V[E[:, 2]]) +10-element BinaryHV with 7 true and 3 false: 1 - 0 - 0 - 0 1 0 + 1 + 1 + 1 0 0 - 0 + 1 + 1 ``` # Extended help @@ -414,17 +419,17 @@ This encoding is based on the following mathematical notation: *Undirected graphs* ```math -\\oplus_{i=1}^{m} S_i \\otimes T_i +\\bigoplus_{i=1}^{E} V_i \\otimes V_j ``` *Directed graphs* ```math -\\oplus_{i=1}^{m} S_i \\otimes \\Pi(T_i) +\\bigoplus_{i=1}^{E} V_i \\otimes \\rho(V_j) ``` -where `K` and `V` are the key and value hypervector collections, `m` is the size of the -hypervector collection, `i` is the position of the entry in the collection, and `\\otimes`, -`\\oplus` and `\\Pi` are the binding, bundling and shift operations. +where \$V\$ is the node hypervector, \$i\$ and \$j\$ refer to the source and target nodes, +\$E\$ is the set of edges between nodes in the graph, and \$\\otimes\$, +\$\\oplus\$ and \$\\rho\$ are the binding, bundling and shift operations. # See also @@ -439,140 +444,3 @@ function graph(source::T, target::T; directed::Bool = false) where {T <: Abstrac @assert length(source) == length(target) "`source` and `target` must be the same length" return hashtable(source, shift.(target, convert(Int, directed))) end - - -""" - level(v::HV, n::Int) where {HV <: AbstractHV} - level(HV::Type{<:AbstractHV}, n::Int; D::Int = 10_000) - -Creates a set of level correlated hypervectors, where the first and last hypervectors are quasi-orthogonal. - -# Arguments -- `v::HV`: Base hypervector -- `m::Int`: Number of levels (alternatively, provide a vector to be encoded) -""" -function level(v::HV, m::Int) where {HV <: AbstractHV} - hvs = [v] - p = 2 / m - while length(hvs) < m - u = last(hvs) - push!(hvs, perturbate(u, p)) - end - return hvs -end - -level(HV::Type{<:AbstractHV}, n::Int; D::Int = 10_000) = level(HV(; D = D), n) -level(HVv, vals::AbstractVector) = level(HVv, length(vals)) -level(HVv, vals::UnitRange) = level(HVv, length(vals)) - - -""" - encodelevel(hvlevels::AbstractVector{<:AbstractHV}, numvalues; testbound=false) - -Generate an encoding function based on `level`, for encoding numerical values. It returns a function -that gives the corresponding hypervector for a given numerical input. - -# Arguments -- hvlevels::AbstractVector{<:AbstractHV}: vector of hypervectors representing the level encoding -- numvalues: the range or vector with the corresponding numerical values -- [testbound=false]: optional keyword argument to check whether the provided value is in bounds - -# Example -```julia -numvalues = range(0, 2pi, 100) -hvlevels = level(BipolarHV(), 100) - -encoder = encodelevel(hvlevels, numvalues) - -encoder(pi/3) # hypervector that best represents this numerical value -``` -""" -function encodelevel(hvlevels::AbstractVector{<:AbstractHV}, numvalues; testbound = false) - @assert length(hvlevels) == length(numvalues) "HV levels do not match numerical values" - # construct the encoder - function encoder(x::Number) - @assert !testbound || minimum(numvalues) ≤ x ≤ maximum(numvalues) "x not in numerical range" - (_, ind) = findmin(v -> abs(x - v), numvalues) - return hvlevels[ind] - end - return encoder -end - -""" - encodelevel(hvlevels::AbstractVector{<:AbstractHV}, a::Number, b::Number; testbound=false) - -See `encodelevel`, same but provide lower (`a`) and upper (`b`) limit of the interval to be encoded. -""" -encodelevel(hvlevels::AbstractVector{<:AbstractHV}, a::Number, b::Number; testbound = false) = encodelevel(hvlevels, range(a, b, length(hvlevels)); testbound) - -encodelevel(HV, numvalues; testbound = false) = encodelevel(level(HV, length(numvalues)), numvalues; testbound) - - -""" - decodelevel(hvlevels::AbstractVector{<:AbstractHV}, numvalues) - -Generate a decoding function based on `level`, for decoding numerical values. It returns a function -that gives the numerical value for a given hypervector, based on similarity matching. - -# Arguments -- hvlevels::AbstractVector{<:AbstractHV}: vector of hypervectors representing the level encoding -- numvalues: the range or vector with the corresponding numerical values - -# Example -```julia -numvalues = range(0, 2pi, 100) -hvlevels = level(BipolarHV(), 100) - -decoder = decodelevel(hvlevels, numvalues) - -decoder(hvlevels[17]) # value that closely matches the corresponding HV -``` -""" -function decodelevel(hvlevels::AbstractVector{<:AbstractHV}, numvalues; testbound = false) - # `testbound` is accepted (and ignored) for symmetry with `encodelevel`, so - # that the generic instance-based methods can forward keywords to both. - @assert length(hvlevels) == length(numvalues) "HV levels do not match numerical values" - # construct the decoder - function decoder(hv::AbstractHV) - (_, ind) = findmax(v -> similarity(v, hv), hvlevels) - return numvalues[ind] - end - return decoder -end - -decodelevel(hvlevels::AbstractVector{<:AbstractHV}, a::Number, b::Number; testbound = false) = decodelevel(hvlevels, range(a, b, length(hvlevels)); testbound) - -decodelevel(HV, numvalues; testbound = false) = decodelevel(level(HV, length(numvalues)), numvalues; testbound) - -""" - convertlevel(hvlevels, numvals..., kwargs...) - convertlevel(HV::AbstractHV, numvals..., kwargs...) - -Creates the `encoder` and `decoder` for a level encoding in one step. See `encodelevel` -and `decodelevel` for their respective documentations. -""" -convertlevel(hvlevels, numvals...; kwargs...) = encodelevel(hvlevels, numvals...; kwargs...), decodelevel(hvlevels, numvals...; kwargs...) - -convertlevel(hv::AbstractHV, numvals...; kwargs...) = encodelevel(hv, numvals...; kwargs...), decodelevel(hv, numvals...; kwargs...) - - -# levels using FHRR - -function level(v::FHRR, m::Int; β = 1 / m) - return [v^(x * β) for x in 1:m] -end - -function level(v::FHRR, vals::Union{AbstractVector{<:Number}, UnitRange}; β = 1 / (maximum(vals) - minimum(vals))) - return [v^(x * β) for x in vals] -end - -function encodelevel(v::FHRR, vals = (0, 1); β = 1 / (maximum(vals) - minimum(vals))) - return x -> v^(β * x) -end - -function decodelevel(v::FHRR, vals = (0, 1); β = 1 / (maximum(vals) - minimum(vals))) - return u -> @.(real(log(u.v) / log(v.v) / β)) |> mean -end - - -convertlevel(v::FHRR, vals = (0, 1); kwargs...) = encodelevel(v, vals; kwargs...), decodelevel(v, vals; kwargs...) diff --git a/src/inference.jl b/src/inference.jl index 63c853c..5cfbdd8 100644 --- a/src/inference.jl +++ b/src/inference.jl @@ -27,6 +27,81 @@ similarity(u::BinaryHV, v::BinaryHV) = sim_jacc(u, v) similarity(u::GradedHV, v::GradedHV) = sim_jacc(u, v) similarity(u::FHRR, v::FHRR) = real(dot(u.v, v.v)) / length(u) +""" + similaritymetric(HV) + similaritymetric(hv::AbstractHV) + +Which metric [`similarity`](@ref) uses for hypervectors of type `HV`, as a `Symbol`. + +`:jaccard` for [`BinaryHV`](@ref) and [`GradedHV`](@ref), whose elements are +non-negative, and `:cosine` for the types that can take negative values -- +[`BipolarHV`](@ref), [`TernaryHV`](@ref), [`RealHV`](@ref), +[`GradedBipolarHV`](@ref) and [`FHRR`](@ref). (For `FHRR` the normalized real part +of the Hermitian inner product *is* the cosine similarity, since every element has +unit modulus.) + +Use together with [`chancesimilarity`](@ref), which tells you what an +*unrelated* pair scores under that metric -- the number you need to judge whether a +given similarity is large. + +# Examples + +```jldoctest +julia> similaritymetric(BinaryHV), similaritymetric(BipolarHV) +(:jaccard, :cosine) + +julia> similaritymetric(BinaryHV(; D = 10)) # also works on an instance +:jaccard +``` + +# See also + +[`similarity`](@ref), [`chancesimilarity`](@ref) +""" +similaritymetric(::Type{<:AbstractHV}) = :cosine +similaritymetric(::Type{BinaryHV}) = :jaccard +similaritymetric(::Type{<:GradedHV}) = :jaccard +similaritymetric(hv::AbstractHV) = similaritymetric(typeof(hv)) + +""" + chancesimilarity(HV) + chancesimilarity(hv::AbstractHV) + +The expected [`similarity`](@ref) between two *independent random* hypervectors of +type `HV`: the value that means "unrelated". + +This is the baseline you need to interpret a similarity score, and it is **not +always zero**. Under cosine (see [`similaritymetric`](@ref)) unrelated hypervectors +score `0.0`, but under Jaccard they score `1/3` -- so a similarity of `0.35` is +strong evidence of a relationship for a [`BipolarHV`](@ref) and no evidence at all +for a [`BinaryHV`](@ref). + +!!! note + The value assumes hypervectors built with the type's default element + distribution. Passing a custom `distr` to [`RealHV`](@ref), [`GradedHV`](@ref) + or [`GradedBipolarHV`](@ref) can shift the baseline: a distribution that is not + centred makes random vectors point in a common direction, raising the chance + level above zero. + +# Examples + +```jldoctest +julia> chancesimilarity(BipolarHV) +0.0 + +julia> chancesimilarity(BinaryHV) # Jaccard: unrelated is 1/3, not 0 +0.3333333333333333 +``` + +# See also + +[`similarity`](@ref), [`similaritymetric`](@ref) +""" +chancesimilarity(::Type{<:AbstractHV}) = 0.0 +chancesimilarity(::Type{BinaryHV}) = 1 / 3 +chancesimilarity(::Type{<:GradedHV}) = 1 / 3 +chancesimilarity(hv::AbstractHV) = chancesimilarity(typeof(hv)) + """ similarity(u::AbstractVector, v::AbstractVector; method::Symbol) diff --git a/src/operations.jl b/src/operations.jl index c57b93b..cdc00bf 100644 --- a/src/operations.jl +++ b/src/operations.jl @@ -171,6 +171,21 @@ deterministic tie-breaking ([`BinaryHV`](@ref), [`BipolarHV`](@ref)), elementwis addition ([`TernaryHV`](@ref), [`RealHV`](@ref)), fuzzy aggregation ([`GradedHV`](@ref), [`GradedBipolarHV`](@ref)) or phasor addition ([`FHRR`](@ref)). +!!! warning "Bundling is m-way, not pairwise" + `bundle` combines **all** its inputs in one go, and for most types the result depends + on how many there are -- a majority vote over `m` inputs, a `1/√m` rescaling, a + renormalization. Bundling is therefore *not* associative for [`BinaryHV`](@ref), + [`BipolarHV`](@ref), [`RealHV`](@ref) and [`FHRR`](@ref) (it is for + [`TernaryHV`](@ref) and the graded types, whose rules happen to be associative). + + Three spellings do the right thing, because each reaches `bundle` with every input at + once: `bundle(hvs)`, `sum(hvs)`, and chained `x + y + z` (Julia parses a chain of `+` + as a single variadic call). + + Folding pairwise does **not**: `(x + y) + z`, `reduce(+, hvs)` and `foldl(+, hvs)` + bundle a bundle. The result is heavily biased towards the last inputs and loses the + defining property that a bundle is equally similar to each of its parts. + # See also [`bind`](@ref), [`similarity`](@ref) @@ -183,6 +198,12 @@ end Base.:+(u::HV, v::AbstractArray...) where {HV <: AbstractHV} = bundle((u, v...)) +# `sum` on a collection of hypervectors means bundling them, all at once. Without this +# method Base would fold pairwise with `+`, which bundles a bundle and biases the result +# towards the last elements (see the warning on `bundle`). +Base.sum(hvs::AbstractVector{<:AbstractHV}) = bundle(hvs) +Base.sum(hvs::Tuple{Vararg{AbstractHV}}) = bundle(hvs) + # BINDING # ------- """ @@ -253,8 +274,23 @@ Base.:/(hv1::HV, hv2::HV) where {HV <: AbstractHV} = unbind(hv1, hv2) # -------- # Shifting / Permutation +""" + shift!(hv::AbstractHV, k::Int) + +Permutes hypervector in-place by a specified number of shifts. + +This operations is used to assign an order to hypervectors. +""" shift!(hv::AbstractHV, k = 1) = (circshift!(hv.v, k); hv) + +""" + shift(hv::AbstractHV, k::Int) + +Permutes hypervector in-place by a specified number of shifts. + +This operations is used to assign an order to hypervectors. +""" function shift(hv::AbstractHV, k = 1) r = similar(hv) r.v .= circshift(hv.v, k) @@ -272,7 +308,18 @@ function shift(hv::V, k = 1) where {V <: Union{BinaryHV, BipolarHV}} return V(circshift!(v, hv.v, k)) end +""" + ρ(hv::AbstractHV, k::Int = 1) + +Alias of [`shift`](@ref). +""" ρ(hv::AbstractHV, k = 1) = shift(hv, k) + +""" + ρ!(hv::AbstractHV, k::Int = 1) + +Alias of [`shift!`](@ref). +""" ρ!(hv::AbstractHV, k = 1) = shift!(hv, k) @@ -392,7 +439,76 @@ function perturbate!(::Type{HVBitVec}, hv::AbstractHV, binargs; rng::AbstractRNG return hv end +""" + perturbate!(hv::AbstractHV, args...; rng::AbstractRNG = Random.GLOBAL_RNG) + +Perturbate hypervectors by randomly flipping values. + +Refer to [`perturbate`](@ref) for example use cases. +""" perturbate!(hv, args...; rng::AbstractRNG = Random.default_rng()) = perturbate!(vectype(hv), hv, args...; rng = rng) + +""" + perturbate(hv::AbstractHV, args...; rng::AbstractRNG = Random.GLOBAL_RNG) + +Perturbate hypervectors by randomly flipping values. + +# Examples + +```julia-repl +julia> v = BinaryHV(ones(10)) +10-element BinaryHV with 10 true and 0 false: + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + +julia> perturbate(v, 1) +10-element BinaryHV with 9 true and 1 false: + 1 + 1 + 1 + 1 + 1 + 0 + 1 + 1 + 1 + 1 + +julia> perturbate(v, 0.5) +10-element BinaryHV with 5 true and 5 false: + 0 + 1 + 1 + 0 + 1 + 0 + 0 + 1 + 0 + 1 + +julia> perturbate(v, [1,2,3]) +10-element BinaryHV with 7 true and 3 false: + 0 + 0 + 0 + 1 + 1 + 1 + 1 + 1 + 1 + 1 +``` +""" perturbate(hv::AbstractHV, args...; rng::AbstractRNG = Random.default_rng(), kwargs...) = perturbate!(copy(hv), args...; rng = rng, kwargs...) # OTHER diff --git a/src/representations.jl b/src/representations.jl index e22b1b2..13e4357 100644 --- a/src/representations.jl +++ b/src/representations.jl @@ -69,9 +69,7 @@ Render a hypervector as a square unicode heatmap of its leading `⌊√D⌋²` e Only available when UnicodePlots is loaded (`using UnicodePlots`); implemented in the UnicodePlotting package extension. -# See also - -[`unicodehistogram`](@ref) +See also [`unicodehistogram`](@ref). """ function unicodeheatmap end @@ -85,8 +83,6 @@ and [`TernaryHV`](@ref); phases for [`FHRR`](@ref)). Only available when UnicodePlots is loaded (`using UnicodePlots`); implemented in the UnicodePlotting package extension. -# See also - -[`unicodeheatmap`](@ref) +See also [`unicodeheatmap`](@ref). """ function unicodehistogram end diff --git a/src/types.jl b/src/types.jl index 2349671..d0415f1 100644 --- a/src/types.jl +++ b/src/types.jl @@ -98,7 +98,92 @@ Base.size(hv::AbstractHV) = size(hv.v) Base.sum(hv::AbstractHV) = sum(hv.v) LinearAlgebra.norm(hv::AbstractHV) = norm(hv.v) + +""" + normalize!(hv::AbstractHV) + +Normalize the input hypervector. + +# Examples + +```julia-repl +julia> v = bundle([encode(TernaryHV, i; D = 8) for i in 1:100]) +8-element TernaryHV{Int64} with 4 positives, 0 zeros, and 4 negatives: + 12 + -2 + -16 + 12 + -4 + 6 + 12 + -12 + +julia> normalize!(v) +8-element TernaryHV{Int64} with 4 positives, 0 zeros, and 4 negatives: + 1 + -1 + -1 + 1 + -1 + 1 + 1 + -1 + +julia> v +8-element TernaryHV{Int64} with 4 positives, 0 zeros, and 4 negatives: + 1 + -1 + -1 + 1 + -1 + 1 + 1 + -1 +""" normalize!(hv::AbstractHV) = hv + +""" + normalize(hv::AbstractHV) + +Return a copy of the normalized version of the input hypervector. + +# Example + +```julia-repl +julia> v = bundle([encode(TernaryHV, i; D = 8) for i in 1:100]) +8-element TernaryHV{Int64} with 4 positives, 0 zeros, and 4 negatives: + 12 + -2 + -16 + 12 + -4 + 6 + 12 + -12 + +julia> normalize(v) +8-element TernaryHV{Int64} with 4 positives, 0 zeros, and 4 negatives: + 1 + -1 + -1 + 1 + -1 + 1 + 1 + -1 + +julia> v +8-element TernaryHV{Int64} with 4 positives, 0 zeros, and 4 negatives: + 12 + -2 + -16 + 12 + -4 + 6 + 12 + -12 +``` +""" normalize(hv::AbstractHV) = (c = copy(hv); normalize!(c); c) # Types that carry an element distribution (`distr` field) override the diff --git a/test/encoding.jl b/test/encoding.jl index b9eb8b2..f4aada2 100644 --- a/test/encoding.jl +++ b/test/encoding.jl @@ -98,49 +98,246 @@ end @test_throws AssertionError graph(hvs[s], hvs[[1, 2, 3]]) end - @testset "levels" begin - numvals = 0:0.1:2pi - levels = level(BinaryHV(; D = 100), numvals) - - @test length(levels) == length(numvals) - @test eltype(levels) <: BinaryHV - - encoder, decoder = convertlevel(levels, numvals) - hv = encoder(1.467) - @test hv isa BinaryHV - x = decoder(hv) - @test 1 ≤ x ≤ 2 - - # regression (TODO §1.4): the instance-based decodelevel/convertlevel - # path used to throw a MethodError caused by kwarg forwarding - dec = decodelevel(BinaryHV(; D = 100, seed = 3), numvals) - @test dec isa Function - enc2, dec2 = convertlevel(BinaryHV(; D = 100, seed = 3), numvals) - @test enc2(1.0) isa BinaryHV - @test minimum(numvals) ≤ dec2(enc2(1.0)) ≤ maximum(numvals) - # keywords also forward through the shared-ladder path - enc3, dec3 = convertlevel(levels, numvals; testbound = true) - @test dec3(enc3(1.467)) isa Number - # NOTE (TODO §1.4b): enc2/dec2 are built over *different* random ladders, - # so decode(encode(x)) ≈ x does not hold on the instance path — flagged, - # not asserted here. + @testset "LevelEncoder" begin + @testset "ladder for all types: $HV" for HV in + (BinaryHV, BipolarHV, TernaryHV, RealHV, GradedHV, GradedBipolarHV, FHRR) + lvl = LevelEncoder(HV, (0, 1), 20; D = 1_000, seed = 17) + @test lvl isa LevelEncoder{<:HV} + @test lvl.base === nothing # the ladder mechanism, even for FHRR + @test length(lvl.levels) == length(lvl.values) == 20 + @test encode(lvl, 0.5) isa HV + # round-trips within one quantization step + step = 1 / 19 + for x in (0.0, 0.31, 0.77, 1.0) + @test abs(decode(lvl, encode(lvl, x)) - x) ≤ step + end + # adjacent levels are more similar than distant ones + @test similarity(lvl.levels[1], lvl.levels[2]) > + similarity(lvl.levels[1], lvl.levels[end]) + end + + @testset "one shared level set (regression TODO §1.4/§1.4b)" begin + # every encode/decode call draws from the level set built at + # construction, so separate calls are comparable by construction — + # the incomparability bug of the old encodelevel/decodelevel + # function family cannot occur + lvl = LevelEncoder(BipolarHV, (0, 2π), 50; D = 2_000, seed = 3) + @test encode(lvl, 1.0) == encode(lvl, 1.0) + a, b, c = encode(lvl, 1.0), encode(lvl, 1.1), encode(lvl, 5.0) + @test similarity(a, b) > similarity(a, c) + # encode and decode share ONE ladder: exact round-trip on the grid + # (the old instance path built two ladders; error was up to 1.0) + for v in lvl.values[[1, 10, 25, 50]] + @test decode(lvl, encode(lvl, v)) == v + end + # seeded construction is fully deterministic + @test LevelEncoder(BipolarHV, (0, 1), 5; D = 100, seed = 9).levels == + LevelEncoder(BipolarHV, (0, 1), 5; D = 100, seed = 9).levels + # the old function family is gone, not just patched + for f in (:level, :encodelevel, :decodelevel, :convertlevel) + @test !isdefined(HyperdimensionalComputing, f) + end + end + + @testset "bandwidth controls similarity decay" begin + ends(lvl) = similarity(lvl.levels[1], lvl.levels[end]) + narrow = LevelEncoder(BipolarHV, (0, 1), 20; D = 2_000, bandwidth = 0.02, seed = 7) + wide = LevelEncoder(BipolarHV, (0, 1), 20; D = 2_000, bandwidth = 0.3, seed = 7) + @test ends(narrow) > ends(wide) + # FPE: β is the analogous knob + slow = LevelEncoder(FHRR, 0:0.1:1; β = 0.05, D = 2_000, seed = 7) + fast = LevelEncoder(FHRR, 0:0.1:1; β = 1.0, D = 2_000, seed = 7) + @test similarity(encode(slow, 0.0), encode(slow, 1.0)) > + similarity(encode(fast, 0.0), encode(fast, 1.0)) + end + + @testset "fractional power encoding (FHRR)" begin + fpe = LevelEncoder(FHRR, 0:0.1:10; D = 1_000, seed = 11) + @test fpe.base isa FHRR + hv = encode(fpe, 3.14) # continuous: no quantization + @test hv isa FHRR + @test all(abs.(hv.v) .≈ 1) + # nearest-neighbour decode lands on the grid, within one step + @test abs(decode(fpe, hv) - 3.14) ≤ 0.1 + # analytic decode is continuous and exact on a clean vector + @test decode(fpe, hv; method = :analytic) ≈ 3.14 + # similarity decays with distance + @test similarity(encode(fpe, 2.0), encode(fpe, 2.2)) > + similarity(encode(fpe, 2.0), encode(fpe, 8.0)) + end + + @testset "precomputed levels constructor" begin + src = LevelEncoder(BinaryHV, (0, 1), 10; D = 500, seed = 5) + lvl = LevelEncoder(src.levels, src.values) + @test decode(lvl, encode(lvl, 0.42)) == decode(src, encode(src, 0.42)) + @test_throws ArgumentError LevelEncoder(src.levels, 1:3) # length mismatch + # analytic decoding is FPE-only + @test_throws ArgumentError decode(lvl, encode(lvl, 0.5); method = :analytic) + @test_throws ArgumentError decode(lvl, encode(lvl, 0.5); method = :nonsense) + end + + @testset "bounds and argument checks" begin + lvl = LevelEncoder(BinaryHV, (0, 1), 10; D = 200, seed = 1) + @test encode(lvl, 3.0) == lvl.levels[end] # snaps to nearest by default + @test_throws DomainError encode(lvl, 3.0; testbound = true) + @test encode(lvl, 0.5; testbound = true) isa BinaryHV + @test_throws ArgumentError LevelEncoder(BinaryHV, (0, 1), 1) + @test_throws ArgumentError LevelEncoder(BinaryHV, (0, 1), 10; bandwidth = 0) + @test_throws ArgumentError LevelEncoder(BinaryHV, (0, 1), 10; bandwidth = 1.5) + end end - @testset "FHRR numbers" begin + @testset "RandomProjection" begin + x = [0.9, -0.2, 0.4] + far = [-1.5, 2.0, -0.3] + + indomain = Dict( + BinaryHV => hv -> all(b -> b == 0 || b == 1, collect(hv)), + BipolarHV => hv -> all(b -> b == -1 || b == 1, collect(hv)), + TernaryHV => hv -> all(b -> b in (-1, 0, 1), collect(hv)), + RealHV => hv -> all(isfinite, collect(hv)), + GradedHV => hv -> all(b -> 0 ≤ b ≤ 1, collect(hv)), + GradedBipolarHV => hv -> all(b -> -1 ≤ b ≤ 1, collect(hv)), + FHRR => hv -> all(z -> abs(z) ≈ 1, hv.v), + ) + + @testset "comparability and locality: $HV" for HV in keys(indomain) + rp = RandomProjection(HV, 3; D = 2_000, seed = 21) + # separate encode calls share R, so they are comparable and + # deterministic ... + @test encode(rp, x) == encode(rp, x) + # ... and nearby feature vectors are more similar than distant ones + @test similarity(encode(rp, x), encode(rp, x .+ 0.02)) > + similarity(encode(rp, x), encode(rp, far)) + end + + @testset "matrix distributions stay in domain: $HV" for HV in keys(indomain) + for matrix in (:gaussian, :bipolar, :sparse_ternary) + rp = RandomProjection(HV, 3; D = 500, matrix, seed = 22) + @test indomain[HV](encode(rp, x)) + end + end + + @testset "supplied matrix" begin + rp = RandomProjection(BipolarHV, 3; D = 200, seed = 6) + # reusing the drawn matrix reproduces identical encodings + rp2 = RandomProjection(BipolarHV, rp.R) + @test encode(rp2, x) == encode(rp, x) + # d and D are inferred from the supplied matrix + @test size(rp2.R) == (200, 3) + # the ternary supplied-matrix path (no target_sparsity keyword) + rpt = RandomProjection(TernaryHV, rp.R; θ = 0.1) + @test rpt.θ == 0.1 + @test encode(rpt, x) isa TernaryHV + # a zero projection row must not throw: z = 0 has a defined image + rpz = RandomProjection(BipolarHV, [0.0 0.0 0.0; 1.0 -2.0 0.5]) + @test collect(encode(rpz, x))[1] == -1 + @test collect(encode(RandomProjection(TernaryHV, [0.0 0.0 0.0]), x))[1] == 0 + end + + @testset "sparse_ternary matrix is sparse and ternary" begin + R = HyperdimensionalComputing.projection_matrix(:sparse_ternary, 1_000, 9; rng = Xoshiro(5)) + @test all(r -> r in (-1.0, 0.0, 1.0), R) + # nonzero density 1/√d = 1/3 + @test count(!iszero, R) / length(R) ≈ 1 / 3 atol = 0.03 + @test_throws ArgumentError RandomProjection(BinaryHV, 3; matrix = :nope) + end + + @testset "ternary data-driven threshold" begin + X = randn(Xoshiro(1), 4, 100) + rp = RandomProjection(TernaryHV, X; target_sparsity = 0.7, D = 1_000, seed = 2) + encoded = encode(rp, X) # per-column, all comparable + @test length(encoded) == 100 + @test encoded[7] == encode(rp, X[:, 7]) + zerofrac = sum(hv -> count(iszero, collect(hv)), encoded) / (1_000 * 100) + @test zerofrac ≈ 0.7 atol = 0.01 + @test !ismutable(rp) # construction from data, not fitting + @test_throws ArgumentError RandomProjection(TernaryHV, X; target_sparsity = 1.2) + end + + @testset "ternary constructor: positional collision" begin + # `RandomProjection(TernaryHV, M)` reads M as a supplied projection + # matrix; `RandomProjection(TernaryHV, M; target_sparsity)` reads M + # as training data. Same positional shape, opposite meanings of the + # matrix — these tests pin the documented reading of each path and + # that a misread cannot silently encode the intended features. + R = randn(Xoshiro(11), 50, 4) + x4 = [0.3, -1.0, 0.8, 0.1] - v = FHRR() + # supplied-matrix path: D and d are read off size(R), and encode + # projects through R itself (R is not treated as data) + rp = RandomProjection(TernaryHV, R; θ = 0.2) + @test size(rp.R) == (50, 4) + z = R * x4 + @test collect(encode(rp, x4)) == Int.(sign.(z) .* (abs.(z) .> 0.2)) - numvals = 0:0.1:10 + # data-driven path: M is d × n data, so the drawn R has D rows and + # size(M, 1) columns, and the sparsity target is hit + X = randn(Xoshiro(12), 4, 60) + rpd = RandomProjection(TernaryHV, X; target_sparsity = 0.3, D = 500, seed = 13) + @test size(rpd.R) == (500, 4) + zerofrac = sum(hv -> count(iszero, collect(hv)), encode(rpd, X)) / (500 * 60) + @test zerofrac ≈ 0.3 atol = 0.01 - encoder, decoder = convertlevel(v, numvals) + # the dangerous misread: data passed WITHOUT target_sparsity is + # taken as a projection matrix (here D = 4, d = 60). That reading + # is coherent on its own terms, but the intended 4-feature vectors + # CANNOT be encoded through it — the mistake surfaces as an + # immediate error at first use, never as silently wrong + # hypervectors of the intended features + oops = RandomProjection(TernaryHV, X) + @test size(oops.R) == (4, 60) + @test_throws DimensionMismatch encode(oops, x4) + @test encode(oops, randn(Xoshiro(14), 60)) isa TernaryHV + # NOTE: a *square* matrix is genuinely ambiguous — both readings + # yield a working encoder — and cannot be caught by tests without + # changing the constructor. Tracked as a design item in TODO.md. + end - x, y, z = 2, 5, 10 + @testset "scalar and vector θ, rethreshold" begin + rp = RandomProjection(TernaryHV, 3; D = 400, θ = 0.5, seed = 8) + rpv = rethreshold(rp, fill(0.5, 400)) + @test rpv.R === rp.R # rethreshold shares the projection matrix + @test encode(rpv, x) == encode(rp, x) # uniform vector θ ≡ scalar θ + # a genuinely per-component θ changes only its components + rp0 = rethreshold(rp, 0.0) + @test encode(rethreshold(rp, zeros(400)), x) == encode(rp0, x) + @test_throws ArgumentError rethreshold(rp, ones(7)) # must have length D + @test_throws ArgumentError RandomProjection(BinaryHV, 3; β = -1) + end - hx, hy, hz = encoder.((x, y, z)) + @testset "FHRR is random Fourier features (shared phase helper)" begin + rff = RandomProjection(FHRR, 3; D = 500, β = 0.5, seed = 4) + # the nonlinearity IS the shared helper, not a copy of it + @test encode(rff, x) == HyperdimensionalComputing.phase_encode(rff.R * x, 0.5) + # ... which LevelEncoder's fractional power path also routes through + fpe = LevelEncoder(FHRR, 0:0.1:1; D = 500, seed = 4) + @test encode(fpe, 0.3) == + HyperdimensionalComputing.phase_encode(angle.(fpe.base.v), fpe.bandwidth * 0.3) + # similarity approximates the Gaussian kernel exp(-β²‖x-y‖²/2) + kern = RandomProjection(FHRR, 3; β = 0.5, seed = 1) + y = x .+ [0.2, 0.0, 0.0] + @test similarity(encode(kern, x), encode(kern, y)) ≈ exp(-(0.5 * 0.2)^2 / 2) atol = 0.02 + end - @test hx isa FHRR - @test similarity(hx, hy) > similarity(hx, hz) + @testset "decode is clean-up, never inversion" begin + X = randn(Xoshiro(3), 3, 20) + rp = RandomProjection(BipolarHV, 3; D = 1_000, seed = 9) + refs = encode(rp, X) + τ, i, nn = decode(rp, refs[13], refs) + @test i == 13 && nn === refs[13] && τ ≈ 1 + # noisy queries still clean up to the right reference + @test decode(rp, encode(rp, X[:, 13] .+ 0.01), refs)[2] == 13 + # no codebook, no decode: a random projection is lossy + @test_throws ArgumentError decode(rp, refs[1]) + @test_throws ArgumentError decode(rp, refs[1], refs; method = :analytic) + @test_throws ArgumentError decode(rp, refs[1], refs; method = :nonsense) + end - @test decoder(hx) < decoder(hy) < decoder(hz) + @testset "argument checks" begin + rp = RandomProjection(BinaryHV, 3; D = 100, seed = 10) + @test_throws DimensionMismatch encode(rp, [1.0, 2.0]) + @test_throws ArgumentError RandomProjection(BinaryHV, 0) + end end end diff --git a/test/inference.jl b/test/inference.jl index df47e21..51763f1 100644 --- a/test/inference.jl +++ b/test/inference.jl @@ -71,8 +71,30 @@ @test δ === similarity end + @testset "similarity metric traits" begin + @test similaritymetric(BinaryHV) == :jaccard + @test similaritymetric(GradedHV) == :jaccard + for HV in [BipolarHV, TernaryHV, RealHV, GradedBipolarHV, FHRR] + @test similaritymetric(HV) == :cosine + end + ## instance forms agree with the type forms + @test similaritymetric(BinaryHV(; D = 10)) == :jaccard + @test chancesimilarity(BinaryHV(; D = 10)) == chancesimilarity(BinaryHV) + + ## the documented baseline must match what unrelated hypervectors actually score + for HV in [BinaryHV, BipolarHV, TernaryHV, RealHV, GradedHV, GradedBipolarHV, FHRR] + measured = mean(similarity(HV(; D = 5_000), HV(; D = 5_000)) for _ in 1:50) + @test isapprox(measured, chancesimilarity(HV); atol = 0.02) + end + + ## FHRR's similarity really is cosine, as the docstring claims + a, b = FHRR(; D = 1000), FHRR(; D = 1000) + av, bv = collect(a), collect(b) + @test similarity(a, b) ≈ real(dot(av, bv)) / (norm(av) * norm(bv)) + end + @testset "Similarity matrix" begin - levels = level(RealHV(; D = 100), 10) + levels = LevelEncoder(RealHV, (0, 1), 10; D = 100).levels M = similarity(levels) @test M isa Matrix @test size(M) == (10, 10) diff --git a/test/operations.jl b/test/operations.jl index 2ef0315..8798ef4 100644 --- a/test/operations.jl +++ b/test/operations.jl @@ -145,6 +145,30 @@ Random.seed!(42) end end + # Bundling is m-way: it combines all inputs at once, and for most types the rule + # depends on how many there are. Chained `+` and `sum` must reach `bundle` with the + # whole collection; folding pairwise bundles a bundle and is NOT equivalent. + @testset "bundle is m-way, not pairwise" begin + for HV in [BinaryHV, BipolarHV, TernaryHV, RealHV, GradedHV, GradedBipolarHV, FHRR] + x = encode(HV, :a; D = 1000) + y = encode(HV, :b; D = 1000) + z = encode(HV, :c; D = 1000) + reference = bundle([x, y, z]) + ## Julia parses a chain of `+` as one variadic call, so this is a 3-way bundle + @test x + y + z == reference + @test sum([x, y, z]) == reference + @test sum((x, y, z)) == reference + end + + ## the semantic property: a bundle is equally similar to each of its parts, + ## whereas pairwise folding is dominated by the last inputs + hvs = [encode(BipolarHV, i; D = 10_000) for i in 1:5] + spread(hv) = (sims = [similarity(hv, h) for h in hvs]; maximum(sims) - minimum(sims)) + @test spread(bundle(hvs)) < 0.1 + @test spread(reduce(+, hvs)) > 0.2 + @test similarity(reduce(+, hvs), last(hvs)) > similarity(reduce(+, hvs), first(hvs)) + end + @testset "unbind" begin # XOR- and multiplication-based types: binding is self-inverse, roundtrip exact for HV in [BinaryHV, BipolarHV, TernaryHV] @@ -207,8 +231,9 @@ Random.seed!(42) hvm = perturbate(hv1, m; rng = Xoshiro(1)) @test all(hvm.v[.!m] .== hv1.v[.!m]) @test all(abs.(hvm.v[m]) .≈ 1) - # `level` for FHRR uses its own ^-based method, not perturbation - @test level(FHRR(; D = 100), 5) isa Vector{<:FHRR} + # the perturbation ladder (LevelEncoder with an explicit level count) + # works for FHRR thanks to the phase-resampling methods above + @test LevelEncoder(FHRR, (0, 1), 5; D = 100).levels isa Vector{<:FHRR} end end diff --git a/test/types.jl b/test/types.jl index 7096d56..9a0de5d 100644 --- a/test/types.jl +++ b/test/types.jl @@ -167,6 +167,34 @@ using Distributions @test occursin("$(count(==(-1), collect(x))) negatives", s) end + # Sharp edges of the `false ↦ +1 / true ↦ -1` storage convention. These + # lock CURRENT, deliberate behaviour — the convention is load-bearing + # (XOR on stored bits IS the ±1 product) and each trap below has caused a + # real bug (the original polarity flip; RandomProjection's bipolar + # nonlinearity). Do not "fix" these; they are documented-by-test so the + # obvious-but-wrong idiom trips a red test instead of shipping silent + # corruption. (Bind identity, construction round-trip and the zero-throws + # are locked in the two BipolarHV testsets above.) + @testset "BipolarHV sharp edges" begin + # `sign.(z)` is NOT a valid real-to-bipolar map: sign(0) = 0, and zero + # has no bipolar state — this is the trap RandomProjection's bipolar + # nonlinearity had to avoid + @test_throws ArgumentError BipolarHV(sign.([0.5, -1.2, 0.0, 2.0])) + + # Bool vectors are RAW STORED BITS, not logical positivity: + # true ↦ -1, false ↦ +1 — the opposite of what "true is positive" + # intuition expects, deliberately, so bind-as-XOR is the exact product + @test collect(BipolarHV([true, false, true])) == [-1, 1, -1] + + # hence the correct thresholding of reals is an explicit ±1 ifelse; + # the "obvious" BipolarHV(z .> 0) hits the raw-bits path and silently + # yields the OPPOSITE polarity — both are pinned so the difference is + # on record + z = [0.5, -1.2, 2.0] # no zeros: only polarity is at issue here + @test collect(BipolarHV(ifelse.(z .> 0, 1, -1))) == [1, -1, 1] + @test collect(BipolarHV(z .> 0)) == [-1, 1, -1] # inverted! never reach for this + end + @testset "BinaryHV" begin hdv = BinaryHV(; D = n)