diff --git a/docs/pages.jl b/docs/pages.jl index a15e17918..ee608e185 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -14,10 +14,12 @@ pages = [ ], "Examples" => Any[ "Building a model to add to ReservoirComputing.jl" => "examples/model_es2n.md", + "Morphing patterns with conceptors" => "examples/conceptors_morphing.md", ], "API Documentation" => Any[ "Layers" => "api/layers.md", "Models" => "api/models.md", + "Conceptors" => "api/conceptors.md", "Utilities" => "api/utils.md", "Train" => "api/train.md", "Predict" => "api/predict.md", diff --git a/docs/src/api/conceptors.md b/docs/src/api/conceptors.md new file mode 100644 index 000000000..79d8012cb --- /dev/null +++ b/docs/src/api/conceptors.md @@ -0,0 +1,58 @@ +# Conceptors + +Conceptors after Jaeger (2014), [Jaeger2014conceptors](@cite). See the +[Morphing patterns with conceptors](@ref) example for an end-to-end walkthrough. + +## Conceptor matrices + +```@docs + conceptor_matrix + correlation_matrix + conceptor_from_states + conceptor_singular_values + quota +``` + +## Aperture adaptation + +```@docs + aperture_adapt + adapt_singular_value + reaperture + attenuation + optimal_aperture +``` + +## Boolean algebra + +```@docs + conceptor_not + conceptor_and + conceptor_or +``` + +## The `Conceptor` wrapper and its library + +```@docs + Conceptor + has_conceptor + get_conceptor + store_conceptor! + set_active_conceptor + active_conceptor +``` + +## Loading, generation, and morphing + +```@docs + load! + generate + morph_conceptor + ridge_map +``` + +## Conceptor-filtered training + +```@docs + store_conceptors! +``` diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml new file mode 100644 index 000000000..bdba183c7 --- /dev/null +++ b/docs/src/assets/Project.toml @@ -0,0 +1,41 @@ +[deps] +CellularAutomata = "878138dc-5b27-11ea-1a71-cb95d38d6b29" +ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" +DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" +DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" +DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +LIBSVM = "b1bec4e5-fd48-53fe-b0cb-9723c09d164b" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +MLJLinearModels = "6ee0df7b-362f-4a72-a706-9e79364fb692" +OrdinaryDiffEqAdamsBashforthMoulton = "89bda076-bce5-4f1c-845f-551c83cdda9a" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +ReservoirComputing = "7c2d2b1e-3dd4-11ea-355a-8f6a8116e294" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" + +[compat] +CellularAutomata = "0.0.6" +ConcreteStructs = "0.2" +DataInterpolations = "6, 7, 8" +DelayDiffEq = "6" +Documenter = "1" +DocumenterCitations = "1" +DocumenterInterLinks = "1" +JLD2 = "0.6" +LIBSVM = "0.8" +LinearSolve = "3" +MLJLinearModels = "0.10" +OrdinaryDiffEqAdamsBashforthMoulton = "2" +OrdinaryDiffEqTsit5 = "1, 2" +Plots = "1" +ReservoirComputing = "0.12.0" +SciMLBase = "2.51, 3" +Static = "1" +StatsBase = "0.34.4" diff --git a/docs/src/examples/conceptors_morphing.md b/docs/src/examples/conceptors_morphing.md new file mode 100644 index 000000000..2a3251742 --- /dev/null +++ b/docs/src/examples/conceptors_morphing.md @@ -0,0 +1,127 @@ +# Morphing patterns with conceptors + +This example reproduces the morphing square from Jaeger (2014), *Controlling +Recurrent Neural Networks by Conceptors* ([Jaeger2014conceptors](@cite), +arXiv:1403.3369, Figure 2). A single reservoir is *loaded* with four patterns; +their [`conceptor`](@ref conceptor_matrix)s are then linearly **morphed** and used +to drive the reservoir autonomously, generating a continuum of patterns that +interpolate between, and extrapolate beyond, the four originals. + +A conceptor is a positive semidefinite matrix ``C = R (R + \alpha^{-2} I)^{-1}`` +built from the correlation matrix ``R`` of a reservoir's driven states. Inserted +into the autonomous update ``x(n) = C\,\tanh(W x(n-1) + b)``, it constrains the +reservoir dynamics to the subspace excited by one loaded pattern, so the reservoir +re-generates that pattern. A linear combination of conceptors morphs between them. + +## The four driving patterns + +Following the original demonstration we use two sines of slightly different +period and two minor variations of a 5-periodic random pattern, all in +``[-0.9, 0.9]``. The two period-5 vectors are the ones from the published figure. + +```@example morphing +using ReservoirComputing +using Random +using Plots + +const PERIOD_1 = 8.8342522 +const PERIOD_2 = 9.8342522 +const DRIVER_LEN = 1500 + +# period-5 patterns, sampled as rp[mod(n,5)+1] (values from Jaeger 2014, Fig. 1B) +const RAND5 = [0.16, 0.9, -0.9, -0.21, -0.55] +const PERTURB5 = [0.10, 0.9, -0.9, -0.65, -0.54] +period5(rp, n) = rp[mod(n, 5) + 1] + +n = 1:DRIVER_LEN +patterns = [ + :s1 => reshape(sin.(2π .* n ./ PERIOD_1), 1, :), + :s2 => reshape(sin.(2π .* n ./ PERIOD_2), 1, :), + :r3 => reshape([period5(RAND5, i) for i in n], 1, :), + :r4 => reshape([period5(PERTURB5, i) for i in n], 1, :), +] +nothing # hide +``` + +## Building and loading the reservoir + +We build a 100-unit [`ESN`](@ref) with the scalings used in the report +(spectral radius 1.5, input scaling 1.5, bias scaling 0.2) and wrap it in a +[`Conceptor`](@ref). [`load!`](@ref) drives the reservoir with each pattern, +stores a conceptor for it, and recomputes the recurrent weights into an +input-internalizing matrix so the reservoir can run autonomously. + +```@example morphing +rng = Xoshiro(3) + +bias_init(r, dims...) = 0.2f0 .* randn(r, Float32, dims...) +input_init(r, dims...) = 1.5f0 .* randn(r, Float32, dims...) +res_init(r, dims...) = rand_sparse(r, dims...; radius = 1.5f0, sparsity = 0.1f0) + +esn = ESN(1, 100, 1; use_bias = true, init_bias = bias_init, + init_input = input_init, init_reservoir = res_init) +concept = Conceptor(esn) +ps = initialparameters(rng, concept) +st = initialstates(rng, concept) + +ps, st = load!(rng, concept, patterns, ps, st; aperture = 4.0, washout = 500) +nothing # hide +``` + +## Morphing across the square + +For a grid of mixing coordinates ``a, b \in \{-0.5, \dots, 1.5\}`` we form the +morphed conceptor ``M = \mu_1 C_{s1} + \mu_2 C_{s2} + \mu_3 C_{r3} + \mu_4 C_{r4}`` +with +```math +\mu_1 = (1-a)\,b,\quad \mu_2 = a\,b,\quad \mu_3 = (1-a)(1-b),\quad \mu_4 = a\,(1-b), +``` +(the coefficients always sum to one), and run the reservoir autonomously under +``M`` with [`generate`](@ref). All panels share one start state so their phases are +comparable. + +```@example morphing +morph_weights(a, b) = (; s1 = (1 - a) * b, s2 = a * b, + r3 = (1 - a) * (1 - b), r4 = a * (1 - b)) + +grid = collect(-0.5:0.25:1.5) +x0 = rand(Xoshiro(4), 100) +prototypes = Set([(0.0, 1.0), (1.0, 1.0), (0.0, 0.0), (1.0, 0.0)]) + +plt = plot(layout = (length(grid), length(grid)), size = (1000, 1000), + legend = false, framestyle = :box, ticks = false, link = :all) + +for (ib, b) in enumerate(grid), (ia, a) in enumerate(grid) + M = morph_conceptor(st, morph_weights(a, b)) + Y, _ = generate(concept, ps, st; conceptor = M, steps = 15, + washout = 190, init_state = x0) + k = (ib - 1) * length(grid) + ia # row ib (top = b = -0.5), col ia + proto = (a, b) in prototypes + plot!(plt, subplot = k, vec(Y); color = :black, linewidth = proto ? 2.5 : 1.2, + ylims = (-1, 1), background_color_subplot = proto ? :gray88 : :white) +end +plt +``` + +The four shaded panels are the loaded prototypes: the two sines sit on the +``b = 1`` row, the two period-5 patterns on the ``b = 0`` row. The inner block +(``a, b \in [0, 1]``) interpolates smoothly between them; the outer panels +extrapolate, producing distorted, higher-frequency variants. + +## Single-pattern recall and the conceptor algebra + +Loading also lets the reservoir recall a single pattern by name, and the stored +conceptors support an aperture adaptation ([`aperture_adapt`](@ref)) and a Boolean +algebra ([`conceptor_and`](@ref), [`conceptor_or`](@ref), [`conceptor_not`](@ref)). + +```@example morphing +# autonomous recall of the first sine +Yrec, _ = generate(concept, ps, st; conceptor = :s1, steps = 100, washout = 200) + +# conceptor algebra: shared vs. combined subspaces of the two sines +Cs1, Cs2 = get_conceptor(st, :s1), get_conceptor(st, :s2) +(quota(conceptor_and(Cs1, Cs2)), quota(Cs1), quota(conceptor_or(Cs1, Cs2))) +``` + +The triple is increasing — `AND` keeps only the shared directions, `OR` spans the +union — illustrating the lattice structure of the conceptor algebra. diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 0f3946ab0..2882e3e35 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -1,3 +1,12 @@ +@article{Jaeger2014conceptors, + title = {Controlling Recurrent Neural Networks by Conceptors}, + url = {https://arxiv.org/abs/1403.3369}, + DOI = {10.48550/arXiv.1403.3369}, + journal = {arXiv preprint arXiv:1403.3369}, + author = {Jaeger, Herbert}, + year = {2014}, +} + @article{Lu2017, title = {Reservoir observers: Model-free inference of unmeasured variables in chaotic systems}, volume = {27}, diff --git a/src/ReservoirComputing.jl b/src/ReservoirComputing.jl index 3bec529aa..34962ccdd 100644 --- a/src/ReservoirComputing.jl +++ b/src/ReservoirComputing.jl @@ -2,12 +2,13 @@ module ReservoirComputing using ArrayInterface: ArrayInterface using ConcreteStructs: @concrete -using LinearAlgebra: eigvals, I, qr, Diagonal, diag, mul!, Symmetric, norm +using LinearAlgebra: eigvals, eigen, I, qr, Diagonal, diag, mul!, Symmetric, norm, + svd, svdvals, nullspace, pinv, tr, checksquare, issymmetric using LuxCore: AbstractLuxLayer, AbstractLuxContainerLayer, AbstractLuxWrapperLayer, setup, apply, replicate import LuxCore: initialparameters, initialstates, statelength, outputsize using NNlib: tanh_fast -using Random: Random, AbstractRNG, randperm +using Random: Random, AbstractRNG, randperm, randn using Static: StaticBool, StaticSymbol, True, False, static, known, StaticInteger using Reexport: Reexport, @reexport using WeightInitializers: WeightInitializers, DeviceAgnostic, PartialFunction, Utils, @@ -59,6 +60,8 @@ include("models/lifesn.jl") include("models/ngrc.jl") include("models/rmnesn.jl") include("models/rmnresesn.jl") +#conceptors +include("conceptors.jl") #extensions include("extensions/reca.jl") @@ -85,6 +88,14 @@ export polynomial_monomials, chebyshev_monomials, predict, QRSolver, resetcarry! export AdditiveEIESN, DeepESN, DelayESN, EIESN, ES2N, ESN, EuSN, HybridESN, InputDelayESN, LIFESN, ResESN, StateDelayESN, SVESM export NGRC export RMNESN, RMNResESN +#conceptors +export Conceptor, correlation_matrix, conceptor_matrix, conceptor_from_states, + conceptor_singular_values, quota +export aperture_adapt, adapt_singular_value, reaperture, attenuation, optimal_aperture +export conceptor_not, conceptor_and, conceptor_or +export has_conceptor, get_conceptor, store_conceptor!, store_conceptors!, + set_active_conceptor, active_conceptor +export ridge_map, load!, generate, morph_conceptor #ext export RECACell, RECA export RandomMapping, RandomMaps diff --git a/src/conceptors.jl b/src/conceptors.jl new file mode 100644 index 000000000..75f8fbcf0 --- /dev/null +++ b/src/conceptors.jl @@ -0,0 +1,1101 @@ +@doc raw""" + correlation_matrix(states) -> Matrix + +State correlation matrix ``R = X X^\top / L`` for a reservoir state collection +`states` of size ``N \times L`` (one state per column). This is the matrix whose +regularized-identity map defines a conceptor. + +# Arguments + +- `states::AbstractMatrix{<:Real}`: Reservoir states arranged as + `(reservoir_dimension, sample_count)`. + +# Returns + +- A dense square matrix with the floating-point element type inferred from `states`. + +# Throws + +- `ArgumentError`: If `states` has no elements. + +# Example + +```julia +states = rand(Float32, 50, 1_000) +correlation = correlation_matrix(states) +``` +""" +function correlation_matrix(states::AbstractMatrix{<:Real}) + isempty(states) && throw(ArgumentError("states must contain at least one state")) + element_type = float(eltype(states)) + state_matrix = element_type.(states) + return (state_matrix * state_matrix') / size(state_matrix, 2) +end + +@doc raw""" + conceptor_matrix(correlation, aperture) -> Matrix + +Conceptor derived from a correlation matrix and the given `aperture` using +Jaeger (2014), Equation 7. The result is +symmetric positive semidefinite with singular values in ``[0, 1)``. + +# Arguments + +- `correlation::AbstractMatrix{<:Real}`: Symmetric state correlation matrix. +- `aperture::Real`: Finite positive aperture. Larger values admit more state-space + directions; smaller values suppress more directions. + +# Returns + +- A dense conceptor matrix with the floating-point element type of `correlation`. + +# Throws + +- `ArgumentError`: If `aperture` is not finite and positive, or if `correlation` + is not symmetric. +- `DimensionMismatch`: If `correlation` is not square. + +# Example + +```julia +correlation = correlation_matrix(states) +conceptor = conceptor_matrix(correlation, 10) +``` +""" +function conceptor_matrix(correlation::AbstractMatrix{<:Real}, aperture::Real) + isfinite(aperture) && aperture > 0 || + throw(ArgumentError("aperture must be finite and positive, got $aperture")) + checksquare(correlation) + issymmetric(correlation) || throw(ArgumentError("correlation must be symmetric")) + element_type = float(eltype(correlation)) + correlation_matrix = element_type.(correlation) + typed_aperture = element_type(aperture) + inverse_aperture_squared = inv(typed_aperture * typed_aperture) + conceptor = if element_type <: Union{Float32, Float64} + decomposition = eigen(Symmetric(correlation_matrix)) + values = max.(decomposition.values, zero(element_type)) + decomposition.vectors * + Diagonal(values ./ (values .+ inverse_aperture_squared)) * + decomposition.vectors' + else + correlation_matrix / + Symmetric(correlation_matrix + inverse_aperture_squared * I) + end + return (conceptor + conceptor') / element_type(2) +end + +""" + conceptor_from_states(states, aperture) -> Matrix + +Convenience composition of [`correlation_matrix`](@ref) and +[`conceptor_matrix`](@ref): the conceptor characterizing a reservoir state cloud +`states` (size `N × L`) at the given `aperture`. + +# Arguments + +- `states::AbstractMatrix{<:Real}`: Reservoir states, one state per column. +- `aperture::Real`: Finite positive aperture. + +# Returns + +- The conceptor associated with the sample correlation of `states`. + +# Example + +```julia +conceptor = conceptor_from_states(states, 10) +``` +""" +function conceptor_from_states(states::AbstractMatrix{<:Real}, aperture::Real) + return conceptor_matrix(correlation_matrix(states), aperture) +end + +""" + conceptor_singular_values(conceptor) -> Vector + +Singular values of a (symmetric) conceptor matrix `C`, in descending order. For a +conceptor these coincide with its eigenvalues and lie in `[0, 1]`. + +# Arguments + +- `conceptor::AbstractMatrix{<:Real}`: Square conceptor matrix. + +# Returns + +- A vector of singular values in descending order. + +# Throws + +- `DimensionMismatch`: If `conceptor` is not square. +""" +function conceptor_singular_values(conceptor::AbstractMatrix{<:Real}) + checksquare(conceptor) + return svdvals(float.(conceptor)) +end + +""" + quota(conceptor) -> Real + +Mean singular value of a conceptor, `tr(C) / N`. Jaeger's "quota" measures the +fraction of reservoir state space the conceptor leaves open (0 = point, 1 = all). + +# Arguments + +- `conceptor::AbstractMatrix{<:Real}`: Square conceptor matrix. + +# Returns + +- The mean singular value. Values near zero indicate a restrictive conceptor; + values near one indicate a permissive conceptor. +""" +function quota(conceptor::AbstractMatrix{<:Real}) + checksquare(conceptor) + return tr(conceptor) / size(conceptor, 1) +end + +# ====================================================================== +# Aperture adaptation and aperture selection +# ====================================================================== + +@doc raw""" + adapt_singular_value(singular_value, aperture_factor) -> Real + +Aperture-adapted value for an input `singular_value` in `[0, 1]` and a +nonnegative `aperture_factor` (Jaeger 2014, Proposition 3, Equation 19). +The boundary values zero and one are fixed points. + +# Arguments + +- `singular_value::Real`: Value to adapt, normally in `[0, 1]`. +- `aperture_factor::Real`: Nonnegative factor applied to the aperture. + +# Returns + +- The adapted singular value, using the floating-point type inferred from + `singular_value`. + +# Throws + +- `ArgumentError`: If `aperture_factor` is negative. +""" +function adapt_singular_value(singular_value::Real, aperture_factor::Real) + aperture_factor >= 0 || + throw(ArgumentError("aperture_factor must be nonnegative, got $aperture_factor")) + element_type = float(typeof(singular_value)) + typed_value = element_type(singular_value) + typed_factor = element_type(aperture_factor) + typed_value <= zero(element_type) && return zero(element_type) + typed_value >= one(element_type) && return one(element_type) + iszero(typed_factor) && return zero(element_type) + isinf(typed_factor) && return one(element_type) + return typed_value / + (typed_value + (one(element_type) - typed_value) / typed_factor^2) +end + +@doc raw""" + aperture_adapt(conceptor, aperture_factor) -> Matrix + +Aperture adaptation of a `conceptor` by a nonnegative `aperture_factor` +(Jaeger 2014, Definition 3). This is equivalent to multiplying its aperture by +the supplied factor. Zero and infinite factors are handled exactly. + +# Arguments + +- `conceptor::AbstractMatrix{<:Real}`: Symmetric square conceptor matrix. +- `aperture_factor::Real`: Nonnegative aperture multiplier. + +# Returns + +- A conceptor with the same principal directions and adapted singular values. + +# Throws + +- `ArgumentError`: If the factor is negative or `conceptor` is not symmetric. +- `DimensionMismatch`: If `conceptor` is not square. + +# Example + +```julia +wider_conceptor = aperture_adapt(conceptor, 2) +projector = aperture_adapt(conceptor, Inf) +``` +""" +function aperture_adapt(conceptor::AbstractMatrix{<:Real}, aperture_factor::Real) + aperture_factor >= 0 || + throw(ArgumentError("aperture_factor must be nonnegative, got $aperture_factor")) + checksquare(conceptor) + issymmetric(conceptor) || throw(ArgumentError("conceptor must be symmetric")) + element_type = float(eltype(conceptor)) + conceptor_matrix = element_type.(conceptor) + typed_factor = element_type(aperture_factor) + if isfinite(typed_factor) && !iszero(typed_factor) + inverse_factor_squared = inv(typed_factor * typed_factor) + adapted = conceptor_matrix / + Symmetric(conceptor_matrix + inverse_factor_squared * (I - conceptor_matrix)) + return (adapted + adapted') / element_type(2) + end + decomposition = eigen(Symmetric(conceptor_matrix)) + adapted_values = adapt_singular_value.( + clamp.(decomposition.values, zero(element_type), one(element_type)), typed_factor + ) + return decomposition.vectors * Diagonal(adapted_values) * decomposition.vectors' +end + +""" + reaperture(conceptor, from_aperture, to_aperture) -> Matrix + +Re-express a conceptor `C` that currently has aperture `from_aperture` so that it +has aperture `to_aperture` (Jaeger 2014, Proposition 5). + +# Arguments + +- `conceptor::AbstractMatrix{<:Real}`: Conceptor formed at `from_aperture`. +- `from_aperture::Real`: Finite positive current aperture. +- `to_aperture::Real`: Nonnegative target aperture; `Inf` is allowed. + +# Returns + +- The aperture-adapted conceptor. +""" +function reaperture( + conceptor::AbstractMatrix{<:Real}, from_aperture::Real, to_aperture::Real + ) + isfinite(from_aperture) && from_aperture > 0 || + throw(ArgumentError("from_aperture must be finite and positive")) + to_aperture >= 0 || throw(ArgumentError("to_aperture must be nonnegative")) + return aperture_adapt(conceptor, to_aperture / from_aperture) +end + +@doc raw""" + attenuation(conceptor, recurrent_weights, bias; steps, washout, init_state, rng) -> Real + +Attenuation ``a_C = E[\|z(n) - x(n)\|^2] / E[\|z(n)\|^2]`` of a loaded reservoir +(recurrent weights `W`, bias `b`) run autonomously under conceptor `C` +(Jaeger 2014, Eq. 23), where ``z(n) = \tanh(W x(n-1) + b)`` is the unconstrained +update and ``x(n) = C z(n)`` is the conceptor-constrained state. The attenuation +is the fraction of reservoir signal energy suppressed by `C`; as a function of +aperture it passes through a minimum at the best-reconstructing aperture. + +# Arguments + +- `conceptor`: Square conceptor matching the recurrent reservoir dimension. +- `recurrent_weights`: Square autonomous recurrent weight matrix. +- `bias`: Reservoir bias vector. + +# Keywords + +- `steps::Int = 500`: Number of samples included in the energy ratio. +- `washout::Int = 200`: Initial autonomous steps excluded from the ratio. +- `init_state = nothing`: Optional initial reservoir state. +- `rng = Random.default_rng()`: Random number generator used when `init_state` is + omitted. + +# Returns + +- The fraction of unconstrained reservoir energy suppressed by the conceptor. +""" +function attenuation( + conceptor::AbstractMatrix{<:Real}, + recurrent_weights::AbstractMatrix{<:Real}, + bias::AbstractVector{<:Real}; + steps::Int = 500, + washout::Int = 200, + init_state::Union{Nothing, AbstractVector} = nothing, + rng::AbstractRNG = Random.default_rng(), + ) + checksquare(conceptor) + checksquare(recurrent_weights) + size(conceptor) == size(recurrent_weights) || + throw(DimensionMismatch("conceptor and recurrent_weights must have equal dimensions")) + size(recurrent_weights, 1) == length(bias) || + throw(DimensionMismatch("bias length must equal the recurrent matrix dimension")) + steps > 0 || throw(ArgumentError("steps must be positive")) + washout >= 0 || throw(ArgumentError("washout must be nonnegative")) + element_type = float( + promote_type(eltype(conceptor), eltype(recurrent_weights), eltype(bias)) + ) + conceptor_matrix = element_type.(conceptor) + recurrent_matrix = element_type.(recurrent_weights) + bias_vector = element_type.(bias) + reservoir_dimension = size(recurrent_matrix, 1) + state = init_state === nothing ? + element_type(0.5) .* randn(rng, element_type, reservoir_dimension) : + element_type.(init_state) + length(state) == reservoir_dimension || + throw(DimensionMismatch("init_state must have length $reservoir_dimension")) + + suppressed_energy = zero(element_type) + total_energy = zero(element_type) + for step in 1:(washout + steps) + unconstrained_state = tanh.(recurrent_matrix * state .+ bias_vector) + state = conceptor_matrix * unconstrained_state + if step > washout + suppressed_energy += sum(abs2, unconstrained_state .- state) + total_energy += sum(abs2, unconstrained_state) + end + end + iszero(total_energy) && + throw(ArgumentError("attenuation is undefined for a zero-energy rollout")) + return suppressed_energy / total_energy +end + +""" + optimal_aperture(correlation, recurrent_weights, bias, apertures; kwargs...) + +Select the aperture minimizing the [`attenuation`](@ref) criterion over a grid +`apertures`. For each candidate aperture, its conceptor is formed and the +loaded reservoir (`W`, `b`) is run autonomously to measure its attenuation. Returns +the minimizing aperture together with the full vector of attenuations (aligned with +`apertures`) so the characteristic trough can be inspected. `kwargs` are forwarded +to [`attenuation`](@ref). + +# Arguments + +- `correlation`: Symmetric state correlation matrix. +- `recurrent_weights`: Loaded recurrent weight matrix. +- `bias`: Reservoir bias vector. +- `apertures`: Nonempty vector of finite positive candidate apertures. + +# Returns + +- `(best_aperture, attenuations)`, where `attenuations` follows the order of + `apertures`. +""" +function optimal_aperture( + correlation::AbstractMatrix{<:Real}, + recurrent_weights::AbstractMatrix{<:Real}, + bias::AbstractVector{<:Real}, + apertures::AbstractVector{<:Real}; + kwargs..., + ) + isempty(apertures) && throw(ArgumentError("apertures must not be empty")) + attenuations = [ + attenuation( + conceptor_matrix(correlation, aperture), recurrent_weights, bias; kwargs... + ) for aperture in apertures + ] + return apertures[argmin(attenuations)], attenuations +end + +# ====================================================================== +# Boolean operations on conceptors (Jaeger 2014, Section 3.9) +# ====================================================================== + +@doc raw""" + conceptor_not(conceptor) -> Matrix + +Logical negation of a conceptor (Jaeger 2014, Equation 28). +Exchanges the roles of the directions the conceptor admits and suppresses. + +# Arguments + +- `conceptor::AbstractMatrix{<:Real}`: Square conceptor matrix. + +# Returns + +- The complementary conceptor. Its singular values are one minus those of + `conceptor`. +""" +function conceptor_not(conceptor::AbstractMatrix{<:Real}) + checksquare(conceptor) + element_type = float(eltype(conceptor)) + conceptor_matrix = element_type.(conceptor) + return Matrix{element_type}(I, size(conceptor_matrix)) - conceptor_matrix +end + +@doc raw""" + conceptor_and(first_conceptor, second_conceptor) -> Matrix + +Logical conjunction of two conceptors (Jaeger 2014, Equation 32). For +full-rank conceptors this equals ``(C^{-1} + B^{-1} - I)^{-1}``; the implementation +uses the range-intersection projector form so that singular conceptors are handled +robustly. The result admits exactly the reservoir directions admitted by *both* +`C` and `B`. + +# Arguments + +- `first_conceptor::AbstractMatrix{<:Real}`: First symmetric conceptor. +- `second_conceptor::AbstractMatrix{<:Real}`: Second symmetric conceptor of the + same size. + +# Returns + +- The conjunction of the two conceptors, including when either input is + rank-deficient. + +# Throws + +- `DimensionMismatch`: If the conceptors are not square or have different sizes. +- `ArgumentError`: If either conceptor is not symmetric. +""" +function conceptor_and( + first_conceptor::AbstractMatrix{<:Real}, + second_conceptor::AbstractMatrix{<:Real} + ) + checksquare(first_conceptor) + checksquare(second_conceptor) + element_type = float(promote_type(eltype(first_conceptor), eltype(second_conceptor))) + first_matrix = element_type.(first_conceptor) + second_matrix = element_type.(second_conceptor) + size(first_matrix) == size(second_matrix) || + throw(DimensionMismatch("conceptors must have equal size")) + issymmetric(first_matrix) || throw(ArgumentError("first_conceptor must be symmetric")) + issymmetric(second_matrix) || throw(ArgumentError("second_conceptor must be symmetric")) + + # Basis of R(C) ∩ R(B) = N(P_{N(C)} + P_{N(B)}), via the null-space projectors. + first_nullspace = nullspace(first_matrix) + second_nullspace = nullspace(second_matrix) + nullspace_projector = + first_nullspace * first_nullspace' + second_nullspace * second_nullspace' + intersection_basis = nullspace(nullspace_projector) + + if size(intersection_basis, 2) == 0 + return zeros(element_type, size(first_matrix)) + end + core = intersection_basis' * (pinv(first_matrix) + pinv(second_matrix) - I) * + intersection_basis + result = intersection_basis * + (core \ Matrix{element_type}(I, size(core))) * + intersection_basis' + return (result + result') / element_type(2) +end + +@doc raw""" + conceptor_or(first_conceptor, second_conceptor) -> Matrix + +Logical disjunction of two conceptors (Jaeger 2014, Equation 31), defined +via De Morgan's law. The result admits +every reservoir direction admitted by `C` or by `B`. + +# Arguments + +- `first_conceptor::AbstractMatrix{<:Real}`: First symmetric conceptor. +- `second_conceptor::AbstractMatrix{<:Real}`: Second symmetric conceptor of the + same size. + +# Returns + +- The disjunction of the two conceptors. +""" +function conceptor_or( + first_conceptor::AbstractMatrix{<:Real}, + second_conceptor::AbstractMatrix{<:Real} + ) + return conceptor_not( + conceptor_and(conceptor_not(first_conceptor), conceptor_not(second_conceptor)) + ) +end + +# ====================================================================== +# The Conceptor reservoir-computer wrapper and its parameter/state plumbing +# ====================================================================== + +@doc raw""" + Conceptor(model) + +Wrap a ReservoirComputing `model` (e.g. an [`ESN`](@ref)) so that conceptor +matrices can be derived from it, stored by name, combined with the conceptor +algebra, and used to constrain autonomous generation (Jaeger 2014). + +The wrapped `model` provides the reservoir update; the `Conceptor` adds a +conceptor library to the model state. Parameters and states are nested under the +`model` field, leaving room for conceptor bookkeeping at the top level. + +# Arguments + +- `model`: Reservoir computer to wrap. Pattern loading expects an ESN-compatible + parameter layout with reservoir, input, bias, and readout weights. + +# Returns + +- A conceptor-enabled reservoir computer used with [`initialparameters`](@ref), + [`initialstates`](@ref), [`load!`](@ref), and [`generate`](@ref). + +# Example + +```julia +model = Conceptor(ESN(1, 100, 1; use_bias = true)) +parameters = initialparameters(rng, model) +states = initialstates(rng, model) +``` +""" +@concrete struct Conceptor <: AbstractReservoirComputer{(:model,)} + model::Any +end + +function initialparameters(rng::AbstractRNG, concept::Conceptor) + return (; model = initialparameters(rng, concept.model), readout = nothing) +end + +function initialstates(rng::AbstractRNG, concept::Conceptor) + return (; + model = initialstates(rng, concept.model), + conceptors = Dict{Symbol, Matrix}(), + apertures = Dict{Symbol, Real}(), + active_conceptor = nothing, + ) +end + +""" + has_conceptor(st, name) -> Bool + +Whether a conceptor called `name` has been stored in the state `st`. + +Returns `true` when `name` is present and `false` otherwise. +""" +has_conceptor(st::NamedTuple, name::Symbol) = haskey(st.conceptors, name) + +""" + get_conceptor(st, name) -> Matrix or nothing + +The stored conceptor matrix called `name`, or `nothing` if it is absent. + +Use [`has_conceptor`](@ref) when absence should be checked separately. +""" +get_conceptor(st::NamedTuple, name::Symbol) = get(st.conceptors, name, nothing) + +""" + store_conceptor!(st, name, conceptor, aperture) -> st + +Record conceptor matrix `C` (formed at `aperture`) under `name` in the state's +conceptor library. Mutates the library dictionaries in place and returns `st`. + +# Arguments + +- `st::NamedTuple`: State returned by `initialstates` for a [`Conceptor`](@ref). +- `name::Symbol`: Name used to retrieve the conceptor later. +- `conceptor::AbstractMatrix{<:Real}`: Square conceptor matrix. +- `aperture::Real`: Finite positive aperture associated with the matrix. + +# Returns + +- The same state object, with its conceptor and aperture dictionaries updated. +""" +function store_conceptor!( + st::NamedTuple, + name::Symbol, + conceptor::AbstractMatrix{<:Real}, + aperture::Real, + ) + checksquare(conceptor) + isfinite(aperture) && aperture > 0 || + throw(ArgumentError("aperture must be finite and positive, got $aperture")) + st.conceptors[name] = Matrix(conceptor) + st.apertures[name] = convert(float(eltype(conceptor)), aperture) + return st +end + +""" + set_active_conceptor(st, name) -> st + +Mark the conceptor `name` (or `nothing`) as the active one, returning an updated +state. The active conceptor is the default constraint for generation. + +Passing `nothing` clears the active selection. A `KeyError` is thrown when a +nonexistent name is selected. +""" +function set_active_conceptor(st::NamedTuple, name::Union{Symbol, Nothing}) + name !== nothing && !has_conceptor(st, name) && throw(KeyError(name)) + return merge(st, (; active_conceptor = name)) +end + +""" + active_conceptor(st) -> Matrix + +The currently active conceptor matrix. Throws if none has been set. + +Use [`set_active_conceptor`](@ref) to select or clear the active conceptor. +""" +function active_conceptor(st::NamedTuple) + name = st.active_conceptor + name === nothing && + throw(ArgumentError("no active conceptor; call set_active_conceptor")) + return st.conceptors[name] +end + +# Resolve a conceptor argument that is either an explicit matrix or a stored name. +resolve_conceptor(::NamedTuple, conceptor::AbstractMatrix{<:Real}) = conceptor +function resolve_conceptor(st::NamedTuple, conceptor::Symbol) + stored_conceptor = get_conceptor(st, conceptor) + stored_conceptor === nothing && throw(KeyError(conceptor)) + return stored_conceptor +end + +""" + collectstates(concept::Conceptor, signal, ps, st; conceptor=nothing) -> (states, st) + +Stream `signal` through the wrapped reservoir and collect its state sequence. When +`conceptor` (a stored `Symbol` or an explicit matrix) is given, the collected +states are filtered through it, `C x(n)`, realizing the conceptor as a state +modifier. +""" +function collectstates( + concept::Conceptor, + signal::AbstractMatrix, + ps::NamedTuple, + st::NamedTuple; + conceptor::Union{Symbol, AbstractMatrix, Nothing} = nothing, + ) + states, st_model = collectstates(concept.model, signal, ps.model, st.model) + new_st = merge(st, (; model = st_model)) + conceptor === nothing && return states, new_st + return resolve_conceptor(new_st, conceptor) * states, new_st +end + +""" + addreadout!(concept::Conceptor, weights, ps, st) -> (ps, st) + +Install trained readout weights `W` into the wrapped model's readout layer +(`ps.model.readout.weight`), keeping the generation and classification paths on a +single readout convention. +""" +function addreadout!( + concept::Conceptor, weights::AbstractMatrix, ps::NamedTuple, st::NamedTuple + ) + return set_readout_weight(ps, weights), st +end + +reservoir_params(ps::NamedTuple) = ps.model.reservoir +readout_params(ps::NamedTuple) = ps.model.readout + +function reservoir_bias(ps::NamedTuple) + rps = reservoir_params(ps) + element_type = eltype(rps.reservoir_matrix) + return haskey(rps, :bias) ? + vec(rps.bias) : zeros(element_type, size(rps.reservoir_matrix, 1)) +end + +function set_reservoir_matrix(ps::NamedTuple, weights::AbstractMatrix) + reservoir = merge(ps.model.reservoir, (; reservoir_matrix = weights)) + return merge(ps, (; model = merge(ps.model, (; reservoir = reservoir)))) +end + +function set_readout_weight(ps::NamedTuple, weights::AbstractMatrix) + readout = merge(ps.model.readout, (; weight = weights)) + return merge(ps, (; model = merge(ps.model, (; readout = readout)))) +end + +@doc raw""" + ridge_map(features, targets, reg) -> Matrix + +Tikhonov-regularized linear map `M` minimizing +``\|M \,\text{features} - \text{targets}\|^2 + \text{reg}\,\|M\|^2``. `features` is +``n_\text{in} \times T``, `targets` is ``n_\text{out} \times T``, and the returned +`M` has size ``n_\text{out} \times n_\text{in}``. + +# Arguments + +- `features`: Feature matrix with observations in columns. +- `targets`: Target matrix with observations in columns. +- `reg::Real`: Nonnegative ridge penalty. + +# Returns + +- A matrix mapping feature columns to target columns. +""" +function ridge_map( + features::AbstractMatrix{<:Real}, + targets::AbstractMatrix{<:Real}, + reg::Real, + ) + element_type = float(promote_type(eltype(features), eltype(targets))) + feature_matrix = element_type.(features) + target_matrix = element_type.(targets) + return train(StandardRidge(element_type, reg), feature_matrix, target_matrix) +end + +# Per-pattern aperture lookup: a scalar applies to every pattern; a dictionary +# supplies a value per name. +aperture_for(aperture::Real, ::Symbol) = aperture +function aperture_for(aperture::AbstractDict{Symbol, <:Real}, name::Symbol) + haskey(aperture, name) || throw(ArgumentError("no aperture provided for :$name")) + return aperture[name] +end + +# loading patterns + +function _drive_pattern( + rng::AbstractRNG, + concept::Conceptor, + signal::AbstractMatrix, + ps::NamedTuple, + st::NamedTuple, + washout::Int, + ) + st_reset = resetcarry!(rng, concept.model, st.model; init_carry = nothing) + states, st_model = collectstates(concept.model, signal, ps.model, st_reset) + signal_length = size(states, 2) + washout < signal_length || + throw(ArgumentError("washout=$washout exceeds signal length=$signal_length")) + current_states = states[:, (washout + 1):signal_length] + lagged_states = states[:, washout:(signal_length - 1)] + aligned_signal = signal[:, (washout + 1):signal_length] + return current_states, lagged_states, aligned_signal, merge(st, (; model = st_model)) +end + +@doc raw""" + load!(rng, concept, named_signals, ps, st; + aperture=10.0, washout=500, reg_recurrent=1e-4, reg_readout=1e-2) + +Load named driving patterns into the reservoir wrapped by `concept` +(Jaeger 2014, Section 3.2). For each `name => signal` pair the reservoir is run +from a cleared carry, the first `washout` states are discarded, and the remaining +states are used to + + * derive and store a conceptor from the state correlation matrix, and + * accumulate data for a single shared input-internalizing recurrent matrix `W` + and readout `W_out`. + +`W` is fitted (ridge, `reg_recurrent`) so that +``W x(n-1) \approx W^* x(n-1) + W_\text{in} p(n)``, absorbing the input drive into +the recurrent weights; `W_out` is fitted (ridge, `reg_readout`) so that +``W_\text{out} x(n) \approx p(n)`` across all patterns. `aperture` is either a +scalar or a dictionary mapping each name to an aperture. + +Returns `(ps, st)` with `reservoir_matrix` and `readout.weight` replaced and the +conceptor library populated. + +# Arguments + +- `rng::AbstractRNG`: Random number generator used when reservoir carries reset. +- `concept::Conceptor`: Conceptor-wrapped reservoir. +- `named_signals`: Iterable of `Symbol => signal` pairs. A vector signal is treated + as one-dimensional; a matrix has one input channel per row and time per column. +- `ps::NamedTuple`: Current parameters. +- `st::NamedTuple`: Current states. + +# Keywords + +- `aperture = 10.0`: One aperture for every pattern, or a dictionary keyed by + pattern name. +- `washout::Int = 500`: Initial samples excluded from fitting. +- `reg_recurrent::Real = 1.0e-4`: Ridge penalty for recurrent weights. +- `reg_readout::Real = 1.0e-2`: Ridge penalty for readout weights. + +# Returns + +- `(parameters, states)`: Updated parameters and a state containing one conceptor + per named signal. + +# Throws + +- `ArgumentError`: If a name is not a `Symbol`, an aperture is missing, or washout + is not shorter than a signal. +""" +function load!( + rng::AbstractRNG, + concept::Conceptor, + named_signals, + ps::NamedTuple, + st::NamedTuple; + aperture = 10.0, + washout::Int = 500, + reg_recurrent::Real = 1.0e-4, + reg_readout::Real = 1.0e-2, + ) + rps = reservoir_params(ps) + element_type = float(eltype(rps.reservoir_matrix)) + original_recurrent_weights = element_type.(rps.reservoir_matrix) + input_weights = element_type.(rps.input_matrix) + + state_sets = Matrix{element_type}[] + lagged_state_sets = Matrix{element_type}[] + signal_sets = Matrix{element_type}[] + st_acc = st + + for (name, signal) in named_signals + name isa Symbol || + throw(ArgumentError("signal name must be a Symbol, got $(typeof(name))")) + signal_matrix = signal isa AbstractMatrix ? + element_type.(signal) : reshape(element_type.(signal), 1, :) + current_states, lagged_states, aligned_signal, st_acc = + _drive_pattern(rng, concept, signal_matrix, ps, st_acc, washout) + pattern_aperture = aperture_for(aperture, name) + store_conceptor!( + st_acc, name, conceptor_from_states(current_states, pattern_aperture), + pattern_aperture + ) + push!(state_sets, current_states) + push!(lagged_state_sets, lagged_states) + push!(signal_sets, aligned_signal) + end + + all_states = reduce(hcat, state_sets) + all_lagged_states = reduce(hcat, lagged_state_sets) + all_signals = reduce(hcat, signal_sets) + + readout_weights = ridge_map(all_states, all_signals, reg_readout) + recurrent_target = + original_recurrent_weights * all_lagged_states .+ input_weights * all_signals + recurrent_weights = ridge_map(all_lagged_states, recurrent_target, reg_recurrent) + + updated_parameters = set_reservoir_matrix(ps, element_type.(recurrent_weights)) + updated_parameters = + set_readout_weight(updated_parameters, element_type.(readout_weights)) + return updated_parameters, st_acc +end + +# ====================================================================== +# Autonomous pattern generation (Jaeger 2014, Section 3.2) +# ====================================================================== + +@doc raw""" + generate(concept, ps, st; conceptor, steps, washout=200, + init_state=nothing, rng=Random.default_rng()) -> (Y, X) + +Run the loaded reservoir autonomously under `conceptor` and read out the observer +signal: + +```math +x(n) = C \tanh(W x(n-1) + b), \qquad y(n) = W_\text{out} x(n). +``` + +`conceptor` is either a stored conceptor `Symbol` or an explicit conceptor matrix +(e.g. the output of [`morph_conceptor`](@ref) or the Boolean operations). The first +`washout` steps let the autonomous orbit settle and are discarded. Returns `(Y, X)`: +the post-washout observer outputs (`out_dims × steps`) and reservoir states +(`res_dims × steps`). The rollout uses the element type of the reservoir parameters. + +# Arguments + +- `concept::Conceptor`: Loaded conceptor model. +- `ps::NamedTuple`: Parameters returned by [`load!`](@ref). +- `st::NamedTuple`: State containing the requested conceptor. + +# Keywords + +- `conceptor`: Stored conceptor name or an explicit square conceptor matrix. +- `steps::Int`: Number of returned time steps; must be positive. +- `washout::Int = 200`: Autonomous steps discarded before recording. +- `init_state = nothing`: Optional reservoir state vector. A random state is used + when omitted. +- `rng = Random.default_rng()`: Random number generator for the initial state. + +# Returns + +- `(outputs, states)`: Observer outputs and reservoir states, with time along + columns. + +# Throws + +- `KeyError`: If a requested conceptor name is absent. +- `DimensionMismatch`: If the conceptor or initial state does not match the + reservoir dimension. +- `ArgumentError`: If `steps` is not positive or `washout` is negative. + +# Example + +```julia +outputs, states = generate( + concept, parameters, model_state; + conceptor = :sine, steps = 500, washout = 100, rng +) +``` +""" +function generate( + concept::Conceptor, + ps::NamedTuple, + st::NamedTuple; + conceptor::Union{Symbol, AbstractMatrix}, + steps::Int, + washout::Int = 200, + init_state::Union{Nothing, AbstractVector} = nothing, + rng::AbstractRNG = Random.default_rng(), + ) + rps = reservoir_params(ps) + element_type = float(eltype(rps.reservoir_matrix)) + recurrent_weights = element_type.(rps.reservoir_matrix) + bias = reservoir_bias(ps) + readout_weights = element_type.(readout_params(ps).weight) + conceptor_matrix = element_type.(resolve_conceptor(st, conceptor)) + reservoir_dimension = size(recurrent_weights, 1) + checksquare(conceptor_matrix) + size(conceptor_matrix, 1) == reservoir_dimension || + throw( + DimensionMismatch( + "conceptor size must match the reservoir dimension $reservoir_dimension" + ) + ) + steps > 0 || throw(ArgumentError("steps must be positive")) + washout >= 0 || throw(ArgumentError("washout must be nonnegative")) + + state = init_state === nothing ? + element_type(0.5) .* randn(rng, element_type, reservoir_dimension) : + element_type.(init_state) + length(state) == reservoir_dimension || + throw(DimensionMismatch("init_state must have length $reservoir_dimension")) + output_dimension = size(readout_weights, 1) + outputs = zeros(element_type, output_dimension, steps) + states = zeros(element_type, reservoir_dimension, steps) + + for step in 1:(washout + steps) + state = conceptor_matrix * tanh.(recurrent_weights * state .+ bias) + if step > washout + output_index = step - washout + @views states[:, output_index] = state + @views outputs[:, output_index] = readout_weights * state + end + end + return outputs, states +end + +# morphing + +@doc raw""" + morph_conceptor(st, weights) -> Matrix + +Linearly combine stored conceptors into a morphed conceptor +(Jaeger 2014). `weights` is an iterable of name-to-weight pairs, a `NamedTuple`, +or a `Dict{Symbol,<:Real}`. Coefficients summing to one +interpolate between the named prototypes; coefficients outside `[0, 1]` extrapolate. + +# Arguments + +- `st::NamedTuple`: State containing stored conceptors. +- `weights`: Named tuple, dictionary, or iterable of name-to-weight pairs. + +# Returns + +- The weighted sum of the named conceptors. + +# Throws + +- `KeyError`: If a named conceptor is absent. +- `ArgumentError`: If `weights` is empty. +""" +function morph_conceptor(st::NamedTuple, weights) + morphed_conceptor = nothing + for (name, weight) in pairs(weights) + stored_conceptor = get_conceptor(st, Symbol(name)) + stored_conceptor === nothing && throw(KeyError(Symbol(name))) + contribution = convert(eltype(stored_conceptor), weight) .* stored_conceptor + morphed_conceptor = morphed_conceptor === nothing ? + contribution : morphed_conceptor .+ contribution + end + morphed_conceptor === nothing && throw(ArgumentError("no weights provided")) + return morphed_conceptor +end + +# ====================================================================== +# Conceptors as state filters for supervised readout training +# ====================================================================== + +# Normalize a target into a row-major matrix (1 × T for a vector target). +_as_matrix(target::AbstractMatrix) = target +_as_matrix(target::AbstractVector) = reshape(target, 1, :) + +""" + store_conceptors!(rng, concept, named_signals, ps, st; + aperture=1.0, init_carry=nothing) -> st + +Derive and store one conceptor per named signal without modifying the reservoir +weights. For each `name => signal` the reservoir is run from a cleared carry and +the conceptor of its state cloud is stored. Use +this to build a conceptor library for classification, where each class pattern +gets its own conceptor. `aperture` is a scalar or a dictionary keyed by name. + +# Returns + +- The updated state. Reservoir parameters are unchanged. + +# See also + +[`load!`](@ref), [`train!`](@ref), [`store_conceptor!`](@ref) +""" +function store_conceptors!( + rng::AbstractRNG, + concept::Conceptor, + named_signals, + ps::NamedTuple, + st::NamedTuple; + aperture = 1.0, + init_carry = nothing, + ) + st_acc = st + for (name, signal) in named_signals + name isa Symbol || + throw(ArgumentError("signal name must be a Symbol, got $(typeof(name))")) + sig = signal isa AbstractMatrix ? signal : reshape(signal, 1, :) + st_reset = resetcarry!(rng, concept.model, st_acc.model; init_carry = init_carry) + states, st_model = collectstates(concept.model, sig, ps.model, st_reset) + st_acc = merge(st_acc, (; model = st_model)) + pattern_aperture = aperture_for(aperture, name) + store_conceptor!( + st_acc, name, conceptor_from_states(states, pattern_aperture), pattern_aperture + ) + end + return st_acc +end + +""" + train!(rng, concept::Conceptor, named_signals, named_targets, ps, st, + train_method=StandardRidge(0.0); + washout=0, return_states=false, init_carry=nothing) -> (ps, st) + +Train a single readout on conceptor-filtered reservoir features. For each +`name => signal` the states are collected through the conceptor named `name` +(stored beforehand via [`store_conceptors!`](@ref)), the first `washout` columns +are dropped, and features and the matching `named_targets` are pooled across +patterns before the readout is fit with `train_method`. With `return_states=true` +the pooled feature matrix is also returned. + +# Arguments + +- `named_signals`: Iterable of named input sequences. +- `named_targets`: Iterable containing one target sequence for every signal name. +- `train_method`: Readout fitting method; defaults to unregularized + [`StandardRidge`](@ref). + +# Keywords + +- `washout::Int = 0`: Initial columns removed from every feature and target. +- `return_states::Bool = false`: Return pooled conceptor-filtered features. +- `init_carry = nothing`: Initial carry passed when resetting each sequence. +- `kwargs...`: Additional options forwarded to `train`. + +# Returns + +- `(parameters, states)` by default. +- `((parameters, states), features)` when `return_states = true`. +""" +function train!( + rng::AbstractRNG, + concept::Conceptor, + named_signals, + named_targets, + ps::NamedTuple, + st::NamedTuple, + train_method = StandardRidge(0.0); + washout::Int = 0, + return_states::Bool = false, + init_carry = nothing, + kwargs..., + ) + targets = Dict{Symbol, AbstractMatrix}( + Symbol(name) => _as_matrix(target) for (name, target) in named_targets + ) + + all_states = AbstractMatrix[] + all_targets = AbstractMatrix[] + for (name, signal) in named_signals + name isa Symbol || + throw(ArgumentError("signal name must be a Symbol, got $(typeof(name))")) + haskey(targets, name) || throw(ArgumentError("no target data for pattern :$name")) + st_reset = resetcarry!(rng, concept.model, st.model; init_carry = init_carry) + st_tmp = merge(st, (; model = st_reset)) + states, _ = collectstates(concept, signal, ps, st_tmp; conceptor = name) + tgt = targets[name] + if washout > 0 + states = states[:, (washout + 1):end] + tgt = tgt[:, (washout + 1):end] + end + push!(all_states, states) + push!(all_targets, tgt) + end + + features = reduce(hcat, all_states) + pooled_targets = reduce(hcat, all_targets) + W_out = train(train_method, features, pooled_targets; kwargs...) + ps2, st2 = addreadout!(concept, W_out, ps, st) + return return_states ? ((ps2, st2), features) : (ps2, st2) +end diff --git a/test/test_conceptors.jl b/test/test_conceptors.jl new file mode 100644 index 000000000..8b2abd92f --- /dev/null +++ b/test/test_conceptors.jl @@ -0,0 +1,191 @@ +using Test +using LinearAlgebra +using Random +using Statistics +using ReservoirComputing + +# A small reproducible reservoir state cloud for algebra tests. +function sample_correlation(rng, state_dimension, sample_count) + states = randn(rng, state_dimension, sample_count) + return correlation_matrix(states) +end + +@testset "Conceptors" begin + rng = Xoshiro(2024) + + @testset "conceptor matrix" begin + correlation = sample_correlation(rng, 20, 200) + conceptor = conceptor_matrix(correlation, 10.0) + @test size(conceptor) == (20, 20) + @test conceptor ≈ conceptor' + eigenvalues = eigvals(Symmetric(conceptor)) + @test all(>=(-1.0e-10), eigenvalues) + @test all(<=(1.0 + 1.0e-10), eigenvalues) + @test conceptor_matrix(correlation, 1.0e8) ≈ I(20) atol = 1.0e-3 + @test norm(conceptor_matrix(correlation, 1.0e-6)) < 1.0e-6 + @test_throws ArgumentError conceptor_matrix(correlation, -1.0) + @test_throws ArgumentError conceptor_matrix(correlation, Inf) + @test_throws ArgumentError correlation_matrix(zeros(3, 0)) + + for element_type in (Float32, Float64, BigFloat) + typed_states = element_type.(randn(rng, 6, 30)) + typed_correlation = correlation_matrix(typed_states) + typed_conceptor = conceptor_matrix(typed_correlation, 3) + @test eltype(typed_correlation) == element_type + @test eltype(typed_conceptor) == element_type + @test typed_conceptor ≈ typed_conceptor' + @test eltype(aperture_adapt(typed_conceptor, 2)) == element_type + end + end + + @testset "aperture adaptation" begin + correlation = sample_correlation(rng, 15, 150) + conceptor = conceptor_matrix(correlation, 5.0) + @test aperture_adapt(conceptor, 1.0) ≈ conceptor + # Repeated adaptations multiply their aperture factors. + @test aperture_adapt(aperture_adapt(conceptor, 2.0), 3.0) ≈ + aperture_adapt(conceptor, 6.0) + quotas = [ + quota(aperture_adapt(conceptor, factor)) for + factor in (0.25, 1.0, 4.0, 16.0) + ] + @test issorted(quotas) + projector = aperture_adapt(conceptor, Inf) + @test projector * projector ≈ projector atol = 1.0e-8 + @test norm(aperture_adapt(conceptor, 0.0)) < 1.0e-8 + @test conceptor_matrix(correlation, 7.0) ≈ + aperture_adapt(conceptor_matrix(correlation, 1.0), 7.0) + end + + @testset "Boolean algebra" begin + first_conceptor = conceptor_matrix(sample_correlation(rng, 12, 120), 8.0) + second_conceptor = conceptor_matrix(sample_correlation(rng, 12, 120), 8.0) + third_conceptor = conceptor_matrix(sample_correlation(rng, 12, 120), 8.0) + + @test conceptor_not(conceptor_not(first_conceptor)) ≈ first_conceptor + @test conceptor_not(first_conceptor) ≈ I(12) - first_conceptor + # De Morgan + @test conceptor_or(first_conceptor, second_conceptor) ≈ + conceptor_not( + conceptor_and( + conceptor_not(first_conceptor), conceptor_not(second_conceptor) + ) + ) + @test conceptor_and(first_conceptor, second_conceptor) ≈ + conceptor_not( + conceptor_or( + conceptor_not(first_conceptor), conceptor_not(second_conceptor) + ) + ) + # commutativity + @test conceptor_and(first_conceptor, second_conceptor) ≈ + conceptor_and(second_conceptor, first_conceptor) + @test conceptor_or(first_conceptor, second_conceptor) ≈ + conceptor_or(second_conceptor, first_conceptor) + # AND shrinks, OR grows (quota ordering) + @test quota(conceptor_and(first_conceptor, second_conceptor)) <= + quota(first_conceptor) + 1.0e-8 + @test quota(first_conceptor) <= + quota(conceptor_or(first_conceptor, second_conceptor)) + 1.0e-8 + # associativity (holds for conceptors, unlike idempotence) + @test conceptor_and( + conceptor_and(first_conceptor, second_conceptor), third_conceptor + ) ≈ conceptor_and( + first_conceptor, conceptor_and(second_conceptor, third_conceptor) + ) atol = 1.0e-6 + @test conceptor_or( + conceptor_or(first_conceptor, second_conceptor), third_conceptor + ) ≈ conceptor_or( + first_conceptor, conceptor_or(second_conceptor, third_conceptor) + ) atol = 1.0e-6 + # neutral / absorbing elements + @test conceptor_and(first_conceptor, Matrix(1.0I, 12, 12)) ≈ + first_conceptor atol = 1.0e-6 + @test conceptor_or(first_conceptor, zeros(12, 12)) ≈ + first_conceptor atol = 1.0e-6 + + # The generalized definition must also work for singular conceptors. + first_singular_conceptor = Diagonal([0.8, 0.0, 0.4]) |> Matrix + second_singular_conceptor = Diagonal([0.5, 0.7, 0.0]) |> Matrix + @test conceptor_and(first_singular_conceptor, second_singular_conceptor) ≈ + Diagonal([1 / (1 / 0.8 + 1 / 0.5 - 1), 0, 0]) + end + + @testset "wrapper / library" begin + rng2 = Xoshiro(7) + esn = ESN(1, 30, 1; use_bias = true) + concept = Conceptor(esn) + ps = initialparameters(rng2, concept) + st = initialstates(rng2, concept) + @test haskey(ps, :model) + @test isempty(st.conceptors) + conceptor = conceptor_matrix(sample_correlation(rng2, 30, 100), 10.0) + store_conceptor!(st, :first_pattern, conceptor, 10.0) + @test has_conceptor(st, :first_pattern) + @test get_conceptor(st, :first_pattern) ≈ conceptor + st = set_active_conceptor(st, :first_pattern) + @test active_conceptor(st) ≈ conceptor + @test_throws KeyError set_active_conceptor(st, :missing) + end + + @testset "load! / generate roundtrip" begin + rng3 = Xoshiro(99) + bias_init(random_generator, dims...) = + 0.2f0 .* randn(random_generator, Float32, dims...) + input_init(random_generator, dims...) = + 1.5f0 .* randn(random_generator, Float32, dims...) + res_init(random_generator, dims...) = + rand_sparse(random_generator, dims...; radius = 1.5f0, sparsity = 0.1f0) + esn = ESN( + 1, 80, 1; use_bias = true, init_bias = bias_init, + init_input = input_init, init_reservoir = res_init + ) + concept = Conceptor(esn) + ps = initialparameters(rng3, concept) + st = initialstates(rng3, concept) + + sample_indices = 1:1000 + period = 8.8342522 + signal = reshape(Float32.(sin.(2pi .* sample_indices ./ period)), 1, :) + ps, st = load!( + rng3, concept, [:sine => signal], ps, st; aperture = 1000.0, washout = 400 + ) + @test has_conceptor(st, :sine) + + outputs, states = generate( + concept, ps, st; conceptor = :sine, steps = 200, washout = 300, rng = rng3 + ) + @test size(outputs) == (1, 200) + @test size(states) == (80, 200) + output = vec(outputs) + best = minimum( + sqrt( + mean( + abs2, + output .- [ + sign * sin(2pi * index / period + phase) for + index in eachindex(output) + ], + ) + ) / std(output) for phase in 0:0.02:2pi, sign in (-1.0, 1.0) + ) + @test best < 0.05 + @test count(>(0.5), conceptor_singular_values(get_conceptor(st, :sine))) < 80 + end + + @testset "morph_conceptor" begin + rng4 = Xoshiro(5) + esn = ESN(1, 25, 1; use_bias = true) + concept = Conceptor(esn) + st = initialstates(rng4, concept) + first_conceptor = conceptor_matrix(sample_correlation(rng4, 25, 100), 10.0) + second_conceptor = conceptor_matrix(sample_correlation(rng4, 25, 100), 10.0) + store_conceptor!(st, :first_pattern, first_conceptor, 10.0) + store_conceptor!(st, :second_pattern, second_conceptor, 10.0) + @test morph_conceptor(st, (; first_pattern = 1.0, second_pattern = 0.0)) ≈ + first_conceptor + @test morph_conceptor(st, (; first_pattern = 0.3, second_pattern = 0.7)) ≈ + 0.3 .* first_conceptor .+ 0.7 .* second_conceptor + @test_throws KeyError morph_conceptor(st, (; missing_name = 1.0)) + end +end