| … | … | @@ -0,0 +1,226 @@ |
| 1 | // Stage E, step 1: Linas's covering matrix and the initial grammatical-MI similarities (see STAGE_E_SPEC.md §0-1). |
| 2 | // |
| 3 | // * Sections from the counts file (output of `count`), exploded into CrossSections as in |
| 4 | // learn shape-vec.scm explode-sections: Section (g, c1..cn) with count N gives, for each i, |
| 5 | // CrossSection [ci.word, <g, c1 .. ($, ci.dir) .. cn>] with count N (set, not added: a duplicate |
| 6 | // (row, column) is an error here). |
| 7 | // * Rows = words (germs of Sections and of CrossSections); columns = ConnectorSeqs and Shapes. |
| 8 | // * Marginals (matrix support.scm) and MM^T (transpose.scm) on that matrix, as in csetmarg.c. |
| 9 | // * rank-words (similarity.scm:17-42): right-count desc, ties right-support desc. |
| 10 | // * Similarities for all pairs i <= j among the top NRANK ranked words (agglo-mi-rank.scm:535, |
| 11 | // loop-upper-diagonal): fmi = log2(sum_d N(a,d)N(b,d) * T / (m(a) m(b))) and |
| 12 | // rmi = fmi + 0.5 log2(m(a)/T * m(b)/T) + Q, Q = -log2(T / S^2) (symmetric-mi.scm:170-264, |
| 13 | // mi-similarity.scm:49-53); m = MM^T count of the row, T = MM^T total count, S = total-sum-left. |
| 14 | // --as-run: reproduce Linas's stack as it runs today, where the dot product sees Shape columns only. The |
| 15 | // covering object's in-right-basis? is built from a Query on (make-pair $x $y), and direct-sum's make-pair |
| 16 | // with two variables always builds a CrossSection (direct-sum.scm:234-243, make-left-variable-pair tests |
| 17 | // the right atom's type, which is a variable), so add-fast-math (fold-api.scm:351-353) drops every |
| 18 | // ConnectorSeq column. Marginals do not go through that filter and still include Sections. |
| 19 | // Output: one line per similarity: <word a>\t<word b>\t<fmi>\t<rmi>. With --top K, the K best pairs with |
| 20 | // fmi > 4 (the seed list of agglo-mi-rank.scm:628-639) are printed to stderr in Linas's log format. |
| 21 | // |
| 22 | // cover <counts.tsv> [--as-run] [--nrank 50] [--top 12] > sims.tsv |
| 23 | |
| 24 | #include <math.h> |
| 25 | #include <stdint.h> |
| 26 | #include <stdio.h> |
| 27 | #include <stdlib.h> |
| 28 | #include <string.h> |
| 29 | |
| 30 | // ---------------------------------------------------------------- interning |
| 31 | |
| 32 | typedef struct { char **keys; size_t n, cap, mask; long *tab; } intern_t; |
| 33 | |
| 34 | static uint64_t fnv (const char *s) { |
| 35 | uint64_t h = 1469598103934665603ULL; |
| 36 | while (*s) h = (h ^ (unsigned char) *s++) * 1099511628211ULL; |
| 37 | return h; |
| 38 | } |
| 39 | |
| 40 | static void rehash (intern_t *t) { |
| 41 | t->mask = t->mask ? t->mask * 2 + 1 : 4095; |
| 42 | free(t->tab); |
| 43 | t->tab = malloc(sizeof(long) * (t->mask + 1)); |
| 44 | for (size_t i = 0; i <= t->mask; i++) t->tab[i] = -1; |
| 45 | for (size_t k = 0; k < t->n; k++) { |
| 46 | size_t i = fnv(t->keys[k]) & t->mask; |
| 47 | while (t->tab[i] >= 0) i = (i + 1) & t->mask; |
| 48 | t->tab[i] = k; |
| 49 | } |
| 50 | } |
| 51 | |
| 52 | static long intern (intern_t *t, const char *s) { |
| 53 | if (!t->mask) rehash(t); |
| 54 | size_t i = fnv(s) & t->mask; |
| 55 | while (t->tab[i] >= 0) { |
| 56 | if (!strcmp(t->keys[t->tab[i]], s)) return t->tab[i]; |
| 57 | i = (i + 1) & t->mask; |
| 58 | } |
| 59 | if (t->n == t->cap) { t->cap = t->cap ? t->cap * 2 : 4096; t->keys = realloc(t->keys, sizeof(char *) * t->cap); } |
| 60 | t->keys[t->n] = strdup(s); |
| 61 | t->tab[i] = t->n++; |
| 62 | if (t->n * 2 > t->mask) rehash(t); |
| 63 | return t->n - 1; |
| 64 | } |
| 65 | |
| 66 | // ---------------------------------------------------------------- matrix |
| 67 | |
| 68 | typedef struct { long r, c; double n; } entry_t; |
| 69 | typedef struct { double l0, l1, l2, lq, su; } norms_t; |
| 70 | |
| 71 | static intern_t rows, cols; |
| 72 | static entry_t *ents; |
| 73 | static size_t nents, capents; |
| 74 | |
| 75 | static void add_entry (long r, long c, double n) { |
| 76 | if (nents == capents) { capents = capents ? capents * 2 : 1 << 18; ents = realloc(ents, sizeof(entry_t) * capents); } |
| 77 | ents[nents++] = (entry_t) { r, c, n }; |
| 78 | } |
| 79 | |
| 80 | 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; } |
| 81 | |
| 82 | static int by_row_col (const void *a, const void *b) { |
| 83 | const entry_t *x = a, *y = b; |
| 84 | if (x->r != y->r) return x->r < y->r ? -1 : 1; |
| 85 | return (x->c > y->c) - (x->c < y->c); |
| 86 | } |
| 87 | |
| 88 | // Parse the "(Connector (WordNode \"w\") (SexNode \"d\"))" items of a ConnectorSeq string. |
| 89 | static int connectors (const char *seq, char words[][1024], char *dirs, int max) { |
| 90 | int n = 0; |
| 91 | const char *p = seq; |
| 92 | while ((p = strstr(p, "(Connector (WordNode \"")) && n < max) { |
| 93 | p += 22; |
| 94 | const char *e = strstr(p, "\") (SexNode \""); |
| 95 | if (!e || e - p >= 1024) return -1; |
| 96 | memcpy(words[n], p, e - p); |
| 97 | words[n][e - p] = 0; |
| 98 | dirs[n] = e[13]; |
| 99 | p = e + 13; |
| 100 | n++; |
| 101 | } |
| 102 | return n; |
| 103 | } |
| 104 | |
| 105 | static double *rank_count, *rank_supp; |
| 106 | |
| 107 | static int by_rank (const void *a, const void *b) { |
| 108 | long x = *(const long *) a, y = *(const long *) b; |
| 109 | if (rank_count[x] != rank_count[y]) return rank_count[x] > rank_count[y] ? -1 : 1; |
| 110 | if (rank_supp[x] != rank_supp[y]) return rank_supp[x] > rank_supp[y] ? -1 : 1; |
| 111 | return (x > y) - (x < y); |
| 112 | } |
| 113 | |
| 114 | typedef struct { long a, b; double fmi, rmi; } sim_t; |
| 115 | |
| 116 | static int by_rmi (const void *x, const void *y) { |
| 117 | const sim_t *p = x, *q = y; |
| 118 | return (p->rmi < q->rmi) - (p->rmi > q->rmi); |
| 119 | } |
| 120 | |
| 121 | int main (int argc, char **argv) { |
| 122 | if (argc < 2) { fprintf(stderr, "usage: cover <counts.tsv> [--as-run] [--nrank 50] [--top 12]\n"); return 1; } |
| 123 | int nrank = 50, top = 12, as_run = 0; |
| 124 | for (int i = 2; i < argc; i++) { |
| 125 | if (!strcmp(argv[i], "--as-run")) as_run = 1; |
| 126 | else if (!strcmp(argv[i], "--nrank") && i + 1 < argc) nrank = atoi(argv[++i]); |
| 127 | else if (!strcmp(argv[i], "--top") && i + 1 < argc) top = atoi(argv[++i]); |
| 128 | } |
| 129 | FILE *f = fopen(argv[1], "r"); |
| 130 | if (!f) { perror(argv[1]); return 1; } |
| 131 | static char line[1 << 16], seq[1 << 16], shape[1 << 16]; |
| 132 | static char cw[64][1024]; |
| 133 | char cd[64]; |
| 134 | const char *pre = "(Section (WordNode \""; |
| 135 | long nsect = 0; |
| 136 | while (fgets(line, sizeof line, f)) { |
| 137 | if (strncmp(line, pre, strlen(pre))) continue; |
| 138 | char *g0 = line + strlen(pre), *g1 = strstr(g0, "\") (ConnectorSeq"); |
| 139 | char *t1 = strchr(line, '\t'), *t2 = t1 ? strchr(t1 + 1, '\t') : NULL; |
| 140 | if (!g1 || !t1 || !t2) continue; |
| 141 | double a, b, n; |
| 142 | if (3 != sscanf(t2 + 1, "%lf %lf %lf", &a, &b, &n)) { fprintf(stderr, "bad count: %s", line); return 1; } |
| 143 | *g1 = 0; |
| 144 | t1[-1] = 0; |
| 145 | strcpy(seq, g1 + 3); |
| 146 | int nc = connectors(seq, cw, cd, 64); |
| 147 | if (nc <= 0) { fprintf(stderr, "bad section: %s\n", seq); return 1; } |
| 148 | add_entry(intern(&rows, g0), intern(&cols, seq), n); |
| 149 | nsect++; |
| 150 | for (int i = 0; i < nc; i++) { // explode-sections: hole at connector i |
| 151 | int k = snprintf(shape, sizeof shape, "(ShapeLink (WordNode \"%s\")", g0); |
| 152 | for (int j = 0; j < nc; j++) { |
| 153 | if (j == i) k += snprintf(shape + k, sizeof shape - k, " (Connector (VariableNode \"$connector-word\") (SexNode \"%c\"))", cd[j]); |
| 154 | else k += snprintf(shape + k, sizeof shape - k, " (Connector (WordNode \"%s\") (SexNode \"%c\"))", cw[j], cd[j]); |
| 155 | } |
| 156 | snprintf(shape + k, sizeof shape - k, ")"); |
| 157 | add_entry(intern(&rows, cw[i]), intern(&cols, shape), n); |
| 158 | } |
| 159 | } |
| 160 | fclose(f); |
| 161 | qsort(ents, nents, sizeof(entry_t), by_row_col); |
| 162 | for (size_t k = 1; k < nents; k++) |
| 163 | if (ents[k].r == ents[k - 1].r && ents[k].c == ents[k - 1].c) { |
| 164 | fprintf(stderr, "duplicate (row, column): %s %s\n", rows.keys[ents[k].r], cols.keys[ents[k].c]); |
| 165 | return 2; |
| 166 | } |
| 167 | fprintf(stderr, "sections %ld, crosses %zu, rows %zu, columns %zu\n", nsect, nents - nsect, rows.n, cols.n); |
| 168 | |
| 169 | // support.scm marginals: right = per row, left = per column. |
| 170 | norms_t *rn = calloc(rows.n, sizeof(norms_t)), *cn = calloc(cols.n, sizeof(norms_t)); |
| 171 | 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); } |
| 172 | double S = 0; // total-sum-left = sum over rows of right-sum (support.scm do-left-totals) |
| 173 | for (size_t r = 0; r < rows.n; r++) S += rn[r].su; |
| 174 | |
| 175 | // transpose.scm set-mmt-marginals: count part (index 1) per row; total = sum over rows. |
| 176 | double *m = calloc(rows.n, sizeof(double)), T = 0; |
| 177 | for (size_t k = 0; k < nents; k++) m[ents[k].r] += cn[ents[k].c].l1 * fabs(ents[k].n); |
| 178 | for (size_t r = 0; r < rows.n; r++) T += m[r]; |
| 179 | double Q = -log2(T / (S * S)); |
| 180 | fprintf(stderr, "S (total-sum-left) %.17g T (MM^T total count) %.17g Q %.17g\n", S, T, Q); |
| 181 | |
| 182 | // rank-words. |
| 183 | long *order = malloc(sizeof(long) * rows.n); |
| 184 | rank_count = malloc(sizeof(double) * rows.n); rank_supp = malloc(sizeof(double) * rows.n); |
| 185 | for (size_t r = 0; r < rows.n; r++) { order[r] = r; rank_count[r] = rn[r].l1; rank_supp[r] = rn[r].l0; } |
| 186 | qsort(order, rows.n, sizeof(long), by_rank); |
| 187 | |
| 188 | // Row start offsets (entries are sorted by row, then column). |
| 189 | size_t *start = calloc(rows.n + 1, sizeof(size_t)); |
| 190 | for (size_t k = 0; k < nents; k++) start[ents[k].r + 1]++; |
| 191 | for (size_t r = 0; r < rows.n; r++) start[r + 1] += start[r]; |
| 192 | |
| 193 | int nr = nrank < (int) rows.n ? nrank : (int) rows.n; |
| 194 | sim_t *sims = malloc(sizeof(sim_t) * nr * (nr + 1) / 2); |
| 195 | size_t ns = 0; |
| 196 | for (int i = 0; i < nr; i++) for (int j = i; j < nr; j++) { |
| 197 | long a = order[i], b = order[j]; |
| 198 | double dot = 0; |
| 199 | size_t p = start[a], q = start[b]; |
| 200 | while (p < start[a + 1] && q < start[b + 1]) { |
| 201 | if (ents[p].c == ents[q].c) { |
| 202 | if (!as_run || !strncmp(cols.keys[ents[p].c], "(ShapeLink", 10)) dot += ents[p].n * ents[q].n; |
| 203 | p++; q++; |
| 204 | } |
| 205 | else if (ents[p].c < ents[q].c) p++; |
| 206 | else q++; |
| 207 | } |
| 208 | double fmi = dot > 0 ? log2(dot * T / (m[a] * m[b])) : -INFINITY; |
| 209 | double rmi = fmi + 0.5 * log2((m[a] / T) * (m[b] / T)) + Q; |
| 210 | sims[ns++] = (sim_t) { a, b, fmi, rmi }; |
| 211 | printf("%s\t%s\t%.17g\t%.17g\n", rows.keys[a], rows.keys[b], fmi, rmi); |
| 212 | } |
| 213 | |
| 214 | // The seed list: fmi > 4.0, a != b, by ranked-MI (agglo-mi-rank.scm:628-639, mi-similarity.scm:70-116). |
| 215 | qsort(sims, ns, sizeof(sim_t), by_rmi); |
| 216 | int shown = 0; |
| 217 | for (size_t k = 0; k < ns && shown < top; k++) |
| 218 | if (sims[k].a != sims[k].b && sims[k].fmi > 4.0) { |
| 219 | 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]); |
| 220 | shown++; |
| 221 | } |
| 222 | fprintf(stderr, "top %d ranked words:", nr < 12 ? nr : 12); |
| 223 | for (int i = 0; i < nr && i < 12; i++) fprintf(stderr, " `%s`", rows.keys[order[i]]); |
| 224 | fprintf(stderr, "\n"); |
| 225 | return 0; |
| 226 | } |