# The marriage model of Gottman, Murray, Swanson, Tyson & Swanson (2002). Each speaking turn in a
# conflict discussion gets a score (coded positivity minus negativity). For wife W and husband H:
#   W[t+1] = a + r1 W[t] + I_HW(H[t])        H[t+1] = b + r2 H[t] + I_WH(W[t+1])
# a/(1 - r1) is the wife's "uninfluenced set point" (where she would settle on her own), r1 her
# emotional inertia, I_HW the husband's influence on her. Steady states are where the two
# null-clines W = (a + I_HW(H))/(1 - r1) and H = (b + I_WH(W))/(1 - r2) cross. Two forms of influence
# function from the book: bilinear (a slope for positive and a steeper one for negative scores) and
# ogive (influence saturates). Parameters are illustrative, not estimates for any couple.
import math

def bilinear(pos, neg):
    return lambda x: pos * x if x >= 0 else neg * x

def ogive(pos, neg, s=1.5):
    return lambda x: pos * math.tanh(x / s) if x >= 0 else neg * math.tanh(x / s)

def run(a, r1, b, r2, IHW, IWH, W0, H0, n=400):
    W, H = W0, H0
    for _ in range(n):
        W = a + r1 * W + IHW(H)
        H = b + r2 * H + IWH(W)
    return round(W, 2), round(H, 2)

def attractors(a, r1, b, r2, IHW, IWH):
    found = {}
    for W0 in range(-12, 13, 2):
        for H0 in range(-12, 13, 2):
            st = run(a, r1, b, r2, IHW, IWH, W0, H0)
            if all(abs(v) < 1e6 for v in st):
                found[st] = found.get(st, 0) + 1
    return found

cases = [
    ("bilinear, negativity twice as influential", 1.0, 0.5, 0.8, 0.5, bilinear(0.2, 0.4), bilinear(0.2, 0.4)),
    ("ogive, mild negative influence",            1.0, 0.5, 0.8, 0.5, ogive(1.0, 1.5), ogive(1.0, 1.5)),
    ("ogive, strong negative influence",          1.0, 0.5, 0.8, 0.5, ogive(1.0, 4.0), ogive(1.0, 4.0)),
    ("same, but set points raised (a = b = 2)",   2.0, 0.5, 2.0, 0.5, ogive(1.0, 4.0), ogive(1.0, 4.0)),
]
for label, a, r1, b, r2, IHW, IWH in cases:
    att = attractors(a, r1, b, r2, IHW, IWH)
    total = sum(att.values())
    desc = "; ".join(f"(W, H) = ({w:+.2f}, {h:+.2f}) from {k / total:.0%} of starts" for (w, h), k in sorted(att.items()))
    print(f"{label}:\n  uninfluenced set points {a / (1 - r1):+.1f}, {b / (1 - r2):+.1f} -> {desc}")
# bilinear, negativity twice as influential:
#   uninfluenced set points +2.0, +1.6 -> (W, H) = (+3.14, +2.86) from 100% of starts
# ogive, mild negative influence:
#   uninfluenced set points +2.0, +1.6 -> (W, H) = (+3.97, +3.58) from 100% of starts
# ogive, strong negative influence:
#   uninfluenced set points +2.0, +1.6 -> (W, H) = (-6.00, -6.39) from 50% of starts; (W, H) = (+3.97, +3.58) from 50% of starts
# same, but set points raised (a = b = 2):
#   uninfluenced set points +4.0, +4.0 -> (W, H) = (-3.91, -3.91) from 34% of starts; (W, H) = (+6.00, +6.00) from 66% of starts

# Why "predicts divorce with 90 % accuracy" needs cross-validation (the point made by Heyman & Slep
# 2001). Synthetic data: 60 couples, 15 of whom later separate. 30 coded interaction measures; only
# 2 are truly (weakly, d = 0.5) related to separating, the other 28 are noise. The analyst picks the
# 5 measures that best separate the groups and classifies each couple by a simple score built from
# them. Accuracy is computed (a) on the same couples used to pick the measures and (b) by
# leave-one-out, picking the measures again without the couple being predicted.
import math, random, statistics as st

def make(rng, n=60, n_sep=15, k=30, real=2, d=0.5):
    y = [1] * n_sep + [0] * (n - n_sep)
    X = [[rng.gauss(d if (j < real and yi) else 0.0, 1) for j in range(k)] for yi in y]
    return X, y

def fit(X, y, top=5):
    g1 = [x for x, yi in zip(X, y) if yi]; g0 = [x for x, yi in zip(X, y) if not yi]
    k = len(X[0]); stats = []
    for j in range(k):
        a, b = [x[j] for x in g1], [x[j] for x in g0]
        sp = math.sqrt((st.variance(a) + st.variance(b)) / 2)
        stats.append(((st.mean(a) - st.mean(b)) / sp, j, st.mean(a), st.mean(b), sp))
    best = sorted(stats, key=lambda s: -abs(s[0]))[:top]
    def score(x): return sum(dz * (x[j] - (m1 + m0) / 2) / sp for dz, j, m1, m0, sp in best)
    scores = sorted(score(x) for x in X)
    cut = scores[len(scores) - sum(y)]          # call the top-scoring couples "will separate"
    return lambda x: 1 if score(x) >= cut else 0, [j for _, j, *_ in best]

rng = random.Random(7)
ins, loo, real_picked = [], [], []
for rep in range(60):
    X, y = make(rng)
    pred, picked = fit(X, y)
    ins.append(sum(pred(x) == yi for x, yi in zip(X, y)) / len(y))
    real_picked.append(sum(j < 2 for j in picked))
    hits = 0
    for i in range(len(y)):
        p, _ = fit(X[:i] + X[i + 1:], y[:i] + y[i + 1:])
        hits += p(X[i]) == y[i]
    loo.append(hits / len(y))
print(f"guessing 'nobody separates' is right {45 / 60:.0%} of the time")
print(f"accuracy on the couples used to choose the measures: {st.mean(ins):.0%}")
print(f"leave-one-out accuracy (measures re-chosen without the predicted couple): {st.mean(loo):.0%}")
print(f"of the 5 measures chosen, on average {st.mean(real_picked):.1f} are among the 2 that are real")
# guessing 'nobody separates' is right 75% of the time
# accuracy on the couples used to choose the measures: 80%
# leave-one-out accuracy (measures re-chosen without the predicted couple): 66%
# of the 5 measures chosen, on average 1.2 are among the 2 that are real