// Stage E: port of Linas's in-group-mi-cluster (learn scm/gram-class/agglo-mi-rank.scm), as it runs today
// (STAGE_E_SPEC.md; similarities over Shape columns only, §5). Built in steps; this file currently does
// step E2a: the first round's seed pair, in-group (optimal-mi-in-group) and jaccard trim, printed in
// Linas's log format so the two logs can be diffed.
//
// Matrix model (shape-vec.scm, direct-sum.scm): rows are words or word-classes; a pair is either a Section
// (row = germ, column = ConnectorSeq) or a CrossSection (row = connector word, column = Shape). Pairs are
// explicit, with their own counts, so later steps can merge them the way gram-majority.scm does.
//
// cluster <counts.tsv> [--nrank 50] [--intended]
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define QUORUM 0.7
#define COMMONALITY 0.2
#define NOISE 4.0
#define MI_CUTOFF 4.0
// ---------------------------------------------------------------- interning
typedef struct { char **keys; size_t n, cap, mask; long *tab; } intern_t;
static uint64_t fnv (const char *s) {
uint64_t h = 1469598103934665603ULL;
while (*s) h = (h ^ (unsigned char) *s++) * 1099511628211ULL;
return h;
}
static void rehash (intern_t *t) {
t->mask = t->mask ? t->mask * 2 + 1 : 4095;
free(t->tab);
t->tab = malloc(sizeof(long) * (t->mask + 1));
for (size_t i = 0; i <= t->mask; i++) t->tab[i] = -1;
for (size_t k = 0; k < t->n; k++) {
size_t i = fnv(t->keys[k]) & t->mask;
while (t->tab[i] >= 0) i = (i + 1) & t->mask;
t->tab[i] = k;
}
}
static long lookup (intern_t *t, const char *s) {
if (!t->mask) return -1;
size_t i = fnv(s) & t->mask;
while (t->tab[i] >= 0) {
if (!strcmp(t->keys[t->tab[i]], s)) return t->tab[i];
i = (i + 1) & t->mask;
}
return -1;
}
static long intern (intern_t *t, const char *s) {
long k = lookup(t, s);
if (k >= 0) return k;
if (!t->mask) rehash(t);
if (t->n == t->cap) { t->cap = t->cap ? t->cap * 2 : 4096; t->keys = realloc(t->keys, sizeof(char *) * t->cap); }
t->keys[t->n] = strdup(s);
size_t i = fnv(s) & t->mask;
while (t->tab[i] >= 0) i = (i + 1) & t->mask;
t->tab[i] = t->n++;
if (t->n * 2 > t->mask) rehash(t);
return t->n - 1;
}
// ---------------------------------------------------------------- atoms and columns
// Rows: words and word-classes, keyed "W <name>" / "C <name>".
static intern_t atoms;
static const char *aname (long a) { return atoms.keys[a] + 2; }
// A column is a ConnectorSeq (shape = 0) or a Shape (shape = 1: point + connectors, one of them the hole).
typedef struct { int shape, nc, hole; long point; long *ca; char *cd; } col_t;
static intern_t colkeys;
static col_t *cols;
static size_t capcols;
static long make_col (int shape, long point, int nc, const long *ca, const char *cd, int hole) {
static char key[1 << 16];
int k = snprintf(key, sizeof key, "%c%ld", shape ? 'X' : 'S', shape ? point : -1L);
for (int i = 0; i < nc; i++) k += snprintf(key + k, sizeof key - k, " %ld%c", i == hole ? -1L : ca[i], cd[i]);
size_t before = colkeys.n;
long c = intern(&colkeys, key);
if ((size_t) c >= capcols) { capcols = capcols ? capcols * 2 : 1 << 16; cols = realloc(cols, sizeof(col_t) * capcols); }
if (colkeys.n > before) {
cols[c] = (col_t) { shape, nc, hole, point, malloc(sizeof(long) * nc), malloc(nc) };
memcpy(cols[c].ca, ca, sizeof(long) * nc);
memcpy(cols[c].cd, cd, nc);
}
return c;
}
// ---------------------------------------------------------------- pairs
typedef struct { long r, c; double n; int alive; } pair_t;
static pair_t *pairs;
static size_t npairs, cappairs, pmask;
static long *ptab;
static uint64_t phash (long r, long c) { return ((uint64_t) r * 0x9e3779b97f4a7c15ULL) ^ ((uint64_t) c * 0xbf58476d1ce4e5b9ULL); }
static void prehash (void) {
pmask = pmask ? pmask * 2 + 1 : (1 << 20) - 1;
free(ptab);
ptab = malloc(sizeof(long) * (pmask + 1));
for (size_t i = 0; i <= pmask; i++) ptab[i] = -1;
for (size_t k = 0; k < npairs; k++) {
size_t i = phash(pairs[k].r, pairs[k].c) & pmask;
while (ptab[i] >= 0) i = (i + 1) & pmask;
ptab[i] = k;
}
}
static long get_pair (long r, long c) {
if (!pmask) return -1;
size_t i = phash(r, c) & pmask;
while (ptab[i] >= 0) {
pair_t *p = &pairs[ptab[i]];
if (p->r == r && p->c == c) return p->alive ? ptab[i] : -1;
i = (i + 1) & pmask;
}
return -1;
}
static long make_pair (long r, long c) {
if (!pmask) prehash();
size_t i = phash(r, c) & pmask;
while (ptab[i] >= 0) {
pair_t *p = &pairs[ptab[i]];
if (p->r == r && p->c == c) { p->alive = 1; return ptab[i]; }
i = (i + 1) & pmask;
}
if (npairs == cappairs) { cappairs = cappairs ? cappairs * 2 : 1 << 18; pairs = realloc(pairs, sizeof(pair_t) * cappairs); }
pairs[npairs] = (pair_t) { r, c, 0, 1 };
ptab[i] = npairs++;
if (npairs * 2 > pmask) prehash();
return npairs - 1;
}
static double pair_count (long r, long c) { long p = get_pair(r, c); return p < 0 ? 0 : pairs[p].n; }
// ---------------------------------------------------------------- marginals (support.scm, transpose.scm)
typedef struct { double l0, l1, l2, lq, su; } norms_t;
static norms_t *rown, *coln;
static double *mmt, T_mmt, S_sum;
static long **rowp; // per row: its live pairs (right-stars), rebuilt by recompute()
static size_t *rowpn, rowp_n; // rowp_n: rows allocated (atoms can be added between recomputes)
static void norm_add (norms_t *m, double n) { m->l0 += 1; m->l1 += fabs(n); m->l2 += n * n; m->lq += sqrt(fabs(n)); m->su += n; }
static void recompute (void) {
free(rown); free(coln); free(mmt);
rown = calloc(atoms.n, sizeof(norms_t)); coln = calloc(colkeys.n, sizeof(norms_t)); mmt = calloc(atoms.n, sizeof(double));
if (rowp) { for (size_t r = 0; r < rowp_n; r++) free(rowp[r]); free(rowp); free(rowpn); }
rowp = calloc(atoms.n, sizeof(long *)); rowpn = calloc(atoms.n, sizeof(size_t)); rowp_n = atoms.n;
size_t *cap = calloc(atoms.n, sizeof(size_t));
for (size_t k = 0; k < npairs; k++) {
if (!pairs[k].alive) continue;
long r = pairs[k].r;
norm_add(&rown[r], pairs[k].n); norm_add(&coln[pairs[k].c], pairs[k].n);
if (rowpn[r] == cap[r]) { cap[r] = cap[r] ? cap[r] * 2 : 8; rowp[r] = realloc(rowp[r], sizeof(long) * cap[r]); }
rowp[r][rowpn[r]++] = k;
}
free(cap);
S_sum = 0; T_mmt = 0;
for (size_t r = 0; r < atoms.n; r++) S_sum += rown[r].su;
for (size_t k = 0; k < npairs; k++) if (pairs[k].alive) mmt[pairs[k].r] += coln[pairs[k].c].l1 * fabs(pairs[k].n);
for (size_t r = 0; r < atoms.n; r++) T_mmt += mmt[r];
}
// ---------------------------------------------------------------- similarities (symmetric-mi.scm, mi-similarity.scm)
static int as_run = 1;
typedef struct { long a, b; double fmi, rmi; int used; } sim_t;
static sim_t *sims;
static size_t nsims, capsims, smask;
static long *stab;
static uint64_t shash (long a, long b) { if (a > b) { long t = a; a = b; b = t; } return phash(a, b); }
static long find_sim (long a, long b) {
if (!smask) return -1;
size_t i = shash(a, b) & smask;
while (stab[i] >= 0) {
sim_t *s = &sims[stab[i]];
if ((s->a == a && s->b == b) || (s->a == b && s->b == a)) return stab[i];
i = (i + 1) & smask;
}
return -1;
}
static void srehash (void) {
smask = smask ? smask * 2 + 1 : 4095;
free(stab);
stab = malloc(sizeof(long) * (smask + 1));
for (size_t i = 0; i <= smask; i++) stab[i] = -1;
for (size_t k = 0; k < nsims; k++) {
size_t i = shash(sims[k].a, sims[k].b) & smask;
while (stab[i] >= 0) i = (i + 1) & smask;
stab[i] = k;
}
}
// make-gram-mi-simmer: compute and store (overwrite) the similarity of a and b.
static void simmer (long a, long b) {
double dot = 0;
for (size_t i = 0; i < rowpn[a]; i++) {
pair_t *p = &pairs[rowp[a][i]];
if (as_run && !cols[p->c].shape) continue; // §5: ConnectorSeq columns fail in-right-basis?
long q = get_pair(b, p->c);
if (q >= 0) dot += p->n * pairs[q].n;
}
double fmi = dot > 0 ? log2(dot * T_mmt / (mmt[a] * mmt[b])) : -INFINITY;
double rmi = fmi + 0.5 * log2((mmt[a] / T_mmt) * (mmt[b] / T_mmt)) + (-log2(T_mmt / (S_sum * S_sum)));
long k = find_sim(a, b);
if (k < 0) {
if (!smask) srehash();
if (nsims == capsims) { capsims = capsims ? capsims * 2 : 4096; sims = realloc(sims, sizeof(sim_t) * capsims); }
sims[nsims] = (sim_t) { a, b, fmi, rmi, 1 };
size_t i = shash(a, b) & smask;
while (stab[i] >= 0) i = (i + 1) & smask;
stab[i] = nsims++;
if (nsims * 2 > smask) srehash();
} else {
sims[k].fmi = fmi; sims[k].rmi = rmi;
}
if (fmi > 6) printf("\tMI(`%s`, `%s`) = %.6f rank-MI = %.6f\n", aname(a), aname(b), fmi, rmi);
}
static double mi_sim (long a, long b) { long k = find_sim(a, b); return k < 0 ? -INFINITY : sims[k].fmi; }
// ---------------------------------------------------------------- ranking (similarity.scm rank-words)
static long *ranked;
static size_t nranked;
// Diagnostic (CLUSTER_STALE_RANK): rank by the counts a row had when it was first ranked, as a test of
// whether Linas's rank-words reads pre-merge marginals.
static double *rank_l1, *rank_l0;
static size_t rank_n;
static int stale_rank;
static double rk1 (long a) { return stale_rank && (size_t) a < rank_n && rank_l1[a] >= 0 ? rank_l1[a] : rown[a].l1; }
static double rk0 (long a) { return stale_rank && (size_t) a < rank_n && rank_l1[a] >= 0 ? rank_l0[a] : rown[a].l0; }
static int by_rank (const void *x, const void *y) {
long a = *(const long *) x, b = *(const long *) y;
if (rk1(a) != rk1(b)) return rk1(a) > rk1(b) ? -1 : 1;
if (rk0(a) != rk0(b)) return rk0(a) > rk0(b) ? -1 : 1;
return (a > b) - (a < b);
}
static void rank_words (void) {
if (stale_rank) { // remember each row's counts the first time it is ranked
if (rank_n < atoms.n) {
rank_l1 = realloc(rank_l1, sizeof(double) * atoms.n); rank_l0 = realloc(rank_l0, sizeof(double) * atoms.n);
for (size_t r = rank_n; r < atoms.n; r++) rank_l1[r] = -1;
rank_n = atoms.n;
}
for (size_t r = 0; r < atoms.n; r++) if (rank_l1[r] < 0 && rowpn[r] > 0) { rank_l1[r] = rown[r].l1; rank_l0[r] = rown[r].l0; }
}
free(ranked);
ranked = malloc(sizeof(long) * atoms.n);
nranked = 0;
for (size_t r = 0; r < atoms.n; r++) if (rowpn[r] > 0) ranked[nranked++] = r;
qsort(ranked, nranked, sizeof(long), by_rank);
}
// loop-upper-diagonal 0..depth: compute sims not already stored.
static void sims_upper (size_t depth) {
if (depth > nranked) depth = nranked;
for (size_t i = 0; i < depth; i++) for (size_t j = i; j < depth; j++)
if (find_sim(ranked[i], ranked[j]) < 0) simmer(ranked[i], ranked[j]);
}
// ---------------------------------------------------------------- in-group (in-group.scm)
static int round_even (double x) { return (int) rint(x); } // Guile `round` = round half to even
// find-in-group: start with (WA WB), fold over candidates, cons accepted ones on the front.
static int find_in_group (long wa, long wb, double lower, double eps, const long *cand, int ncand, long *grp) {
double minscore = mi_sim(wa, wb) - eps;
int n = 2;
grp[0] = wa; grp[1] = wb; // list order: grp[0] is the head
for (int c = 0; c < ncand; c++) {
long w = cand[c];
if (w == wa || w == wb) continue;
int tight = n == 2 ? 2 : round_even(QUORUM * n);
int maxfail = n - tight, fail = 0, ok = 1;
for (int m = 0; m < n && ok; m++) {
double s = mi_sim(w, grp[m]);
if (s < lower) fail += maxfail + 999;
if (s < minscore) fail += 1;
ok = fail <= maxfail;
}
if (ok) { memmove(grp + 1, grp, sizeof(long) * n); grp[0] = w; n++; }
}
return n;
}
// optimal-mi-in-group (tightness 0.7, step 0.1, max 8.5, lower 1.0, max-jump 2.5, window 1.0, max 12).
static int optimal_in_group (long wa, long wb, const long *cand, int ncand, long *out) {
long grp[4096];
int nsteps = (int) rint(8.5 / 0.1), slots = (int) rint(1.0 / 0.1);
int window[64];
for (int i = 0; i < slots; i++) window[i] = 2;
double eps = 0;
for (int N = 1; N <= nsteps; N++) {
eps = N * 0.1;
int sz = find_in_group(wb, wa, 1.0, eps, cand, ncand, grp);
int prev = window[0];
memmove(window, window + 1, sizeof(int) * (slots - 1));
window[slots - 1] = sz;
if (!(sz - prev < 2.5 && sz < 12)) break;
}
int n = find_in_group(wa, wb, 1.0, eps - 0.1, cand, ncand, grp);
for (int i = 0; i < n; i++) out[i] = grp[n - 1 - i]; // reverse
return n;
}
// ---------------------------------------------------------------- jaccard trim (jaccard.scm, gram-majority.scm)
// count-shared-conseq: (shared, total) over columns with count > NOISE in some member.
static void shared_conseq (const long *g, int n, double *shared, double *total) {
int thresh = n <= 2 ? 2 : round_even(QUORUM * n), low = thresh - 1;
static long *seen; static size_t capseen;
static int *votes; static long *mark; static size_t capmark; static long epoch;
if (capmark < colkeys.n) { capmark = colkeys.n * 2; mark = realloc(mark, sizeof(long) * capmark); votes = realloc(votes, sizeof(int) * capmark); memset(mark, 0, sizeof(long) * capmark); }
epoch++;
size_t ns = 0;
for (int m = 0; m < n; m++) for (size_t i = 0; i < rowpn[g[m]]; i++) {
pair_t *p = &pairs[rowp[g[m]][i]];
if (!(p->n > NOISE)) continue;
if (mark[p->c] != epoch) {
mark[p->c] = epoch; votes[p->c] = 0;
if (ns == capseen) { capseen = capseen ? capseen * 2 : 1024; seen = realloc(seen, sizeof(long) * capseen); }
seen[ns++] = p->c;
}
}
int sh = 0;
for (size_t k = 0; k < ns; k++) {
int v = 0;
for (int m = 0; m < n; m++) if (pair_count(g[m], seen[k]) > NOISE) v++;
if (v > low) sh++;
}
*shared = sh; *total = ns;
}
static int trim_greedy (long *g, int n) {
double sh, tot;
shared_conseq(g, n, &sh, &tot);
double cm = sh / tot;
printf("Start: size=%d overlap = %.0f of %.0f disjuncts, commonality= %4.2f%%\n", n, sh, tot, cm * 100);
for (;;) {
if (n == 2 || COMMONALITY <= cm) break;
double best = cm, bsh = sh, btot = tot;
long bg[4096]; int bn = n, changed = 0;
for (int N = 0; N < n; N++) {
long m[4096]; int k = 0;
for (int i = 0; i < n; i++) if (i != N) m[k++] = g[i];
double s2, t2;
shared_conseq(m, k, &s2, &t2);
double c2 = s2 / t2;
if (best <= c2) { best = c2; bsh = s2; btot = t2; memcpy(bg, m, sizeof(long) * k); bn = k; changed = 1; }
if (COMMONALITY < c2) break;
}
if (!changed) break;
printf("Better: size=%d overlap = %.0f of %.0f disjuncts, commonality= %4.2f%%\n", bn, bsh, btot, best * 100);
memcpy(g, bg, sizeof(long) * bn); n = bn; cm = best; sh = bsh; tot = btot;
}
printf("Best: size=%d overlap = %.0f of %.0f disjuncts, commonality= %4.2f%%\n", n, sh, tot, cm * 100);
return n;
}
// ---------------------------------------------------------------- input
static int connectors (const char *seq, char words[][1024], char *dirs, int max) {
int n = 0;
const char *p = seq;
while ((p = strstr(p, "(Connector (WordNode \"")) && n < max) {
p += 22;
const char *e = strstr(p, "\") (SexNode \"");
if (!e || e - p >= 1024) return -1;
memcpy(words[n], p, e - p);
words[n][e - p] = 0;
dirs[n] = e[13];
p = e + 13;
n++;
}
return n;
}
static long word (const char *w) {
static char key[2048];
snprintf(key, sizeof key, "W %s", w);
return intern(&atoms, key);
}
static void load_sections (const char *path) {
FILE *f = fopen(path, "r");
if (!f) { perror(path); exit(1); }
static char line[1 << 16];
static char cw[64][1024];
char cd[64];
long ca[64];
const char *pre = "(Section (WordNode \"";
while (fgets(line, sizeof line, f)) {
if (strncmp(line, pre, strlen(pre))) continue;
char *g0 = line + strlen(pre), *g1 = strstr(g0, "\") (ConnectorSeq");
char *t1 = strchr(line, '\t'), *t2 = t1 ? strchr(t1 + 1, '\t') : NULL;
if (!g1 || !t1 || !t2) continue;
double a, b, n;
if (3 != sscanf(t2 + 1, "%lf %lf %lf", &a, &b, &n)) { fprintf(stderr, "bad count: %s", line); exit(1); }
*g1 = 0;
int nc = connectors(g1 + 3, cw, cd, 64);
if (nc <= 0) { fprintf(stderr, "bad section\n"); exit(1); }
long g = word(g0);
for (int i = 0; i < nc; i++) ca[i] = word(cw[i]);
long p = make_pair(g, make_col(0, -1, nc, ca, cd, -1)); // may realloc `pairs`
pairs[p].n = n;
}
fclose(f);
}
// explode-sections (shape-vec.scm:547-568): a CrossSection per connector, count copied from the Section.
static void explode_sections (void) {
size_t n0 = npairs;
for (size_t k = 0; k < n0; k++) {
if (!pairs[k].alive || cols[pairs[k].c].shape) continue;
col_t *s = &cols[pairs[k].c];
long ca[64]; char cd[64]; int nc = s->nc;
memcpy(ca, s->ca, sizeof(long) * nc); memcpy(cd, s->cd, nc);
long g = pairs[k].r; double n = pairs[k].n;
for (int i = 0; i < nc; i++) {
long x = make_pair(ca[i], make_col(1, g, nc, ca, cd, i));
pairs[x].n = n;
}
}
}
// ---------------------------------------------------------------- membership (gram-class-api.scm)
static int is_class (long a) { return atoms.keys[a][0] == 'C'; }
typedef struct { long w, cls; double n; } memb_t;
static memb_t *membs;
static size_t nmembs, capmembs;
static long find_memb (long w, long cls) {
for (size_t k = 0; k < nmembs; k++) if (membs[k].w == w && membs[k].cls == cls) return k;
return -1;
}
static long make_memb (long w, long cls) {
long k = find_memb(w, cls);
if (k >= 0) return k;
if (nmembs == capmembs) { capmembs = capmembs ? capmembs * 2 : 256; membs = realloc(membs, sizeof(memb_t) * capmembs); }
membs[nmembs] = (memb_t) { w, cls, 0 };
return nmembs++;
}
static int is_member (long a, long cls) { return a >= 0 && find_memb(a, cls) >= 0; }
// make-class-node (gram-majority.scm:86-109): members' names joined by spaces, ".i" until unique.
static long make_class_node (const long *g, int n) {
static char key[1 << 16];
int k = snprintf(key, sizeof key, "C ");
for (int i = 0; i < n; i++) k += snprintf(key + k, sizeof key - k, i ? " %s" : "%s", aname(g[i]));
while (lookup(&atoms, key) >= 0) k += snprintf(key + k, sizeof key - k, ".i");
return intern(&atoms, key);
}
// ---------------------------------------------------------------- sections and cross-sections (shape-vec.scm)
static int is_section (long p) { return !cols[pairs[p].c].shape; }
// make-cross-sections (shape-vec.scm:501-521): create (if needed) and return the crosses of a Section.
static int make_cross_sections (long sect, long *out) {
col_t s = cols[pairs[sect].c];
long ca[64]; char cd[64];
memcpy(ca, s.ca, sizeof(long) * s.nc); memcpy(cd, s.cd, s.nc);
long germ = pairs[sect].r;
for (int i = 0; i < s.nc; i++) {
long c = make_col(1, germ, s.nc, ca, cd, i);
out[i] = make_pair(ca[i], c);
}
return s.nc;
}
// get-cross-sections (shape-vec.scm:475-498): only those that exist.
static int get_cross_sections (long sect, long *out) {
col_t s = cols[pairs[sect].c];
static char key[1 << 16];
int n = 0;
for (int i = 0; i < s.nc; i++) {
int k = snprintf(key, sizeof key, "X%ld", pairs[sect].r);
for (int j = 0; j < s.nc; j++) k += snprintf(key + k, sizeof key - k, " %ld%c", j == i ? -1L : s.ca[j], s.cd[j]);
long c = lookup(&colkeys, key);
long x = c < 0 ? -1 : get_pair(s.ca[i], c);
if (x >= 0) out[n++] = x;
}
return n;
}
// The Section a CrossSection covers: germ goes back into the hole.
static long section_col_of (long cross) {
col_t sh = cols[pairs[cross].c];
long ca[64]; char cd[64];
memcpy(ca, sh.ca, sizeof(long) * sh.nc); memcpy(cd, sh.cd, sh.nc);
ca[sh.hole] = pairs[cross].r;
return make_col(0, -1, sh.nc, ca, cd, -1);
}
static long make_section (long cross) { long c = section_col_of(cross); return make_pair(cols[pairs[cross].c].point, c); }
static long get_section (long cross) { long c = section_col_of(cross); return get_pair(cols[pairs[cross].c].point, c); }
// flatten-section / flatten-cross / make-flat (shape-vec.scm:310-399).
static long make_flat (long cls, long pnt) {
col_t c = cols[pairs[pnt].c];
long ca[64]; char cd[64];
memcpy(ca, c.ca, sizeof(long) * c.nc); memcpy(cd, c.cd, c.nc);
long germ = pairs[pnt].r, point = c.point;
for (int i = 0; i < c.nc; i++) if (i != c.hole && is_member(ca[i], cls)) ca[i] = cls;
if (c.shape) { if (is_member(point, cls)) point = cls; ca[c.hole] = -1; }
if (is_member(germ, cls)) germ = cls;
long col = make_col(c.shape, c.shape ? point : -1, c.nc, ca, cd, c.shape ? c.hole : -1);
return make_pair(germ, col);
}
// is-nonflat? (shape-vec.scm:405-411): does a connector of the pair already name the class?
static int is_nonflat (long cls, long pnt) {
col_t c = cols[pairs[pnt].c];
for (int i = 0; i < c.nc; i++) if (i != c.hole && c.ca[i] == cls) return 1;
return 0;
}
// ---------------------------------------------------------------- counts (count-api.scm, shape-project.scm)
static void rebalance_count (long sect, double cnt) {
long xs[64];
pairs[sect].n = cnt;
int n = make_cross_sections(sect, xs);
for (int i = 0; i < n; i++) pairs[xs[i]].n = cnt;
}
static void rebalance_merge (long mrg, long donor) {
int is_sect = is_section(donor);
long m = is_sect ? mrg : make_section(mrg);
long d = is_sect ? donor : make_section(donor);
long xs[64];
make_cross_sections(m, xs);
rebalance_count(m, pairs[mrg].n);
rebalance_count(d, pairs[donor].n);
}
static double accumulate_count (long acc, long donor, double frac) {
double moved = frac * pairs[donor].n;
if (!(moved < 1.0e-10)) {
pairs[acc].n += moved;
pairs[donor].n -= moved;
rebalance_merge(acc, donor);
}
return moved;
}
// ---------------------------------------------------------------- merge-majority (gram-majority.scm:113-477)
static long cur_cls;
static long *cur_g;
static int cur_n, cur_thresh;
static int vote_to_accept (long col) {
int v = 0;
for (int i = 0; i < cur_n; i++) if (get_pair(cur_g[i], col) >= 0) v++;
return cur_thresh <= v;
}
static double do_merge (long w, long col, int accept) {
long sect = get_pair(w, col);
if (sect < 0) return 0;
int full = accept || pairs[sect].n <= NOISE || is_nonflat(cur_cls, sect) || is_class(w);
if (!full) return 0; // FRAC = 0
long flat = make_flat(cur_cls, sect);
double moved = accumulate_count(flat, sect, 1.0);
membs[make_memb(w, cur_cls)].n += moved;
return 1.0;
}
static void rebalance_dj (long col) {
for (int i = 0; i < cur_n; i++) {
long sect = get_pair(cur_g[i], col);
if (sect >= 0) { long flat = make_flat(cur_cls, sect); rebalance_merge(flat, sect); }
}
}
static char *done, *shape_done;
static size_t capdone;
static void mark (char **set, long col) {
if ((size_t) col >= capdone) {
size_t nc = colkeys.n * 2 + 1024;
done = realloc(done, nc); shape_done = realloc(shape_done, nc);
memset(done + capdone, 0, nc - capdone); memset(shape_done + capdone, 0, nc - capdone);
capdone = nc;
}
(*set)[col] = 1;
}
static int is_marked (char *set, long col) { return (size_t) col < capdone && set[col]; }
// make-once-predicate: false the first time it is asked about an atom, true afterwards.
static int once_shape_done (long col) { if (is_marked(shape_done, col)) return 1; mark(&shape_done, col); return 0; }
static void move_member_counts (long from) {
double old = 0, nw;
for (size_t k = 0; k < nmembs; k++) if (membs[k].cls == from) old += membs[k].n;
long d = find_memb(from, cur_cls);
nw = d < 0 ? 0 : membs[d].n;
double fract = nw / old;
for (size_t k = 0; k < nmembs; k++) {
if (membs[k].cls != from) continue;
double x = membs[k].n * fract;
long w = membs[k].w;
membs[make_memb(w, cur_cls)].n += x;
membs[find_memb(w, from)].n -= x;
}
if (d >= 0) { membs[d] = membs[--nmembs]; }
}
static void merge_majority (long cls, long *g, int n) {
for (int i = 0; i < n; i++) make_memb(g[i], cls);
cur_cls = cls; cur_g = g; cur_n = n;
cur_thresh = n <= 2 ? 2 : round_even(QUORUM * n);
free(done); free(shape_done); done = shape_done = NULL; capdone = 0;
mark(&done, 0); done[0] = 0;
// get-all-djs: the columns of all members' pairs, first-seen order (Linas: atom-set order).
size_t ndj = 0, capdj = 1024;
long *djs = malloc(sizeof(long) * capdj);
char *seen = calloc(colkeys.n + 1, 1);
for (int i = 0; i < n; i++) for (size_t k = 0; k < rowpn[g[i]]; k++) {
long c = pairs[rowp[g[i]][k]].c;
if (!pairs[rowp[g[i]][k]].alive || seen[c]) continue;
seen[c] = 1;
if (ndj == capdj) { capdj *= 2; djs = realloc(djs, sizeof(long) * capdj); }
djs[ndj++] = c;
}
free(seen);
int scnt = 0, mscnt = 0;
for (size_t k = 0; k < ndj; k++) {
long dj = djs[k];
if (cols[dj].shape) continue;
int have = vote_to_accept(dj);
double merged = 0;
for (int i = 0; i < n; i++) merged += do_merge(g[i], dj, have);
if (0 < merged) {
mark(&done, dj); // record-done: the DJ and the Shapes of its members' Sections
for (int i = 0; i < n; i++) {
long sect = get_pair(g[i], dj);
if (sect < 0) continue;
long xs[64];
int nx = make_cross_sections(sect, xs);
for (int j = 0; j < nx; j++) mark(&done, pairs[xs[j]].c);
}
mscnt++;
}
rebalance_dj(dj);
mark(&done, dj);
scnt++;
}
printf("------ merge-majority: Merge %d of %d sections\n", mscnt, scnt);
// Left-over CrossSections (gram-majority.scm:333-422).
int nleft = 0, mshcnt = 0;
long *left = malloc(sizeof(long) * (ndj + 1));
for (size_t k = 0; k < ndj; k++) if (!is_marked(done, djs[k])) left[nleft++] = djs[k];
if (getenv("CLUSTER_SHUFFLE")) { // diagnostic: does the left-over order (Linas: atom-set order) matter?
uint64_t z = strtoull(getenv("CLUSTER_SHUFFLE"), NULL, 10) * 0x9e3779b97f4a7c15ULL + 1;
for (int k = nleft - 1; k > 0; k--) {
z ^= z << 13; z ^= z >> 7; z ^= z << 17;
int j = (int) (z % (uint64_t) (k + 1));
long t = left[k]; left[k] = left[j]; left[j] = t;
}
}
for (int k = 0; k < nleft; k++) {
long shp = left[k];
if (once_shape_done(shp)) continue;
long alt = -1;
if (vote_to_accept(shp)) alt = shp;
else {
// get-alt-shapes: Shapes of the other crosses of the members' Sections, not yet done (marks them).
long cand[4096]; int nc = 0;
for (int i = 0; i < n; i++) {
long x = get_pair(g[i], shp);
if (x < 0) continue;
long sect = get_section(x);
if (sect < 0) continue;
long xs[64];
int nx = make_cross_sections(sect, xs);
for (int j = 0; j < nx; j++) {
long c = pairs[xs[j]].c;
int dup = 0;
for (int q = 0; q < nc; q++) if (cand[q] == c) dup = 1;
if (!dup && nc < 4096) cand[nc++] = c;
}
}
for (int q = 0; q < nc && alt < 0; q++) if (!once_shape_done(cand[q]) && vote_to_accept(cand[q])) alt = cand[q];
}
int have = alt >= 0;
long mrg = have ? alt : shp;
double mcnt = 0;
for (int i = 0; i < n; i++) mcnt += do_merge(g[i], mrg, have);
rebalance_dj(mrg);
if (0 < mcnt) mshcnt++;
}
printf("------ merge-majority: Remaining %d of %d cross\n", mshcnt, nleft);
for (int i = 0; i < n; i++) if (is_class(g[i])) move_member_counts(g[i]);
free(left); free(djs);
}
// ---------------------------------------------------------------- cleanup (shape-project.scm:310-387)
static void del_sect (long sect) {
long xs[64];
int nx = get_cross_sections(sect, xs);
for (int i = 0; i < nx; i++) pairs[xs[i]].alive = 0;
pairs[sect].alive = 0;
}
static void remove_empty (long row) {
size_t n0 = npairs;
for (size_t k = 0; k < n0; k++) {
if (!pairs[k].alive || pairs[k].r != row || !(pairs[k].n < 1.0e-10)) continue;
if (is_section(k)) del_sect(k);
else { long s = get_section(k); if (s >= 0) del_sect(s); }
}
}
// ---------------------------------------------------------------- main
int main (int argc, char **argv) {
if (argc < 2) { fprintf(stderr, "usage: cluster <counts.tsv> [--nrank 50] [--loops 5] [--intended] [--dump state.tsv]\n"); return 1; }
int nrank = 50, loops = 5;
const char *dump_path = NULL;
for (int i = 2; i < argc; i++) {
if (!strcmp(argv[i], "--nrank") && i + 1 < argc) nrank = atoi(argv[++i]);
else if (!strcmp(argv[i], "--loops") && i + 1 < argc) loops = atoi(argv[++i]);
else if (!strcmp(argv[i], "--intended")) as_run = 0;
else if (!strcmp(argv[i], "--dump") && i + 1 < argc) dump_path = argv[++i];
}
setvbuf(stdout, NULL, _IOLBF, 0);
stale_rank = getenv("CLUSTER_STALE_RANK") != NULL;
load_sections(argv[1]);
explode_sections();
recompute();
rank_words();
sims_upper(nrank);
printf("Done setting up similarity to %d\n", nrank);
for (int iter = 1; iter <= loops; iter++) {
// top-ranked-mi-pairs: fmi > 4, a != b, by ranked-MI; print the first 12.
printf("------ Round %d Next in line:\n", iter);
size_t nc = 0;
long *cand = malloc(sizeof(long) * (nsims + 1));
for (size_t k = 0; k < nsims; k++) if (sims[k].a != sims[k].b && MI_CUTOFF < sims[k].fmi) cand[nc++] = k;
for (size_t i = 1; i < nc; i++) { // insertion sort, rmi descending (stable)
long x = cand[i]; size_t j = i;
while (j > 0 && sims[cand[j - 1]].rmi < sims[x].rmi) { cand[j] = cand[j - 1]; j--; }
cand[j] = x;
}
for (size_t i = 0; i < nc && i < 12; i++)
printf("ranked-MI = %.4f MI = %.4f (`%s`, `%s`)\n", sims[cand[i]].rmi, sims[cand[i]].fmi, aname(sims[cand[i]].a), aname(sims[cand[i]].b));
if (!nc) { printf("no pair above the MI cutoff\n"); free(cand); break; }
long wa = sims[cand[0]].a, wb = sims[cand[0]].b;
free(cand);
printf("------ Start MI-based merge %d with seed pair `%s` and `%s`\n", iter, aname(wa), aname(wb));
rank_words();
int ncand = nranked < (size_t) (nrank + 3 * iter) ? (int) nranked : nrank + 3 * iter;
long grp[4096];
int n = optimal_in_group(wa, wb, ranked, ncand, grp);
printf("Initial in-group size=%d:", n);
for (int i = 0; i < n; i++) printf(" `%s`", aname(grp[i]));
printf("\n");
n = trim_greedy(grp, n);
printf("In-group size=%d:", n);
for (int i = 0; i < n; i++) printf(" `%s`", aname(grp[i]));
printf("\n");
long cls = make_class_node(grp, n);
merge_majority(cls, grp, n);
printf("------ Merged into `%s`\n", aname(cls));
if (getenv("CLUSTER_SHUFFLE")) { // diagnostic: class row and member counts after the merge
double cs = 0, cx = 0, ns = 0, nx = 0, mc = 0;
for (size_t k = 0; k < npairs; k++) if (pairs[k].alive && pairs[k].r == cls) { if (cols[pairs[k].c].shape) { cx += pairs[k].n; nx++; } else { cs += pairs[k].n; ns++; } }
for (size_t k = 0; k < nmembs; k++) if (membs[k].cls == cls) mc += membs[k].n;
printf("class: sections %.0f (count %.0f) crosses %.0f (count %.0f) member-count %.0f\n", ns, cs, nx, cx, mc);
}
recompute(); // marginals of everything touched (full recompute is equivalent)
remove_empty(cls);
for (int i = 0; i < n; i++) remove_empty(grp[i]);
recompute();
// recomp-all-sim for the in-group (PRECISE-SIM #f): every stored sim that involves a member.
for (int i = 0; i < n; i++) {
long w = grp[i];
int has = 0;
for (size_t k = 0; k < nsims && !has; k++) if (sims[k].a == w || sims[k].b == w) has = 1;
if (!has) continue;
simmer(w, w);
size_t n0 = nsims;
for (size_t k = 0; k < n0; k++) {
long o = sims[k].a == w ? sims[k].b : sims[k].b == w ? sims[k].a : -1;
if (o < 0 || o == w) continue;
int earlier = 0;
for (int j = 0; j < i; j++) if (grp[j] == o) earlier = 1;
if (!earlier) simmer(o, w);
}
}
simmer(cls, cls);
for (size_t r = 0; r < atoms.n; r++) if (is_class(r) && rowpn[r] > 0) simmer(cls, r);
// expand-universe: sims for missing pairs within the top NRANK + 2*iter.
rank_words();
sims_upper(nrank + 2 * iter);
if (getenv("CLUSTER_DEBUG")) { // diagnostic: where does the class rank, and its sims with members
for (size_t i = 0; i < nranked; i++) if (ranked[i] == cls) printf("debug: class rank %zu of %zu (count %.0f)\n", i, nranked, rown[cls].l1);
for (int i = 0; i < n; i++) { long k = find_sim(cls, grp[i]); printf("debug: sim(class, %s) = %s\n", aname(grp[i]), k < 0 ? "none" : ""); if (k >= 0) printf("debug: fmi %.4f rmi %.4f\n", sims[k].fmi, sims[k].rmi); }
{ long q = grp[0]; // the first member's remaining pairs: column, count, does the class have it?
for (size_t k = 0; k < rowpn[q] && k < 8; k++) {
pair_t *p = &pairs[rowp[q][k]];
printf("debug: %s pair col %ld [%s] n=%.0f alive=%d; class has it: %d\n", aname(q), p->c, colkeys.keys[p->c], p->n, p->alive, get_pair(cls, p->c) >= 0);
}
size_t nzero = 0; for (size_t k = 0; k < rowpn[q]; k++) if (pairs[rowp[q][k]].n < 1e-10) nzero++;
printf("debug: %s zero-count pairs %zu of %zu\n", aname(q), nzero, rowpn[q]);
}
for (int i = 0; i < n; i++) { // shared columns between the class and each member, by kind
int sh_x = 0, sh_s = 0; double dot_x = 0;
for (size_t k = 0; k < rowpn[cls]; k++) {
pair_t *p = &pairs[rowp[cls][k]];
long q = get_pair(grp[i], p->c);
if (q < 0) continue;
if (cols[p->c].shape) { sh_x++; dot_x += p->n * pairs[q].n; } else sh_s++;
}
printf("debug: class vs %s: shared shapes %d (dot %.0f), shared sections %d\n", aname(grp[i]), sh_x, dot_x, sh_s);
}
for (int i = 0; i < n; i++) {
size_t rk = nranked; for (size_t j = 0; j < nranked; j++) if (ranked[j] == grp[i]) rk = j;
double secs = 0, xs = 0;
for (size_t k = 0; k < rowpn[grp[i]]; k++) { pair_t *p = &pairs[rowp[grp[i]][k]]; if (cols[p->c].shape) xs += p->n; else secs += p->n; }
printf("debug: member %s rank %zu count %.0f (sections %.0f, crosses %.0f) pairs %zu\n", aname(grp[i]), rk, rown[grp[i]].l1, secs, xs, rowpn[grp[i]]);
}
}
}
// Final state (level2/PREREG-LV.md): atoms (id, is_class, name), memberships, live Section pairs.
if (dump_path) {
FILE *f = fopen(dump_path, "w");
if (!f) { perror(dump_path); return 1; }
for (size_t a = 0; a < atoms.n; a++) fprintf(f, "A\t%zu\t%d\t%s\n", a, is_class(a), aname(a));
for (size_t k = 0; k < nmembs; k++) fprintf(f, "M\t%ld\t%ld\t%.6f\n", membs[k].w, membs[k].cls, membs[k].n);
for (size_t k = 0; k < npairs; k++) {
if (!pairs[k].alive || pairs[k].n <= 0 || cols[pairs[k].c].shape) continue;
col_t c = cols[pairs[k].c];
fprintf(f, "S\t%ld\t%.6f\t", pairs[k].r, pairs[k].n);
for (int i = 0; i < c.nc; i++) fprintf(f, i ? " %ld%c" : "%ld%c", c.ca[i], c.cd[i]);
fprintf(f, "\n");
}
fclose(f);
}
return 0;
}