/* model.c — see model.h. */
#include "model.h"
#include <math.h>
#include <stdlib.h>
#include <string.h>
/* A relation forms only when its joint distribution diverges from the factored
* marginals by MORE than sampling noise would produce at its observed count.
* Gate = G-test significance, not a raw MI threshold. */
#define MINOBS 5 /* need a few observations for the test to be valid */
#define CHI2_CRIT 10.828 /* chi-squared, df = NUM_CLASSES-1 = 1, p < 0.001 */
void rel_init(Rel *r, int nf0, int nf1, double alpha) {
r->alpha = alpha; r->nf0 = nf0; r->nf1 = nf1;
r->m0 = calloc((size_t)nf0 * NUM_CLASSES, sizeof *r->m0);
r->s0 = calloc((size_t)nf0, sizeof *r->s0);
r->m1 = calloc((size_t)nf1 * NUM_CLASSES, sizeof *r->m1);
r->s1 = calloc((size_t)nf1, sizeof *r->s1);
r->j = calloc((size_t)nf0 * nf1 * NUM_CLASSES, sizeof *r->j);
r->sj = calloc((size_t)nf0 * nf1, sizeof *r->sj);
}
void rel_free(Rel *r) {
free(r->m0); free(r->s0); free(r->m1); free(r->s1); free(r->j); free(r->sj);
memset(r, 0, sizeof *r);
}
void rel_update(Rel *r, int f0, int f1, int klass) {
r->m0[(size_t)f0 * NUM_CLASSES + klass]++; r->s0[f0]++;
r->m1[(size_t)f1 * NUM_CLASSES + klass]++; r->s1[f1]++;
size_t p = (size_t)f0 * r->nf1 + f1;
r->j[p * NUM_CLASSES + klass]++; r->sj[p]++;
}
static void dirichlet(const long *row, long seen, double a, double out[NUM_CLASSES]) {
double denom = (double)seen + a * NUM_CLASSES;
for (int c = 0; c < NUM_CLASSES; c++) out[c] = ((double)row[c] + a) / denom;
}
void pred_ext(const Rel *r, int f0, double out[NUM_CLASSES]) {
dirichlet(&r->m0[(size_t)f0 * NUM_CLASSES], r->s0[f0], r->alpha, out);
}
void pred_dir(const Rel *r, int f1, double out[NUM_CLASSES]) {
dirichlet(&r->m1[(size_t)f1 * NUM_CLASSES], r->s1[f1], r->alpha, out);
}
/* the "no relation" prediction: average of the two marginals. */
void pred_marg(const Rel *r, int f0, int f1, double out[NUM_CLASSES]) {
double e[NUM_CLASSES], d[NUM_CLASSES];
pred_ext(r, f0, e); pred_dir(r, f1, d);
for (int c = 0; c < NUM_CLASSES; c++) out[c] = 0.5 * (e[c] + d[c]);
}
static void pred_joint_raw(const Rel *r, int f0, int f1, double out[NUM_CLASSES]) {
size_t p = (size_t)f0 * r->nf1 + f1;
dirichlet(&r->j[p * NUM_CLASSES], r->sj[p], r->alpha, out);
}
/* KL( joint || marginal ) — how much the pair's joint distribution cannot be
* reconstructed from its parts. This is the emergence / interaction signal. */
double rel_gain(const Rel *r, int f0, int f1) {
double j[NUM_CLASSES], m[NUM_CLASSES], kl = 0.0;
pred_joint_raw(r, f0, f1, j);
pred_marg(r, f0, f1, m);
for (int c = 0; c < NUM_CLASSES; c++)
if (j[c] > 1e-12) kl += j[c] * log(j[c] / m[c]);
return kl;
}
/* G-test statistic = 2 * sum_c O_c ln(O_c / E_c), with E_c = N * P_marginal(c).
* Asymptotically chi-squared(df = NUM_CLASSES-1) under "no relation". Because
* it multiplies KL by N, a thinly-observed pair cannot pass however peaked it
* looks — this is the small-sample-bias correction the raw MI threshold lacked. */
static double rel_G(const Rel *r, int f0, int f1) {
size_t p = (size_t)f0 * r->nf1 + f1;
long N = r->sj[p];
if (N == 0) return 0.0;
double m[NUM_CLASSES];
pred_marg(r, f0, f1, m);
const long *O = &r->j[p * NUM_CLASSES];
double G = 0.0;
for (int c = 0; c < NUM_CLASSES; c++)
if (O[c] > 0) G += (double)O[c] * log((double)O[c] / ((double)N * m[c]));
return 2.0 * G;
}
int rel_emerged(const Rel *r, int f0, int f1) {
return pair_seen(r, f0, f1) >= MINOBS && rel_G(r, f0, f1) >= CHI2_CRIT;
}
/* gated prediction: use the joint iff the relation has emerged, else marginals. */
void pred_rel(const Rel *r, int f0, int f1, double out[NUM_CLASSES]) {
if (rel_emerged(r, f0, f1)) pred_joint_raw(r, f0, f1, out);
else pred_marg(r, f0, f1, out);
}
int argmax(const double p[NUM_CLASSES]) {
int b = 0;
for (int c = 1; c < NUM_CLASSES; c++) if (p[c] > p[b]) b = c;
return b;
}
double infogain_pair(const Rel *r, int f0, int f1) {
size_t p = (size_t)f0 * r->nf1 + f1;
return 1.0 / ((double)r->sj[p] + r->alpha * NUM_CLASSES);
}
long pair_seen(const Rel *r, int f0, int f1) {
return r->sj[(size_t)f0 * r->nf1 + f1];
}