/* engram_cognition.c — THE ONE OPERATION. See engram_cognition.h.
 * Pure over its inputs (think/warp/express); persistence is additive/supersede
 * only. stdlib + libm + engram_store/reason/geometry. Touches no live daemon. */
#include "engram_cognition.h"
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>

/* ── small helpers ──────────────────────────────────────────────────────────── */
static double vdot(const float* a, const float* b, int dim) {
    double s = 0; for (int i = 0; i < dim; i++) s += (double)a[i] * (double)b[i]; return s;
}
static double vnorm(const float* a, int dim) { return sqrt(vdot(a, a, dim)); }
static char* dupstr(const char* s) {
    if (!s) return NULL; size_t n = strlen(s) + 1; char* p = malloc(n);
    if (p) memcpy(p, s, n); return p;
}
static double clampd(double x, double lo, double hi){ return x<lo?lo:(x>hi?hi:x); }

/* ═══════════════════════════════════════════════ Stance lifecycle ════════════ */
int cog_stance_init(CogStance* s, const char* id, const char* faculty,
                    const char* anchor_region, const char* for_whom,
                    const GeoDescriptor* region) {
    if (!s || !region) return -1;
    memset(s, 0, sizeof *s);
    s->id = dupstr(id); s->faculty = dupstr(faculty);
    s->anchor_region = dupstr(anchor_region); s->for_whom = dupstr(for_whom);
    s->dim = region->dim;
    s->n_axes = region->n_axes > COG_MAX_AXES ? COG_MAX_AXES : region->n_axes;
    for (int k = 0; k < COG_MAX_AXES; k++) s->axis_gain[k] = 1.0;
    s->ext_floor = 1.0; s->drop_frac = 0.5; s->assoc_floor = 0.2;
    s->bias_dir = NULL;
    s->reliability = 0.5;   /* uninformed prior on our own track record */
    return 0;
}
void cog_stance_set_frozen_defaults(CogStance* s) {
    if (!s) return;
    for (int k = 0; k < COG_MAX_AXES; k++) s->axis_gain[k] = 1.0;
    s->ext_floor = 1.0; s->drop_frac = 0.5; s->assoc_floor = 0.2;
    free(s->bias_dir); s->bias_dir = NULL;
}
void cog_stance_free(CogStance* s) {
    if (!s) return;
    free(s->id); free(s->faculty); free(s->anchor_region); free(s->for_whom);
    free(s->bias_dir);
    s->id = s->faculty = s->anchor_region = s->for_whom = NULL; s->bias_dir = NULL;
}
int cog_is_keystone(const CogKeystoneSet* ks, const CogStance* s) {
    if (!s) return 0;
    if (s->keystone) return 1;
    if (!ks || !s->id) return 0;
    for (int i = 0; i < ks->n; i++)
        if (ks->ids[i] && (strcmp(ks->ids[i], s->id) == 0 ||
            (s->anchor_region && strcmp(ks->ids[i], s->anchor_region) == 0))) return 1;
    return 0;
}

/* ═══════════════════════════════════════════════ warped fit (think step 2) ══ */
int cog_warped_fit(const GeoDescriptor* g, const float* x,
                   const CogStance* st, GeoFit* out) {
    if (!g || !x || !out || g->dim <= 0 || !g->centroid) return -1;
    double ext_floor = (st && st->ext_floor > 0) ? st->ext_floor : 1.0;
    int dim = g->dim;
    double rr = 0;
    float* r = malloc((size_t)dim * sizeof(float));
    if (!r) return -1;
    for (int i = 0; i < dim; i++) { double d = (double)x[i] - (double)g->centroid[i]; r[i] = (float)d; rr += d * d; }
    double maha2 = 0, ss_in = 0;
    for (int k = 0; k < g->n_axes; k++) {
        const float* ax = g->axes[k].axis; if (!ax) continue;
        double proj = vdot(r, ax, dim);
        double gain = (st && k < st->n_axes && st->axis_gain[k] > 0) ? st->axis_gain[k] : 1.0;
        double den = g->axes[k].extent * gain; if (den < ext_floor) den = ext_floor;
        maha2 += (proj / den) * (proj / den);
        ss_in += proj * proj;
    }
    double ortho2 = rr - ss_in; if (ortho2 < 0) ortho2 = 0;
    double dist2 = maha2 + ortho2 / (ext_floor * ext_floor);
    out->mahalanobis = sqrt(maha2); out->ortho_residual = sqrt(ortho2);
    out->distance = sqrt(dist2);    out->score = 1.0 / (1.0 + dist2);
    free(r);
    return 0;
}

/* ═══════════════════════════════════════════════ think (the ONE operation) ══ */
void engram_gradient_free(GeoGradient* g) {
    if (!g) return; free(g->direction); g->direction = NULL;
}
int engram_think(const GeoDescriptor* region, const float* anchor,
                 const CogStance* stance, GeoGradient* out) {
    if (!region || !out || region->dim <= 0 || !region->centroid) return -1;
    int dim = region->dim;
    memset(out, 0, sizeof *out);
    out->dim = dim;
    const float* x = anchor ? anchor : region->centroid;   /* re-origin (step 1) */

    GeoFit f;
    if (cog_warped_fit(region, x, stance, &f) != 0) return -1;   /* fit (step 2) */

    /* step 3 — emit a GRADIENT: warped steepest DESCENT of the fit distance². */
    double ext_floor = (stance && stance->ext_floor > 0) ? stance->ext_floor : 1.0;
    float* grad = calloc((size_t)dim, sizeof(float));   /* ∇ dist² wrt x         */
    float* r    = malloc((size_t)dim * sizeof(float));
    out->direction = malloc((size_t)dim * sizeof(float));
    if (!grad || !r || !out->direction) { free(grad); free(r); free(out->direction); out->direction = NULL; return -1; }
    for (int i = 0; i < dim; i++) r[i] = (float)((double)x[i] - (double)region->centroid[i]);
    /* in-subspace: Σ_k 2 (proj/den²) a_k ; also accumulate Σ proj a_k for ortho part */
    float* proj_sum = calloc((size_t)dim, sizeof(float));
    if (!proj_sum) { free(grad); free(r); free(out->direction); out->direction = NULL; return -1; }
    for (int k = 0; k < region->n_axes; k++) {
        const float* ax = region->axes[k].axis; if (!ax) continue;
        double proj = vdot(r, ax, dim);
        double gain = (stance && k < stance->n_axes && stance->axis_gain[k] > 0) ? stance->axis_gain[k] : 1.0;
        double den = region->axes[k].extent * gain; if (den < ext_floor) den = ext_floor;
        double coef = 2.0 * proj / (den * den);
        for (int i = 0; i < dim; i++) { grad[i] += (float)(coef * ax[i]); proj_sum[i] += (float)(proj * ax[i]); }
    }
    /* orthogonal: (2 r − 2 Σ proj a_k) / ext_floor² */
    double inv_f2 = 1.0 / (ext_floor * ext_floor);
    for (int i = 0; i < dim; i++)
        grad[i] += (float)((2.0 * (double)r[i] - 2.0 * (double)proj_sum[i]) * inv_f2);
    free(proj_sum);

    /* steering = −grad (descent), seeded by the stance's bias_dir. */
    for (int i = 0; i < dim; i++) out->direction[i] = -grad[i];
    if (stance && stance->bias_dir) {
        double gn = vnorm(grad, dim), bn = vnorm(stance->bias_dir, dim);
        if (bn > 1e-12) {
            double scale = (gn > 1e-12 ? gn : 1.0);   /* seed at the gradient's scale */
            for (int i = 0; i < dim; i++)
                out->direction[i] += (float)(scale * (double)stance->bias_dir[i] / bn);
        }
    }
    double dn = vnorm(out->direction, dim);
    if (dn > 1e-12) for (int i = 0; i < dim; i++) out->direction[i] /= (float)dn;
    else            for (int i = 0; i < dim; i++) out->direction[i] = 0.0f;   /* at rest */

    out->spread     = f.distance;                       /* spiked (0) .. diffuse   */
    out->confidence = stance ? stance->reliability : 0.5;
    out->magnitude  = f.score;                          /* the read's membership   */
    out->anchor_id  = region->hub_id;                   /* borrowed vantage id     */
    out->n_support  = region->n_members;
    out->stance_id  = stance ? stance->id : NULL;
    free(grad); free(r);
    return 0;
}

/* EXPRESSION — the ONLY collapse to a point (a separate faculty from think). */
int engram_express(const GeoGradient* g, const float* anchor, float* out_point) {
    if (!g || !anchor || !out_point || !g->direction) return -1;
    double commit = clampd(g->confidence, 0.0, 1.0);    /* confident => commit far */
    for (int i = 0; i < g->dim; i++)
        out_point[i] = anchor[i] + g->direction[i] * (float)commit;
    return 0;
}

/* ═══════════════════════════════════════════════ Stance serialization ════════ */
/* Compact line schema "STNC1" (mirrors the reify "GEO1" precedent). */
char* cog_stance_to_metadata(const CogStance* s) {
    if (!s) return NULL;
    size_t cap = 256 + (size_t)s->n_axes * 24 + (size_t)(s->bias_dir ? s->dim * 16 : 0);
    char* buf = malloc(cap); if (!buf) return NULL;
    size_t o = 0;
    o += (size_t)snprintf(buf + o, cap - o, "%s\n", COG_STANCE_META_MAGIC);
    o += (size_t)snprintf(buf + o, cap - o, "f %s\n", s->faculty ? s->faculty : "-");
    o += (size_t)snprintf(buf + o, cap - o, "r %s\n", s->anchor_region ? s->anchor_region : "-");
    o += (size_t)snprintf(buf + o, cap - o, "w %s\n", s->for_whom ? s->for_whom : "-");
    o += (size_t)snprintf(buf + o, cap - o, "k %d\n", s->keystone);
    o += (size_t)snprintf(buf + o, cap - o, "d %d %d\n", s->dim, s->n_axes);
    o += (size_t)snprintf(buf + o, cap - o, "s %.9g %.9g %.9g\n", s->ext_floor, s->drop_frac, s->assoc_floor);
    o += (size_t)snprintf(buf + o, cap - o, "g");
    for (int k = 0; k < s->n_axes; k++) o += (size_t)snprintf(buf + o, cap - o, " %.9g", s->axis_gain[k]);
    o += (size_t)snprintf(buf + o, cap - o, "\n");
    o += (size_t)snprintf(buf + o, cap - o, "c %lld %.9g %.9g %.9g %.9g\n",
                          (long long)s->n_trials, s->brier_sum, s->reliability, s->ema_error, s->last_error);
    if (s->bias_dir) {
        o += (size_t)snprintf(buf + o, cap - o, "b");
        for (int i = 0; i < s->dim; i++) o += (size_t)snprintf(buf + o, cap - o, " %.9g", (double)s->bias_dir[i]);
        o += (size_t)snprintf(buf + o, cap - o, "\n");
    }
    (void)o;
    return buf;
}
int cog_stance_to_node(const CogStance* s, StoreNode* out) {
    if (!s || !out) return -1;
    memset(out, 0, sizeof *out);
    out->id        = dupstr(s->id);
    out->node_type = dupstr(COG_STANCE_NODE_TYPE);
    out->content   = dupstr(s->faculty ? s->faculty : "stance");
    out->label     = dupstr(s->faculty ? s->faculty : "stance");
    out->metadata  = cog_stance_to_metadata(s);
    out->importance = s->reliability;    /* cached denormalized readout (§2.1)     */
    out->confidence = s->reliability;
    out->temporal_decay_rate = 0.0;
    return (out->id && out->node_type && out->metadata) ? 0 : -1;
}
static int parse_floats(const char* line, double* out, int max) {
    int n = 0; const char* p = line;
    while (*p && n < max) {
        while (*p == ' ') p++;
        if (!*p) break;
        char* end; double v = strtod(p, &end);
        if (end == p) break;
        out[n++] = v; p = end;
    }
    return n;
}
int cog_stance_from_node(const StoreNode* n, CogStance* out) {
    if (!n || !out || !n->metadata) return -1;
    memset(out, 0, sizeof *out);
    for (int k = 0; k < COG_MAX_AXES; k++) out->axis_gain[k] = 1.0;
    out->ext_floor = 1.0; out->drop_frac = 0.5; out->assoc_floor = 0.2; out->reliability = 0.5;
    out->id = dupstr(n->id);
    /* verify magic on first line */
    const char* m = n->metadata;
    if (strncmp(m, COG_STANCE_META_MAGIC, strlen(COG_STANCE_META_MAGIC)) != 0) return -1;
    char* copy = dupstr(m); if (!copy) return -1;
    for (char* line = strtok(copy, "\n"); line; line = strtok(NULL, "\n")) {
        if (line[0] == '\0' || line[1] != ' ') {
            if (line[0] == 'g' || line[0] == 'b') { /* vector lines: tag then values */ }
            else continue;
        }
        char tag = line[0];
        const char* rest = line + 1; while (*rest == ' ') rest++;
        if (tag == 'f')      { free(out->faculty); out->faculty = (strcmp(rest, "-") ? dupstr(rest) : NULL); }
        else if (tag == 'r') { free(out->anchor_region); out->anchor_region = (strcmp(rest, "-") ? dupstr(rest) : NULL); }
        else if (tag == 'w') { free(out->for_whom); out->for_whom = (strcmp(rest, "-") ? dupstr(rest) : NULL); }
        else if (tag == 'k') { out->keystone = atoi(rest); }
        else if (tag == 'd') { int a=0,b=0; sscanf(rest, "%d %d", &a, &b); out->dim = a; out->n_axes = b > COG_MAX_AXES ? COG_MAX_AXES : b; }
        else if (tag == 's') { double v[3]={1,0.5,0.2}; parse_floats(rest, v, 3); out->ext_floor=v[0]; out->drop_frac=v[1]; out->assoc_floor=v[2]; }
        else if (tag == 'g') { double v[COG_MAX_AXES]; int c=parse_floats(rest, v, COG_MAX_AXES); for(int k=0;k<c;k++) out->axis_gain[k]=v[k]; }
        else if (tag == 'c') { double v[5]={0,0,0.5,0,0}; parse_floats(rest, v, 5); out->n_trials=(int64_t)v[0]; out->brier_sum=v[1]; out->reliability=v[2]; out->ema_error=v[3]; out->last_error=v[4]; }
        else if (tag == 'b') { if (out->dim>0){ out->bias_dir=calloc((size_t)out->dim,sizeof(float)); double v[4096]; int c=parse_floats(rest,v,out->dim<4096?out->dim:4096); for(int i=0;i<c;i++) out->bias_dir[i]=(float)v[i]; } }
    }
    free(copy);
    return 0;
}

/* ═══════════════════════════════════════════════ grounding as a RELATION ═════ */
static int put_edge(EngramPagedStore* s, const char* id, const char* from, const char* to,
                    const char* relation, double weight, const char* meta) {
    StoreEdge e; memset(&e, 0, sizeof e);
    e.id = (char*)id; e.from_id = (char*)from; e.to_id = (char*)to;
    e.relation = (char*)relation; e.weight = weight; e.confidence = weight;
    e.metadata = (char*)meta;
    return store_put_edge(s, &e);
}
int cog_ground_edge(EngramPagedStore* s, const char* claim_id,
                    const char* evidence_id, double grounding, const char* for_whom) {
    if (!s || !claim_id || !evidence_id) return -1;
    char id[512], meta[256];
    snprintf(id, sizeof id, "gb-%s-%s-%s", claim_id, evidence_id, for_whom ? for_whom : "global");
    snprintf(meta, sizeof meta, "for_whom=%s", for_whom ? for_whom : "-");
    return put_edge(s, id, claim_id, evidence_id, COG_GROUNDED_BY_RELATION, grounding, meta);
}
int cog_salient_edge(EngramPagedStore* s, const char* node_id,
                     const char* observer_id, double salience) {
    if (!s || !node_id || !observer_id) return -1;
    char id[512];
    snprintf(id, sizeof id, "st-%s-%s", node_id, observer_id);
    return put_edge(s, id, node_id, observer_id, COG_SALIENT_TO_RELATION, salience, NULL);
}
int cog_assert_gate(EngramPagedStore* s, const char* claim_id,
                    const char* for_whom, double floor) {
    if (!s || !claim_id) return -1;
    if (!(floor > 0)) floor = 0.5;
    StoreEdge* edges = NULL; size_t n = 0;
    if (store_get_edges_from(s, claim_id, &edges, &n) < 0) return -1;
    double best = 0.0; int found = 0;
    for (size_t i = 0; i < n; i++) {
        if (!edges[i].relation || strcmp(edges[i].relation, COG_GROUNDED_BY_RELATION) != 0) continue;
        /* grounded-for-whom: match observer if requested; global (for_whom=-) always counts */
        int match = 1;
        if (for_whom && edges[i].metadata) {
            const char* fw = strstr(edges[i].metadata, "for_whom=");
            if (fw) { fw += 9; if (strcmp(fw, for_whom) != 0 && strcmp(fw, "-") != 0) match = 0; }
        }
        if (match) { found = 1; if (edges[i].weight > best) best = edges[i].weight; }
    }
    store_edges_free(edges, n);
    if (!found) return 0;                 /* ungrounded => refuse assertion (still held) */
    return (best >= floor) ? 1 : 0;
}

/* ═══════════════════════════════════════════════ THE CORRESPONDENCE-LOOP ═════ */
int engram_correspondence_beat(const GeoDescriptor* region, const float* anchor,
                               double outcome_y, CogStance* stance,
                               int learn, double max_step, CogBeatResult* out) {
    if (!region || !stance || !out) return -1;
    memset(out, 0, sizeof *out);
    if (stance->keystone) { learn = 0; out->wrote_keystone = 1; }  /* §6: never write a keystone */

    GeoGradient g;
    if (engram_think(region, anchor, stance, &g) != 0) return -1;   /* PREDICTION */
    double p = g.magnitude;
    double y = clampd(outcome_y, 0.0, 1.0);
    double err = fabs(p - y);
    out->correspondence = 1.0 - err;
    out->error = err;
    out->brier = (p - y) * (p - y);

    if (learn) {
        /* refine warp: gradient descent of (p−y)² wrt each axis_gain.
         *   p = 1/(1+D²);  ∂p/∂gain_k = 2 p² proj_k² / (ext_k² gain_k³)  (>=0)
         *   ∂(err²)/∂gain_k = 2 (p−y) ∂p/∂gain_k
         * step = −lr · ∂(err²)/∂gain_k, bounded to ±max_step (metastability). */
        int dim = region->dim;
        const float* x = anchor ? anchor : region->centroid;
        float* r = malloc((size_t)dim * sizeof(float));
        if (r) {
            for (int i = 0; i < dim; i++) r[i] = (float)((double)x[i] - (double)region->centroid[i]);
            double lr = 0.5;
            double bound = (max_step > 0) ? max_step : 0.05;   /* bounded update rate */
            for (int k = 0; k < region->n_axes && k < stance->n_axes; k++) {
                const float* ax = region->axes[k].axis; if (!ax) continue;
                double proj = vdot(r, ax, dim);
                double ext = region->axes[k].extent; if (ext < 1e-9) ext = 1e-9;
                double gain = stance->axis_gain[k]; if (gain < 1e-6) gain = 1e-6;
                double dp_dgain = 2.0 * p * p * (proj * proj) / (ext * ext * gain * gain * gain);
                double dErr_dgain = 2.0 * (p - y) * dp_dgain;
                double step = -lr * dErr_dgain;
                step = clampd(step, -bound, bound);
                stance->axis_gain[k] = clampd(gain + step, 0.1, 50.0);
            }
            free(r);
        }
        /* calibration */
        stance->n_trials += 1;
        stance->brier_sum += out->brier;
        stance->last_error = err;
        stance->ema_error = (stance->n_trials == 1) ? err : 0.9 * stance->ema_error + 0.1 * err;
        double mean_brier = stance->brier_sum / (double)stance->n_trials;
        stance->reliability = clampd(1.0 - sqrt(mean_brier), 0.0, 1.0);
    }
    out->reliability = stance->reliability;
    engram_gradient_free(&g);
    return 0;
}