#version 300 es

// @secao     8. Superposição — a matriz de Gram
// @titulo    Superposição — quem mais responde quando uma feature acende
// @subtitulo G = WWᵀ · n = 20 features em d dimensões
// @fonte     "A superposição é uma característica, não um bug"
// @doi       10.5281/zenodo.21865370
// @status    EXACT
// @nota      Cada célula (i,j) é o produto interno wᵢ·wⱼ entre duas features, ambos de norma 1. A diagonal é 1 em toda a matriz (a feature responde a si mesma). O que aparece fora dela é a interferência: a feature j respondendo quando a i acende. Reduza d e o bloco cresce.
// @leitura   As linhas seguem a ordem de importância: a linha 0 é a mais importante. Em d ≥ n a matriz é uma diagonal perfeita — uma direção dedicada por feature, leitura trivial. Conforme d cai, os pontos fora da diagonal se multiplicam e se agrupam em blocos. A varredura acende uma linha por vez: tudo o que compartilha a direção dela acende junto. Em d = 1 a matriz inteira fica acesa, porque G ≡ 1: com uma só direção todas as features são o mesmo vetor e nenhuma pode ser distinguida.
// @honesto   Esta matriz NÃO são os pesos de um LLM e não foi aprendida por gradient descent. É uma alocação que eu construí por regra explícita: as primeiras features recebem eₚ dedicado; as demais caem numa direção por hash e se inclinam dentro de um arco, tanto mais aglomeradas quanto menos importantes. Duas ressalvas que a imagem mostra: colisões de hash podem fazer uma feature menos importante cair na direção de uma dedicada, então a "dedicação" das primeiras não é absoluta; e esta construção não é o ótimo de nenhuma perda — é apenas uma geometria plausível. O que é exato é a álgebra: dado W, G = WWᵀ e a leitura linear x̂ = Gx. Que um modelo treinado chegue a essa geometria bloco-semi-ortogonal é o resultado empírico de Elhage et al. (2022) e Scherlis et al. (2022); eu não o reproduzo aqui.
// @custo     A conta é exata: G_ij = wᵢ·wⱼ é o produto interno literal de dois vetores unitários, sem aproximação. O shader avalia essa expressão por célula a cada frame.
// @param1    dimensões d | 1.0 | 20.0 | 8.0
// @param2    decaimento de importância | 0.0 | 3.0 | 1.0

precision highp float;

uniform vec2  uRes;
uniform float uTime;
uniform float uParam1;   // d — quantas direções o espaço oferece
uniform float uParam2;   // expoente do decaimento de importância

out vec4 fragColor;

const int N = 20;        // n — quantas features precisam de representação
const float MEIA_VOLTA = 1.57079632679;

// ─────────────────────────────────────────────────────────────
// O MODELO
//
//   n features, d dimensões, n > d. Cada feature i tem um
//   vetor UNITÁRIO wᵢ ∈ R^d. O que a rede consegue ler de uma
//   ativação x é
//
//       x̂ = Wᵀ(Wx) = G x ,   G_ij = wᵢ·wⱼ
//
//   Como ‖wᵢ‖ = 1 para toda feature, G_ii = 1 sempre: a feature
//   responde a si mesma com força total. Fora da diagonal, G_ij é
//   exatamente o quanto a feature j vaza na leitura da feature i.
//   Se as wᵢ fossem ortonormais, G = I e x̂ = x — leitura perfeita.
//
//   Alocação (regra declarada, não aprendida):
//     · i < d   → direção dedicada e_i            (monossemântica)
//     · i ≥ d   → direção hash(i) mod d           (compartilhada)
//
//   Dentro de um bloco, as features se inclinam num arco de 0 a
//   π/2: wᵢ = cosθᵢ·eₚ + sinθᵢ·e_{p+1}. O ângulo θᵢ é a fração da
//   importância do bloco que vem ANTES de i. Consequência direta:
//   quanto menos importantes as features, mais elas se acumulam
//   no mesmo ângulo e mais interferem entre si — que é a
//   estratificação descrita em Scherlis et al. (2022).
// ─────────────────────────────────────────────────────────────

uint hash(uint x) {
    x ^= x >> 16; x *= 0x7feb352du;
    x ^= x >> 15; x *= 0x846ca68bu;
    x ^= x >> 16;
    return x;
}

int direcao(int i, int d) {
    if (i < d) return i;
    return int(hash(uint(i) + 0x9e3779b9u) % uint(d));
}

// Importância relativa da feature i. Só a razão entre importâncias
// importa; a escala some quando normalizamos o arco.
float importancia(int i, float decaimento) {
    return pow(1.0 - float(i) / float(N), decaimento);
}

// Ângulo de inclinação de i dentro do seu bloco, em [0, π/2).
// É a importância acumulada ANTES de i, sobre a do bloco inteiro.
// Com decaimento = 0 todas empatam e o arco fica uniforme; com
// decaimento alto as últimas quase não se separam.
float angulo(int i, int d, float decaimento) {
    // Em d = 1 não existe e_{p+1} distinta de e_p para onde inclinar:
    // só há uma direção e todas as features são o mesmo vetor. O arco
    // é vazio, não circular.
    if (d < 2) return 0.0;

    int p = direcao(i, d);
    float antes = 0.0;
    float bloco = 0.0;
    for (int j = 0; j < N; j++) {
        if (direcao(j, d) != p) continue;
        float imp = importancia(j, decaimento);
        bloco += imp;
        if (j < i) antes += imp;
    }
    if (bloco <= 0.0) return 0.0;
    return MEIA_VOLTA * (antes / bloco);
}

// G_ij = wᵢ·wⱼ para vetores unitários na base canônica de R^d.
//
//   wᵢ = cosθᵢ·e_{pᵢ} + sinθᵢ·e_{qᵢ},  qᵢ = (pᵢ+1) mod d
//
//   wᵢ·wⱼ = e_{pᵢ}·e_{pⱼ}cosθᵢcosθⱼ + e_{pᵢ}·e_{qⱼ}cosθᵢsinθⱼ
//         + e_{qᵢ}·e_{pⱼ}sinθᵢcosθⱼ + e_{qᵢ}·e_{qⱼ}sinθᵢsinθⱼ
//
// Cada termo só sobrevive quando os índices coincidem. Em d = 1 o
// ângulo é zero (ver angulo), e o resultado é G ≡ 1: com uma só
// dimensão todas as features são o mesmo vetor e não há como
// distinguir nenhuma delas — o caso mais extremo de superposição.
float gram(int i, int j, int d, float decaimento) {
    int pi = direcao(i, d);
    int pj = direcao(j, d);
    int qi = (pi + 1) % d;
    int qj = (pj + 1) % d;

    float ti = angulo(i, d, decaimento);
    float tj = angulo(j, d, decaimento);

    float g = 0.0;
    if (pi == pj) g += cos(ti) * cos(tj);
    if (pi == qj) g += cos(ti) * sin(tj);
    if (qi == pj) g += sin(ti) * cos(tj);
    if (qi == qj) g += sin(ti) * sin(tj);
    return g;
}

// ─────────────────────────────────────────────────────────────
// Rampa de cor: vazio → azul → âmbar → branco
// ─────────────────────────────────────────────────────────────
vec3 rampa(float v) {
    v = clamp(v, 0.0, 1.0);
    vec3 c0 = vec3(0.015, 0.020, 0.032);
    vec3 c1 = vec3(0.055, 0.120, 0.260);
    vec3 c2 = vec3(0.560, 0.380, 0.140);
    vec3 c3 = vec3(1.000, 0.960, 0.880);
    if (v < 0.34) return mix(c0, c1, v / 0.34);
    if (v < 0.72) return mix(c1, c2, (v - 0.34) / 0.38);
    return mix(c2, c3, (v - 0.72) / 0.28);
}

void main() {
    vec2 uv = gl_FragCoord.xy / uRes;

    int d = clamp(int(uParam1 + 0.5), 1, N);
    float decaimento = max(uParam2, 0.0);

    int col = clamp(int(uv.x * float(N)), 0, N - 1);
    int row = clamp(N - 1 - int(uv.y * float(N)), 0, N - 1);

    float g = gram(row, col, d, decaimento);
    float forca = clamp(g, 0.0, 1.0);
    bool naDiagonal = (row == col);

    // ── varredura: uma feature por vez é a que está ativa ──
    //    posição contínua, para o realce não pular de linha em linha
    float fase = mod(uTime * 0.62, float(N) + 5.0) - 2.5;
    float dist = float(row) - fase;
    float ativo = exp(-dist * dist * 0.30);

    // ── cor base: a diagonal é 1 por construção ──
    vec3 cor = rampa(forca) * 0.62;
    if (naDiagonal) {
        cor = mix(vec3(0.30, 0.62, 0.78), vec3(0.92, 0.98, 1.00), ativo * 0.7);
    }

    // A linha ativa mostra o conjunto de resposta daquela feature.
    float vazamento = naDiagonal ? 0.0 : forca;
    cor += vec3(0.95, 0.42, 0.16) * vazamento * ativo * 1.45;
    cor += vec3(0.85, 0.30, 0.10) * ativo * 0.05;

    // ── separação entre células ──
    vec2 cell = fract(uv * float(N));
    float borda = smoothstep(0.015, 0.075, min(min(cell.x, 1.0 - cell.x),
                                               min(cell.y, 1.0 - cell.y)));
    cor *= borda;

    // ── vinheta fraca ──
    vec2 c = uv - 0.5;
    cor *= 1.0 - dot(c, c) * 0.42;

    // ── grão mínimo: evita banding no degradê escuro ──
    //    (uint/float não existe em GLSL ES 3.0 — converte antes de dividir)
    float grao = float(hash(uint(gl_FragCoord.x) * 1973u
                          + uint(gl_FragCoord.y) * 9277u) % 255u) / 255.0;
    cor += grao * 0.008 - 0.004;

    fragColor = vec4(cor, 1.0);
}