// Stage E, step 1: Linas's covering matrix and the initial grammatical-MI similarities (see STAGE_E_SPEC.md §0-1).
//
// * Sections from the counts file (output of `count`), exploded into CrossSections as in
// learn shape-vec.scm explode-sections: Section (g, c1..cn) with count N gives, for each i,
// CrossSection [ci.word, <g, c1 .. ($, ci.dir) .. cn>] with count N (set, not added: a duplicate
// (row, column) is an error here).
// * Rows = words (germs of Sections and of CrossSections); columns = ConnectorSeqs and Shapes.
// * Marginals (matrix support.scm) and MM^T (transpose.scm) on that matrix, as in csetmarg.c.
// * rank-words (similarity.scm:17-42): right-count desc, ties right-support desc.
// * Similarities for all pairs i <= j among the top NRANK ranked words (agglo-mi-rank.scm:535,
// loop-upper-diagonal): fmi = log2(sum_d N(a,d)N(b,d) * T / (m(a) m(b))) and
// rmi = fmi + 0.5 log2(m(a)/T * m(b)/T) + Q, Q = -log2(T / S^2) (symmetric-mi.scm:170-264,
// mi-similarity.scm:49-53); m = MM^T count of the row, T = MM^T total count, S = total-sum-left.
// --as-run: reproduce Linas's stack as it runs today, where the dot product sees Shape columns only. The
// covering object's in-right-basis? is built from a Query on (make-pair $x $y), and direct-sum's make-pair
// with two variables always builds a CrossSection (direct-sum.scm:234-243, make-left-variable-pair tests
// the right atom's type, which is a variable), so add-fast-math (fold-api.scm:351-353) drops every
// ConnectorSeq column. Marginals do not go through that filter and still include Sections.
// Output: one line per similarity: <word a>\t<word b>\t<fmi>\t<rmi>. With --top K, the K best pairs with
// fmi > 4 (the seed list of agglo-mi-rank.scm:628-639) are printed to stderr in Linas's log format.
//
// cover <counts.tsv> [--as-run] [--nrank 50] [--top 12] > sims.tsv
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
// ---------------------------------------------------------------- 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 intern (intern_t *t, const char *s) {
if (!t->mask) rehash(t);
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;
}
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);
t->tab[i] = t->n++;
if (t->n * 2 > t->mask) rehash(t);
return t->n - 1;
}
// ---------------------------------------------------------------- matrix
typedef struct { long r, c; double n; } entry_t;
typedef struct { double l0, l1, l2, lq, su; } norms_t;
static intern_t rows, cols;
static entry_t *ents;
static size_t nents, capents;
static void add_entry (long r, long c, double n) {
if (nents == capents) { capents = capents ? capents * 2 : 1 << 18; ents = realloc(ents, sizeof(entry_t) * capents); }
ents[nents++] = (entry_t) { r, c, n };
}
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 int by_row_col (const void *a, const void *b) {
const entry_t *x = a, *y = b;
if (x->r != y->r) return x->r < y->r ? -1 : 1;
return (x->c > y->c) - (x->c < y->c);
}
// Parse the "(Connector (WordNode \"w\") (SexNode \"d\"))" items of a ConnectorSeq string.
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 double *rank_count, *rank_supp;
static int by_rank (const void *a, const void *b) {
long x = *(const long *) a, y = *(const long *) b;
if (rank_count[x] != rank_count[y]) return rank_count[x] > rank_count[y] ? -1 : 1;
if (rank_supp[x] != rank_supp[y]) return rank_supp[x] > rank_supp[y] ? -1 : 1;
return (x > y) - (x < y);
}
typedef struct { long a, b; double fmi, rmi; } sim_t;
static int by_rmi (const void *x, const void *y) {
const sim_t *p = x, *q = y;
return (p->rmi < q->rmi) - (p->rmi > q->rmi);
}
int main (int argc, char **argv) {
if (argc < 2) { fprintf(stderr, "usage: cover <counts.tsv> [--as-run] [--nrank 50] [--top 12]\n"); return 1; }
int nrank = 50, top = 12, as_run = 0;
for (int i = 2; i < argc; i++) {
if (!strcmp(argv[i], "--as-run")) as_run = 1;
else if (!strcmp(argv[i], "--nrank") && i + 1 < argc) nrank = atoi(argv[++i]);
else if (!strcmp(argv[i], "--top") && i + 1 < argc) top = atoi(argv[++i]);
}
FILE *f = fopen(argv[1], "r");
if (!f) { perror(argv[1]); return 1; }
static char line[1 << 16], seq[1 << 16], shape[1 << 16];
static char cw[64][1024];
char cd[64];
const char *pre = "(Section (WordNode \"";
long nsect = 0;
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); return 1; }
*g1 = 0;
t1[-1] = 0;
strcpy(seq, g1 + 3);
int nc = connectors(seq, cw, cd, 64);
if (nc <= 0) { fprintf(stderr, "bad section: %s\n", seq); return 1; }
add_entry(intern(&rows, g0), intern(&cols, seq), n);
nsect++;
for (int i = 0; i < nc; i++) { // explode-sections: hole at connector i
int k = snprintf(shape, sizeof shape, "(ShapeLink (WordNode \"%s\")", g0);
for (int j = 0; j < nc; j++) {
if (j == i) k += snprintf(shape + k, sizeof shape - k, " (Connector (VariableNode \"$connector-word\") (SexNode \"%c\"))", cd[j]);
else k += snprintf(shape + k, sizeof shape - k, " (Connector (WordNode \"%s\") (SexNode \"%c\"))", cw[j], cd[j]);
}
snprintf(shape + k, sizeof shape - k, ")");
add_entry(intern(&rows, cw[i]), intern(&cols, shape), n);
}
}
fclose(f);
qsort(ents, nents, sizeof(entry_t), by_row_col);
for (size_t k = 1; k < nents; k++)
if (ents[k].r == ents[k - 1].r && ents[k].c == ents[k - 1].c) {
fprintf(stderr, "duplicate (row, column): %s %s\n", rows.keys[ents[k].r], cols.keys[ents[k].c]);
return 2;
}
fprintf(stderr, "sections %ld, crosses %zu, rows %zu, columns %zu\n", nsect, nents - nsect, rows.n, cols.n);
// support.scm marginals: right = per row, left = per column.
norms_t *rn = calloc(rows.n, sizeof(norms_t)), *cn = calloc(cols.n, sizeof(norms_t));
for (size_t k = 0; k < nents; k++) { norm_add(&rn[ents[k].r], ents[k].n); norm_add(&cn[ents[k].c], ents[k].n); }
double S = 0; // total-sum-left = sum over rows of right-sum (support.scm do-left-totals)
for (size_t r = 0; r < rows.n; r++) S += rn[r].su;
// transpose.scm set-mmt-marginals: count part (index 1) per row; total = sum over rows.
double *m = calloc(rows.n, sizeof(double)), T = 0;
for (size_t k = 0; k < nents; k++) m[ents[k].r] += cn[ents[k].c].l1 * fabs(ents[k].n);
for (size_t r = 0; r < rows.n; r++) T += m[r];
double Q = -log2(T / (S * S));
fprintf(stderr, "S (total-sum-left) %.17g T (MM^T total count) %.17g Q %.17g\n", S, T, Q);
// rank-words.
long *order = malloc(sizeof(long) * rows.n);
rank_count = malloc(sizeof(double) * rows.n); rank_supp = malloc(sizeof(double) * rows.n);
for (size_t r = 0; r < rows.n; r++) { order[r] = r; rank_count[r] = rn[r].l1; rank_supp[r] = rn[r].l0; }
qsort(order, rows.n, sizeof(long), by_rank);
// Row start offsets (entries are sorted by row, then column).
size_t *start = calloc(rows.n + 1, sizeof(size_t));
for (size_t k = 0; k < nents; k++) start[ents[k].r + 1]++;
for (size_t r = 0; r < rows.n; r++) start[r + 1] += start[r];
int nr = nrank < (int) rows.n ? nrank : (int) rows.n;
sim_t *sims = malloc(sizeof(sim_t) * nr * (nr + 1) / 2);
size_t ns = 0;
for (int i = 0; i < nr; i++) for (int j = i; j < nr; j++) {
long a = order[i], b = order[j];
double dot = 0;
size_t p = start[a], q = start[b];
while (p < start[a + 1] && q < start[b + 1]) {
if (ents[p].c == ents[q].c) {
if (!as_run || !strncmp(cols.keys[ents[p].c], "(ShapeLink", 10)) dot += ents[p].n * ents[q].n;
p++; q++;
}
else if (ents[p].c < ents[q].c) p++;
else q++;
}
double fmi = dot > 0 ? log2(dot * T / (m[a] * m[b])) : -INFINITY;
double rmi = fmi + 0.5 * log2((m[a] / T) * (m[b] / T)) + Q;
sims[ns++] = (sim_t) { a, b, fmi, rmi };
printf("%s\t%s\t%.17g\t%.17g\n", rows.keys[a], rows.keys[b], fmi, rmi);
}
// The seed list: fmi > 4.0, a != b, by ranked-MI (agglo-mi-rank.scm:628-639, mi-similarity.scm:70-116).
qsort(sims, ns, sizeof(sim_t), by_rmi);
int shown = 0;
for (size_t k = 0; k < ns && shown < top; k++)
if (sims[k].a != sims[k].b && sims[k].fmi > 4.0) {
fprintf(stderr, "ranked-MI = %.4f MI = %.4f (`%s`, `%s`)\n", sims[k].rmi, sims[k].fmi, rows.keys[sims[k].a], rows.keys[sims[k].b]);
shown++;
}
fprintf(stderr, "top %d ranked words:", nr < 12 ? nr : 12);
for (int i = 0; i < nr && i < 12; i++) fprintf(stderr, " `%s`", rows.keys[order[i]]);
fprintf(stderr, "\n");
return 0;
}