"""Numerical checks of the derivations in Bishop (2006), Chapter 9 (Mixture Models and EM).

Each check prints PASS or FAIL with the numbers behind it; the last line counts them.
Run from this folder: python3 checks_ch09.py (a few seconds). The algorithms checked are
the ones in ch09_mixture_em.py, which also produce the figures.
"""
import itertools

import numpy as np

import common
import ch09_mixture_em as m

RESULTS = []


def check(name, ok, detail=""):
    RESULTS.append(bool(ok))
    print("%s  %s%s" % ("PASS" if ok else "FAIL", name, ("  [" + detail + "]") if detail else ""))


def num_grad(f, x, h=1e-6):
    """Central finite-difference gradient of a scalar function of an array."""
    x = np.array(x, dtype=float)
    g = np.zeros_like(x)
    for i in np.ndindex(x.shape):
        e = np.zeros_like(x)
        e[i] = h
        g[i] = (f(x + e) - f(x - e)) / (2 * h)
    return g


X = common.old_faithful(standardize=True)
N, D = X.shape


def random_gmm(rng, K):
    """Random but sensible mixture parameters for the standardised Old Faithful data."""
    mu = X[rng.choice(N, K, replace=False)]
    A = rng.normal(size=(K, D, D)) * 0.3
    Sigma = np.einsum("kij,klj->kil", A, A) + 0.2 * np.eye(D)
    pi = rng.dirichlet(np.ones(K) * 3)
    return pi, mu, Sigma


def Q(pi, mu, Sigma, gamma):
    """Expected complete-data log likelihood (9.40): sum_n sum_k gamma_nk {ln pi_k + ln N(x_n|mu_k, Sigma_k)}."""
    return float((gamma * m.gmm_log_terms(X, pi, mu, Sigma)).sum())


# ----------------------------------------------------------------------------- K-means
print("== K-means (Section 9.1)")
res = m.kmeans(X, m.BOOK_CENTRES)
J = np.array(res["J"])
check("J (9.1) never increases over the half steps from the book's start (Fig. 9.2)", np.all(np.diff(J) <= 1e-12),
      "J = " + ", ".join("%.2f" % v for v in J))
check("book's start: J after the first E step is about 1050 (Fig. 9.2)", 1000 < J[0] < 1100, "J = %.2f" % J[0])
rng = np.random.default_rng(1)
worst, iters = 0.0, []
for _ in range(200):
    K = rng.integers(2, 7)
    r = m.kmeans(X, X[rng.choice(N, K, replace=False)])
    worst = max(worst, float(np.max(np.diff(r["J"]))))
    iters.append(r["n_m_steps"])
check("J never increases and K-means stops after finitely many steps (Exercise 9.1), 200 random starts, K = 2..6",
      worst <= 1e-10 and max(iters) < 500, "largest increase %.1e, most M steps %d" % (worst, max(iters)))
mu = X[rng.choice(N, 3, replace=False)]
r = m.kmeans_e_step(X, mu)
mu_new = m.kmeans_m_step(X, r, mu)
g = num_grad(lambda v: m.distortion(X, r, v.reshape(3, D)), mu_new.ravel())
check("the M step (9.4) is a stationary point of J for fixed r_nk (9.3)", np.abs(g).max() < 1e-5, "max |dJ/dmu| = %.1e" % np.abs(g).max())

# ----------------------------------------------------------------------------- Gaussian mixtures
print("== EM for Gaussian mixtures (Section 9.2)")
pi, mu, Sigma = random_gmm(np.random.default_rng(2), 3)
gamma, _ = m.gmm_e_step(X, pi, mu, Sigma)
for k in range(3):
    Sk_inv = np.linalg.inv(Sigma[k])
    f = lambda v: m.gmm_loglik(X, pi, np.vstack([mu[:k], v[None], mu[k + 1:]]), Sigma)
    g_num = num_grad(f, mu[k])
    g_formula = Sk_inv @ (gamma[:, k, None] * (X - mu[k])).sum(0)
    check("gradient of ln p(X) w.r.t. mu_%d equals sum_n gamma_nk Sigma_k^-1 (x_n - mu_k), (9.16) with the erratum" % (k + 1),
          np.allclose(g_num, g_formula, rtol=1e-6, atol=1e-6), "max diff %.1e" % np.abs(g_num - g_formula).max())

pi_new, mu_new, Sigma_new = m.gmm_m_step(X, gamma)
g_mu = num_grad(lambda v: Q(pi_new, v.reshape(mu_new.shape), Sigma_new, gamma), mu_new.ravel())
check("(9.17): the new means are a stationary point of Q (9.40)", np.abs(g_mu).max() < 1e-5, "max |dQ/dmu| = %.1e" % np.abs(g_mu).max())


def Q_sym(v, k):
    """Q as a function of the 3 free entries (a, b, c) of the symmetric Sigma_k = [[a, b], [b, c]]."""
    S = Sigma_new.copy()
    S[k] = [[v[0], v[1]], [v[1], v[2]]]
    return Q(pi_new, mu_new, S, gamma)


mx = max(np.abs(num_grad(lambda v: Q_sym(v, k), [Sigma_new[k][0, 0], Sigma_new[k][0, 1], Sigma_new[k][1, 1]])).max()
         for k in range(3))
check("(9.19): the new covariances are a stationary point of Q (symmetric perturbations)", mx < 1e-4, "max |dQ/dSigma| = %.1e" % mx)
g_pi = num_grad(lambda v: Q(v, mu_new, Sigma_new, gamma), pi_new)
check("(9.21)-(9.22): with pi_k = N_k/N, dQ/dpi_k = N for every k, i.e. the Lagrange condition holds with lambda = -N",
      np.allclose(g_pi, N, rtol=1e-6), "dQ/dpi = %s" % np.round(g_pi, 4).tolist())
rng = np.random.default_rng(3)
q0, worst = Q(pi_new, mu_new, Sigma_new, gamma), -np.inf
for _ in range(500):
    p2 = rng.dirichlet(pi_new * 200)
    m2 = mu_new + 0.05 * rng.normal(size=mu_new.shape)
    A = 0.05 * rng.normal(size=(3, D, D))
    S2 = np.array([Sigma_new[k] + A[k] @ A[k].T - 0.5 * (A[k] + A[k].T) * 0.2 for k in range(3)])
    if np.all(np.linalg.eigvalsh(S2) > 0):
        worst = max(worst, Q(p2, m2, S2, gamma) - q0)
check("the M step is a maximum of Q: 500 random nearby parameter values all give a smaller Q", worst < 0, "largest Q - Q* = %.3g" % worst)

worst, n_runs = 0.0, 0
for seed in range(40):
    rng = np.random.default_rng(100 + seed)
    K = 2 + seed % 5
    pi, mu, Sigma = random_gmm(rng, K)
    h = m.em_gmm(X, pi, mu, Sigma, n_iter=150, var_floor=1e-3)
    worst = min(worst, float(np.min(np.diff([s["loglik"] for s in h]))))
    n_runs += 1
check("EM never decreases ln p(X) (Section 9.4): %d runs on Old Faithful, K = 2..6, 150 iterations each" % n_runs,
      worst > -1e-8, "largest decrease %.1e" % -worst)
h = m.em_gmm(X, [0.5, 0.5], m.BOOK_CENTRES, np.stack([m.BOOK_COV * np.eye(2)] * 2), n_iter=60)
ll = np.array([s["loglik"] for s in h])
check("Fig. 9.8 run: ln p rises from %.2f to %.4f and is flat (to 1e-8) from L = 20 on" % (ll[0], ll[-1]),
      np.all(np.diff(ll) > -1e-10) and abs(ll[20] - ll[-1]) < 1e-8, "L = 1, 2, 5, 20: %s" % np.round(ll[[1, 2, 5, 20]], 4).tolist())

# ln p = L + KL for an arbitrary q (9.70)
pi, mu, Sigma = random_gmm(np.random.default_rng(4), 3)
T = m.gmm_log_terms(X, pi, mu, Sigma)
gamma, ll = m.gmm_e_step(X, pi, mu, Sigma)
qz = np.random.default_rng(5).dirichlet(np.ones(3), size=N)
L = float((qz * (T - np.log(qz))).sum())
KL = float((qz * (np.log(qz) - np.log(gamma))).sum())
check("ln p(X|theta) = L(q, theta) + KL(q||p) (9.70)-(9.72) for a random q", abs(L + KL - ll) < 1e-8 * abs(ll),
      "ln p = %.6f, L = %.6f, KL = %.6f" % (ll, L, KL))
check("KL(q||p) >= 0, so L is a lower bound; KL = 0 when q is the posterior", KL > 0 and abs(float((gamma * (T - np.log(gamma))).sum()) - ll) < 1e-8 * abs(ll))

# the posterior factorises over n (9.38), (9.75), and E[z_nk] = gamma(z_nk) (9.39): brute force for 4 points
Xs = X[[0, 1, 5, 9]]
K = 3
Ts = np.stack([np.log(pi[k]) + common.gauss_logpdf(Xs, mu[k], Sigma[k]) for k in range(K)], 1)
joint = {z: float(np.exp(sum(Ts[n, z[n]] for n in range(4)))) for z in itertools.product(range(K), repeat=4)}
Z = sum(joint.values())
gam_s = np.exp(Ts - common.logsumexp(Ts, axis=1)[:, None])
fact = max(abs(joint[z] / Z - np.prod([gam_s[n, z[n]] for n in range(4)])) for z in joint)
Ez = np.array([[sum(p for z, p in joint.items() if z[n] == k) / Z for k in range(K)] for n in range(4)])
check("p(Z|X) = prod_n p(z_n|x_n) by enumerating all 3^4 values of Z (9.38), (9.75)", fact < 1e-12, "max diff %.1e" % fact)
check("E[z_nk] under p(Z|X) equals the responsibility gamma(z_nk) (9.39)", np.abs(Ez - gam_s).max() < 1e-12)

# incremental EM: (9.78), (9.79)
gamma, _ = m.gmm_e_step(X, pi, mu, Sigma)
g2 = gamma.copy()
mm = 17
g2[mm] = np.random.default_rng(6).dirichlet(np.ones(3))
Nk_old = gamma.sum(0)
mu_old = (gamma.T @ X) / Nk_old[:, None]
Nk_new = Nk_old + g2[mm] - gamma[mm]
mu_inc = mu_old + ((g2[mm] - gamma[mm]) / Nk_new)[:, None] * (X[mm] - mu_old)
mu_batch = (g2.T @ X) / g2.sum(0)[:, None]
check("incremental EM: updating one point's responsibilities with (9.78), (9.79) equals the batch formula (9.17)",
      np.abs(mu_inc - mu_batch).max() < 1e-12 and np.allclose(Nk_new, g2.sum(0)))

# singularities (9.15), with the erratum sigma_j^D
x0 = X[40]
sig = np.array([1e-1, 1e-2, 1e-3, 1e-4])
dens = [common.gauss_logpdf(x0[None], x0, s ** 2 * np.eye(2))[0] for s in sig]
check("(9.15) with the erratum: ln N(x_n|x_n, sigma^2 I) = -(D/2) ln(2 pi) - D ln sigma in D = 2",
      np.allclose(dens, -np.log(2 * np.pi) - 2 * np.log(sig)))
# once sigma is below the distance to the nearest other point, each factor 10 adds D ln 10 = 4.605 to ln p
sig = 10.0 ** -np.arange(3, 10)
base = np.log(0.5) + common.gauss_logpdf(X, X.mean(0), np.cov(X.T, bias=True))
lls_mix = np.array([float(common.logsumexp(np.stack([base, np.log(0.5) + common.gauss_logpdf(X, x0, s ** 2 * np.eye(2))], 1), axis=1).sum()) for s in sig])
lls_one = np.array([float(common.gauss_logpdf(X, x0, s ** 2 * np.eye(2)).sum()) for s in sig])
check("a component collapsing onto x_n sends ln p(X) to +infinity (+2 ln 10 per factor 10 in sigma); a single Gaussian sends it to -infinity",
      np.allclose(np.diff(lls_mix), 2 * np.log(10), atol=1e-6) and np.all(np.diff(lls_one) < 0),
      "mixture, sigma = 1e-3 .. 1e-9: %s; single: %s" % (np.round(lls_mix, 2).tolist(), ["%.1e" % v for v in lls_one[:3]]))

# K-means as the eps -> 0 limit of EM with Sigma_k = eps I (Section 9.3.2)
mu0 = m.BOOK_CENTRES
r = m.kmeans_e_step(X, mu0)
R = np.eye(2)[r]
for eps in (1e-1, 1e-2, 1e-3):
    g_eps, _ = m.gmm_e_step(X, np.array([0.3, 0.7]), mu0, np.stack([eps * np.eye(2)] * 2))
    print("      eps = %g: max |gamma - r| = %.2e" % (eps, np.abs(g_eps - R).max()))
d2 = ((X[:, None, :] - mu0[None]) ** 2).sum(-1)
print("      the point nearest the bisector of the starting centres has |d_1 - d_2| = %.4f" % np.abs(d2[:, 0] - d2[:, 1]).min())
eps = 1e-6
g_eps, _ = m.gmm_e_step(X, np.array([0.3, 0.7]), mu0, np.stack([eps * np.eye(2)] * 2))
check("(9.42): as eps -> 0 the responsibilities become the K-means assignments r_nk (9.2), whatever pi",
      np.abs(g_eps - R).max() < 1e-12, "eps = 1e-6: max |gamma - r| = %.1e" % np.abs(g_eps - R).max())
_, mu_em, _ = m.gmm_m_step(X, g_eps, shared_var=eps)
check("... and the EM update of the means (9.17) becomes the K-means update (9.4)", np.abs(mu_em - m.kmeans_m_step(X, r, mu0)).max() < 1e-12)
hist = m.em_gmm(X, [0.5, 0.5], mu0, np.stack([eps * np.eye(2)] * 2), n_iter=6, shared_var=eps)
km = m.kmeans(X, mu0, record=True)
same = all(np.abs(hist[t]["mu"] - km["steps"][2 * t - 1][2]).max() < 1e-10 for t in range(1, 5))
check("EM with Sigma_k = 1e-6 I follows the K-means path of Fig. 9.1 step for step (4 iterations)", same)
J0 = m.distortion(X, r, mu0)
vals = []
for eps in (1e-2, 1e-4, 1e-6):
    g_eps, _ = m.gmm_e_step(X, np.array([0.5, 0.5]), mu0, np.stack([eps * np.eye(2)] * 2))
    EQ = Q(np.array([0.5, 0.5]), mu0, np.stack([eps * np.eye(2)] * 2), g_eps)
    vals.append(2 * eps * EQ)
check("(9.43): 2 eps E_Z[ln p(X, Z)] -> -J as eps -> 0 (with Exercise 9.11)", abs(vals[-1] + J0) < 1e-3 * J0,
      "2 eps E = %s, -J = %.4f" % (np.round(vals, 4).tolist(), -J0))

# ----------------------------------------------------------------------------- Bernoulli mixtures
print("== Mixtures of Bernoulli distributions (Section 9.3.3)")
rng = np.random.default_rng(7)
Kb, Db = 3, 4
pib = np.array([0.5, 0.3, 0.2])
mub = rng.uniform(0.05, 0.95, size=(Kb, Db))
mean_f = pib @ mub                                                              # (9.49)
cov_f = sum(pib[k] * (np.diag(mub[k] * (1 - mub[k])) + np.outer(mub[k], mub[k])) for k in range(Kb)) - np.outer(mean_f, mean_f)  # (9.50)
Mc = 400000
z = rng.choice(Kb, size=Mc, p=pib)
xs = (rng.uniform(size=(Mc, Db)) < mub[z]).astype(float)
se = np.sqrt(0.25 / Mc)
check("mean of a Bernoulli mixture (9.49) by Monte Carlo (400000 draws)", np.abs(xs.mean(0) - mean_f).max() < 5 * se,
      "max diff %.4f (5 se = %.4f)" % (np.abs(xs.mean(0) - mean_f).max(), 5 * se))
cov_mc = np.cov(xs.T, bias=True)
check("covariance of a Bernoulli mixture (9.50) by Monte Carlo: not diagonal, unlike (9.46)", np.abs(cov_mc - cov_f).max() < 0.005,
      "max diff %.4f; largest off-diagonal %.3f" % (np.abs(cov_mc - cov_f).max(), np.abs(cov_f - np.diag(np.diag(cov_f))).max()))
# EM: monotone, and the M step (9.59), (9.60) is a stationary point of (9.55)
Xb, yb = m.digits_234(60)
mu0b = np.random.default_rng(8).uniform(0.25, 0.75, size=(3, 784))
pi_b, mu_b, g_b, lls_b = m.bernoulli_mixture_em(Xb, np.full(3, 1 / 3), mu0b, n_iter=30)
check("EM for the Bernoulli mixture never decreases ln p(X|mu, pi) (9.51) (180 digits, 30 iterations)", np.all(np.diff(lls_b) > -1e-8),
      "from %.1f to %.1f" % (lls_b[0], lls_b[-1]))
check("ln p(X|mu, pi) <= 0 for binary data: the likelihood is bounded (Exercise 9.17), no singularities", max(lls_b) <= 0)
Xsm = (np.random.default_rng(9).uniform(size=(50, 6)) < 0.4).astype(float)
pi_s, mu_s = np.array([0.4, 0.6]), np.random.default_rng(10).uniform(0.2, 0.8, size=(2, 6))
gs, _ = m.bernoulli_e_step(Xsm, pi_s, mu_s)
pi_n, mu_n = m.bernoulli_m_step(Xsm, gs)
EQb = lambda mu_: float((gs * m.bernoulli_log_terms(Xsm, pi_n, mu_)).sum())
gq = num_grad(lambda v: EQb(v.reshape(mu_n.shape)), mu_n.ravel())
check("(9.59): mu_k = xbar_k is a stationary point of the expected complete-data log likelihood (9.55)", np.abs(gq).max() < 1e-5,
      "max |dE/dmu| = %.1e" % np.abs(gq).max())
gp = num_grad(lambda v: float((gs * m.bernoulli_log_terms(Xsm, v, mu_n)).sum()), pi_n)
check("(9.60): with pi_k = N_k/N the Lagrange condition dE/dpi_k = N holds for every k", np.allclose(gp, len(Xsm), rtol=1e-6))

# ----------------------------------------------------------------------------- EM for Bayesian linear regression (Section 9.3.4)
print("== EM for Bayesian linear regression (Section 9.3.4)")
rng = np.random.default_rng(11)
Nr = 25
xr = rng.uniform(0, 1, Nr)
t = np.sin(2 * np.pi * xr) + rng.normal(0, 0.3, Nr)
centres = np.linspace(0, 1, 9)
Phi = np.hstack([np.ones((Nr, 1)), np.exp(-(xr[:, None] - centres[None]) ** 2 / (2 * 0.1 ** 2))])
M = Phi.shape[1]
lam0 = np.linalg.eigvalsh(Phi.T @ Phi)


def posterior(alpha, beta):
    A = alpha * np.eye(M) + beta * Phi.T @ Phi                                 # (3.81)
    S = np.linalg.inv(A)
    return beta * S @ Phi.T @ t, S                                              # (3.53), (3.54)


def log_evidence(alpha, beta):
    """(3.86)."""
    mN, S = posterior(alpha, beta)
    E = beta / 2 * np.sum((t - Phi @ mN) ** 2) + alpha / 2 * mN @ mN
    return M / 2 * np.log(alpha) + Nr / 2 * np.log(beta) - E + 0.5 * np.linalg.slogdet(S)[1] - Nr / 2 * np.log(2 * np.pi)


a_ev, b_ev = 1.0, 1.0
for _ in range(5000):                                                           # evidence re-estimation (3.92), (3.95)
    mN, S = posterior(a_ev, b_ev)
    gam = np.sum(b_ev * lam0 / (a_ev + b_ev * lam0))                           # (3.91)
    a_ev, b_ev = gam / (mN @ mN), (Nr - gam) / np.sum((t - Phi @ mN) ** 2)
a_em, b_em, evs = 1.0, 1.0, []
for _ in range(20000):                                                          # EM: (9.63) and Exercise 9.21
    mN, S = posterior(a_em, b_em)
    evs.append(log_evidence(a_em, b_em))
    a_em = M / (mN @ mN + np.trace(S))
    b_em = Nr / (np.sum((t - Phi @ mN) ** 2) + np.trace(Phi @ S @ Phi.T))
check("EM (9.63) and the evidence re-estimation (3.92), (3.95) reach the same alpha and beta",
      abs(a_em / a_ev - 1) < 1e-6 and abs(b_em / b_ev - 1) < 1e-6,
      "alpha %.6f vs %.6f, beta %.6f vs %.6f" % (a_em, a_ev, b_em, b_ev))
check("the log evidence (3.86) never decreases along the EM iterations", np.all(np.diff(evs) > -1e-10), "from %.4f to %.4f" % (evs[0], evs[-1]))
mN, S = posterior(a_ev, b_ev)
gam = np.sum(b_ev * lam0 / (a_ev + b_ev * lam0))
check("(9.64)-(9.65) at the fixed point: alpha m^T m = gamma = M - alpha Tr(S_N)",
      abs(a_ev * mN @ mN - gam) < 1e-8 and abs(M - a_ev * np.trace(S) - gam) < 1e-8, "gamma = %.6f" % gam)

# ----------------------------------------------------------------------------- the bound of Section 9.4 (Fig. 9.14, em-bound-demo.html)
print("== The lower bound L(q, theta) (Section 9.4, Fig. 9.14)")
for name in m.EM_BOUND_SCENARIOS:
    sc = m.EM_BOUND_SCENARIOS[name]
    x = m.em_bound_data(name)
    for free in ("mu2", "pi"):
        model = sc["model"]
        grid = np.linspace(-3.5, 3.5, 141) if free == "mu2" else np.linspace(0.01, 0.99, 99)
        th0 = 0.3 if free == "mu2" else 0.3
        q = m.mix1d_posterior(x, th0, model, free)
        gap = [m.mix1d_loglik(x, th, model, free) - m.mix1d_bound(x, th, q, model, free) for th in grid]
        idn = max(abs(m.mix1d_loglik(x, th, model, free) - m.mix1d_bound(x, th, q, model, free) - m.mix1d_kl(x, th, q, model, free)) for th in grid)
        d_ll = num_grad(lambda v: m.mix1d_loglik(x, v[0], model, free), [th0], 1e-6)[0]
        d_b = num_grad(lambda v: m.mix1d_bound(x, v[0], q, model, free), [th0], 1e-6)[0]
        th1 = m.mix1d_m_step(x, q, model, free)
        d_b1 = num_grad(lambda v: m.mix1d_bound(x, v[0], q, model, free), [th1], 1e-6)[0]
        gain_ll = m.mix1d_loglik(x, th1, model, free) - m.mix1d_loglik(x, th0, model, free)
        gain_b = m.mix1d_bound(x, th1, q, model, free) - m.mix1d_bound(x, th0, q, model, free)
        ok = (min(gap) > -1e-9 and idn < 1e-9 and abs(d_ll - d_b) < 1e-5 * max(1, abs(d_ll)) and abs(d_b1) < 1e-5
              and gain_ll >= gain_b - 1e-12 and gain_b > 0)
        check("%s data, theta = %s: L <= ln p, ln p = L + KL, equal slopes at theta_old (Exercise 9.25), the M step maximises L, "
              "ln p rises at least as much as L" % (name, free), ok,
              "slopes %.5f / %.5f, gains %.4f >= %.4f" % (d_ll, d_b, gain_ll, gain_b))
x = m.em_bound_data("two")
mod = m.EM_BOUND_SCENARIOS["two"]["model"]
left, right = m.mix1d_em(x, -1.0, mod, "mu2"), m.mix1d_em(x, -0.5, mod, "mu2")
check("'two clusters' data: EM from theta = -1 and -0.5 ends at different local maxima (em-bound-demo preset 5)",
      left[-1] < -1.5 and right[-1] > 1.0,
      "%.4f (ln p %.3f) and %.4f (ln p %.3f)" % (left[-1], m.mix1d_loglik(x, left[-1], mod, "mu2"), right[-1], m.mix1d_loglik(x, right[-1], mod, "mu2")))
xo, mo = m.em_bound_data("overlap"), m.EM_BOUND_SCENARIOS["overlap"]["model"]
xs_, ms_ = m.em_bound_data("separated"), m.EM_BOUND_SCENARIOS["separated"]["model"]
slow, fast = m.mix1d_em(xo, 0.1, mo, "pi", tol=1e-6), m.mix1d_em(xs_, 0.1, ms_, "pi", tol=1e-6)
check("EM for a mixing weight is slow when the components overlap, fast when they are separated (em-bound-demo preset 4)",
      len(slow) > 5 * len(fast), "%d versus %d iterations to |change| < 1e-6" % (len(slow) - 1, len(fast) - 1))

print("\n%d checks: %d PASS, %d FAIL" % (len(RESULTS), sum(RESULTS), len(RESULTS) - sum(RESULTS)))