"""Numerical checks of the derivations in Bishop (2006), Chapter 7 (Sparse Kernel Machines).

Each check prints PASS or FAIL with the numbers behind it; the last line counts them.
Run from this folder: python3 checks_ch07.py (a few seconds). The algorithms checked are the
ones in ch07_sparse_kernel_machines.py that also produce the figures. The SVM solver is an SMO
written from scratch: its solutions are tested against the KKT conditions (7.41)-(7.42) and
(7.57)-(7.59) themselves, against the primal objective, and against a generic constrained
optimiser (scipy SLSQP) on the same dual. The RVM is tested through the evidence (7.85) computed
by two routes, monotone evidence, finite-difference derivatives and the sparsity analysis (7.96)-(7.101).
"""
import numpy as np
from scipy.optimize import minimize

import common
import ch07_sparse_kernel_machines 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 ""))


rng = np.random.default_rng(0)

# ----------------------------------------------------------------------------- 7.1 hard margin
print("== Maximum margin classifier (Section 7.1)")
Xb = rng.normal([-1.2, -1.0], 0.6, size=(15, 2))
Xr = rng.normal([1.2, 1.0], 0.6, size=(15, 2))
X = np.vstack([Xb, Xr])
t = np.r_[-np.ones(15), np.ones(15)]
h = m.svc_fit(X, t)                                        # linear kernel, hard margin
w = (h.a * t) @ X                                           # (7.29)
f = X @ w + h.b
tf = t * f
sv = h.a > 1e-8
check("hard margin: alpha_n >= 0 and sum_n alpha_n t_n = 0 (7.31)",
      h.a.min() >= 0 and abs(h.a @ t) < 1e-8, "min alpha %.1e, |sum a t| = %.1e" % (h.a.min(), abs(h.a @ t)))
check("hard margin: primal feasibility t_n y(x_n) >= 1 (7.5)", tf.min() > 1 - 1e-7, "min t y = %.9f" % tf.min())
check("hard margin: complementary slackness alpha_n (t_n y(x_n) - 1) = 0 (7.25)",
      np.abs(h.a * (tf - 1)).max() < 1e-8, "max |a (t y - 1)| = %.1e" % np.abs(h.a * (tf - 1)).max())
check("only support vectors have alpha > 0: they lie exactly on the margin, all others strictly outside",
      np.all(np.abs(tf[sv] - 1) < 1e-7) and np.all(tf[~sv] > 1 + 1e-6), "%d support vectors of %d" % (sv.sum(), len(t)))
check("dual objective (7.32) equals primal 1/2||w||^2 (7.3) at the optimum",
      abs(h.dual - 0.5 * w @ w) < 1e-8, "dual %.10f primal %.10f" % (h.dual, 0.5 * w @ w))
dist = tf / np.linalg.norm(w)
check("margin = min_n t_n y(x_n)/||w|| = 1/||w|| (7.2)-(7.4)", abs(dist.min() - 1 / np.linalg.norm(w)) < 1e-8 and abs(h.margin - 1 / np.linalg.norm(w)) < 1e-12,
      "margin %.8f, 1/||w|| %.8f" % (dist.min(), 1 / np.linalg.norm(w)))
check("kernel form ||w||^2 = sum a_n a_m t_n t_m k(x_n, x_m) equals the explicit w'w", abs(h.wnorm2 - w @ w) < 1e-9)
# a generic constrained optimiser on the same dual
Q = (t[:, None] * t[None, :]) * (X @ X.T)
res = minimize(lambda a: 0.5 * a @ Q @ a - a.sum(), np.zeros(30), jac=lambda a: Q @ a - 1, bounds=[(0, None)] * 30,
               constraints=[{"type": "eq", "fun": lambda a: a @ t, "jac": lambda a: t}], method="SLSQP", options=dict(ftol=1e-14, maxiter=500))
check("SMO dual optimum agrees with scipy SLSQP on the same dual (7.32)-(7.34)", abs(-res.fun - h.dual) < 1e-7 and np.abs(res.x - h.a).max() < 1e-4,
      "objectives %.10f vs %.10f, max |a - a_slsqp| = %.1e" % (h.dual, -res.fun, np.abs(res.x - h.a).max()))
hs = m.svc_fit(X[sv], t[sv])
Z = rng.normal(0, 1.5, size=(200, 2))
check("only the support vectors matter: refitting on them alone gives the same y(x) everywhere (7.13)",
      np.abs(hs.decision_function(Z) - h.decision_function(Z)).max() < 1e-7, "max |dy| = %.1e" % np.abs(hs.decision_function(Z) - h.decision_function(Z)).max())
Xg, tg = m.blob_ring_data()
for g in (0.8, 3.0):
    hk = m.svc_fit(Xg, tg, m.Kernel("rbf", g))
    tfk = tg * hk.f_train
    check("Gaussian kernel (gamma = %g) hard margin: KKT, margin 1/||w||, dual = primal" % g,
          abs(hk.a @ tg) < 1e-8 and tfk.min() > 1 - 1e-7 and np.abs(hk.a * (tfk - 1)).max() < 1e-8 and abs(hk.dual - 0.5 * hk.wnorm2) < 1e-8
          and abs(hk.margin - 1 / np.sqrt(hk.wnorm2)) < 1e-12, "%d SVs, margin %.4f, KKT gap tolerance met in %d iterations" % (hk.sv.sum(), hk.margin, hk.iters))

# ----------------------------------------------------------------------------- 7.1.1 soft margin
print("== Overlapping class distributions (Section 7.1.1)")
Xo, to = m.ripley_like(60, seed=2)
for kern, C in ((m.Kernel("linear"), 1.0), (m.Kernel("rbf", 0.45), 3.0), (m.Kernel("rbf", 2.0), 0.5)):
    s = m.svc_fit(Xo, to, kern, C)
    f = s.f_train
    free = (s.a > 1e-7) & (s.a < C - 1e-7)
    upper = s.a >= C - 1e-7
    zero = s.a <= 1e-7
    tf = to * f
    ok_free = np.all(np.abs(tf[free] - 1) < 1e-6)               # a < C  =>  xi = 0 and t y = 1 (7.41, 7.42, 7.45)
    ok_up = np.all(tf[upper] <= 1 + 1e-6)                         # a = C  =>  xi >= 0, t y = 1 - xi
    ok_zero = np.all(tf[zero] >= 1 - 1e-6)
    cs = np.abs(s.a * (tf - 1 + s.xi)).max()
    name = "%s kernel, C = %g" % (kern.kind + ("" if kern.kind == "linear" else " gamma=%g" % kern.gamma), C)
    check("soft margin KKT (%s): 0<a<C on the margin, a = C inside or beyond it, a = 0 outside; a_n(t y - 1 + xi) = 0; sum a t = 0; 0 <= a <= C (7.43)-(7.45)" % name,
          ok_free and ok_up and ok_zero and cs < 1e-7 and abs(s.a @ to) < 1e-8 and s.a.min() >= 0 and s.a.max() <= C + 1e-12,
          "free %d, at bound %d, zero %d; max |a(ty-1+xi)| = %.1e" % (free.sum(), upper.sum(), zero.sum(), cs))
    check("soft margin (%s): dual objective (7.43) = primal 1/2||w||^2 + C sum xi (7.40)" % name, abs(s.dual - s.primal) < 1e-6 * max(1, abs(s.primal)),
          "dual %.9f primal %.9f gap %.1e" % (s.dual, s.primal, abs(s.dual - s.primal)))
    if kern.kind == "linear":
        w_ = (s.a * to) @ Xo
        pr = 0.5 * w_ @ w_ + C * np.maximum(0, 1 - to * (Xo @ w_ + s.b)).sum()
        check("... with w, b taken explicitly (7.29), (7.37): primal objective evaluated independently", abs(pr - s.dual) < 1e-6, "primal %.9f dual %.9f" % (pr, s.dual))
Xs, ts = Xo[::3], to[::3]
s = m.svc_fit(Xs, ts, m.Kernel("rbf", 0.5), 2.0)
Q = (ts[:, None] * ts[None, :]) * s.K
res = minimize(lambda a: 0.5 * a @ Q @ a - a.sum(), np.zeros(len(ts)), jac=lambda a: Q @ a - 1, bounds=[(0, 2.0)] * len(ts),
               constraints=[{"type": "eq", "fun": lambda a: a @ ts, "jac": lambda a: ts}], method="SLSQP", options=dict(ftol=1e-14, maxiter=500))
check("soft-margin SMO agrees with scipy SLSQP on the box-constrained dual", abs(-res.fun - s.dual) < 1e-7, "objectives %.10f vs %.10f" % (s.dual, -res.fun))
nu = m.nu_svc_fit(Xo, to, 0.6, m.Kernel("rbf", 0.45))
check("nu-SVM properties (7.50)-(7.52): fraction of margin errors <= nu <= fraction of support vectors",
      (nu.xi > 1e-7).mean() <= 0.6 + 1e-9 and nu.sv.mean() >= 0.6 - 0.5 / len(to) - 1e-9,
      "margin errors %.3f, SVs %.3f, nu = 0.6 (C = %.3g)" % ((nu.xi > 1e-7).mean(), nu.sv.mean(), nu.C))

# ----------------------------------------------------------------------------- 7.1.4 regression
print("== SVM for regression (Section 7.1.4)")
Xr_, tr_ = m.sinusoid_data(25, 0.2, seed=4)
k = m.Kernel("rbf", 15.0)
for C, eps in ((5.0, 0.15), (50.0, 0.05), (1.0, 0.3)):
    r = m.svr_fit(Xr_, tr_, k, C, eps)
    f = r.f_train
    z = tr_ - f
    a, ah = r.a_up, r.a_lo
    tube = np.abs(z) < eps - 1e-6
    out = np.abs(z) > eps + 1e-6
    k1 = np.abs(a * (eps + r.xi + f - tr_)).max()     # a_n (eps + xi_n + y_n - t_n) = 0 (7.57)
    k2 = np.abs(ah * (eps + r.xi_hat - f + tr_)).max()                                           # a^_n (eps + xi^_n - y_n + t_n) = 0 (7.58)
    ok = (abs(coef_sum := r.coef.sum()) < 1e-8 and np.abs(a * ah).max() < 1e-10 and a.min() >= 0 and ah.min() >= 0 and a.max() <= C + 1e-12 and ah.max() <= C + 1e-12
          and k1 < 1e-7 and k2 < 1e-7)
    check("epsilon-SVR KKT (C = %g, eps = %g): sum (a - a^) = 0 (7.60), a a^ = 0, 0 <= a, a^ <= C, (7.57)-(7.58) complementary slackness" % (C, eps), ok,
          "|sum coef| %.1e, max slack %.1e / %.1e" % (abs(coef_sum), k1, k2))
    check("epsilon-SVR tube (C = %g, eps = %g): points strictly inside have a = a^ = 0; points outside have a or a^ equal to C (7.65-style)" % (C, eps),
          np.all((a[tube] < 1e-8) & (ah[tube] < 1e-8)) and np.all(np.maximum(a, ah)[out] > C - 1e-6), "%d inside, %d outside, %d SVs" % (tube.sum(), out.sum(), r.sv.sum()))
    check("epsilon-SVR (C = %g, eps = %g): dual (7.61) = primal (7.55)" % (C, eps), abs(r.dual - r.primal) < 1e-6 * max(1, abs(r.primal)), "dual %.9f primal %.9f" % (r.dual, r.primal))
r = m.svr_fit(Xr_, tr_, k, 5.0, 0.15)
r2 = m.svr_fit(Xr_[r.sv], tr_[r.sv], k, 5.0, 0.15)
xs = np.linspace(0, 1, 100)[:, None]
check("SVR predictions use only the support vectors: refitting on them gives the same curve (7.64)", np.abs(r.predict(xs) - r2.predict(xs)).max() < 1e-6,
      "max |dy| = %.1e" % np.abs(r.predict(xs) - r2.predict(xs)).max())
nv = m.nu_svr_fit(Xr_, tr_, k, 5.0, 0.3)
outside = ((nv.xi + nv.xi_hat) > 1e-7).sum()
check("nu-SVR (7.68)-(7.72): at most nu N points outside the tube, at least nu N support vectors", outside <= 0.3 * 25 + 1e-9 and nv.sv.sum() >= 0.3 * 25 - 0.5,
      "outside %d, SVs %d, nu N = %.1f, epsilon = %.4f" % (outside, nv.sv.sum(), 7.5, nv.eps))

# ----------------------------------------------------------------------------- 7.2 RVM regression
print("== RVM for regression (Section 7.2.1)")
Xs_, ts_ = m.sinusoid_data(10, 0.2, seed=3)
kk = m.Kernel("rbf", 15.0)
Phi = m.design_matrix(Xs_, kk)
alpha = rng.uniform(0.5, 4, Phi.shape[1])
beta = 20.0
d = m.log_evidence_direct(Phi, ts_, alpha, beta)
wb = m.log_evidence_woodbury(Phi, ts_, alpha, beta)
check("marginal likelihood by two routes: direct N x N form (7.85)-(7.86) = M x M posterior form (3.86)", abs(d - wb) < 1e-9, "%.10f vs %.10f" % (d, wb))
# posterior mean is the mode of ln p(w|t)
mu, S = m.rvm_posterior(Phi, ts_, alpha, beta)
grad = -beta * Phi.T @ (Phi @ mu - ts_) - alpha * mu
check("posterior mean (7.82) is the stationary point of the log posterior, covariance (7.83) the inverse Hessian", np.abs(grad).max() < 1e-9, "max |grad| = %.1e" % np.abs(grad).max())
# d ln p / d alpha_i = 1/2 (1/alpha_i - m_i^2 - Sigma_ii): the stationarity behind (7.87)
fd = np.zeros(len(alpha))
for i in range(len(alpha)):
    e = np.zeros(len(alpha))
    e[i] = 1e-6
    fd[i] = (m.log_evidence_direct(Phi, ts_, alpha + e, beta) - m.log_evidence_direct(Phi, ts_, alpha - e, beta)) / 2e-6
an = 0.5 * (1 / alpha - mu ** 2 - np.diag(S))
check("d ln p(t|alpha, beta)/d alpha_i = (1/alpha_i - m_i^2 - Sigma_ii)/2 by finite differences: the source of (7.87)", np.abs(fd - an).max() < 1e-6, "max |difference| %.1e" % np.abs(fd - an).max())
g = 1 - alpha * np.diag(S)
fdb = (m.log_evidence_direct(Phi, ts_, alpha, beta + 1e-5) - m.log_evidence_direct(Phi, ts_, alpha, beta - 1e-5)) / 2e-5
anb = 0.5 * (len(ts_) / beta - np.sum((ts_ - Phi @ mu) ** 2) - np.trace(S @ Phi.T @ Phi))
check("d ln p/d beta = (N/beta - ||t - Phi m||^2 - tr(Sigma Phi'Phi))/2 by finite differences: the source of (7.88)", abs(fdb - anb) < 1e-6, "%.8f vs %.8f" % (fdb, anb))
em = m.rvm_regression(Phi, ts_, update="em", n_iter=300, prune=1e12)
dtr = np.diff(em["trace"])
check("EM updates of alpha and beta: ln p(t|alpha, beta) never decreases", dtr.min() > -1e-9, "%d iterations, min step %.1e, total gain %.3f" % (len(dtr), dtr.min(), em["trace"][-1] - em["trace"][0]))
mk = m.rvm_regression(Phi, ts_, update="mackay")
dtm = np.diff(mk["trace"])
check("re-estimation (7.87)-(7.89): ln p(t|alpha, beta) increases at every update (before pruning removes columns)", dtm[:5].min() > 0 and mk["trace"][-1] > mk["trace"][0],
      "%d iterations, %d of %d steps decrease (largest decrease %.1e, caused by pruning columns), total gain %.3f" % (len(dtm), (dtm < -1e-9).sum(), len(dtm), dtm.min(), mk["trace"][-1] - mk["trace"][0]))
check("RVM: few relevance vectors (sparse)", len(mk["active"]) <= 4, "%d of %d basis functions (including the bias) survive; the nu-SVR needs 7 of 10" % (len(mk["active"]), Phi.shape[1]))
# fixed point of (7.87): d ln p/d alpha = 0 on the survivors
P = Phi[:, mk["active"]]
mu_k, S_k = m.rvm_posterior(P, ts_, mk["alpha"], mk["beta"])
check("at convergence alpha_i = gamma_i / m_i^2 (7.87) holds on the surviving basis functions",
      np.abs(mk["alpha"] - (1 - mk["alpha"] * np.diag(S_k)) / mu_k ** 2).max() < 1e-4 * mk["alpha"].max(), "max relative gap %.1e" % np.abs(mk["alpha"] / ((1 - mk["alpha"] * np.diag(S_k)) / mu_k ** 2) - 1).max())

print("== Sparsity analysis (Section 7.2.2)")
alpha = rng.uniform(0.5, 4, Phi.shape[1])
beta = 15.0
ok_all = True
worst = 0
for i in (1, 3, 5, 8):
    others = np.arange(Phi.shape[1]) != i
    Cm = np.eye(10) / beta + (Phi[:, others] / alpha[others]) @ Phi[:, others].T
    ph = Phi[:, i]
    s_i = ph @ np.linalg.solve(Cm, ph)                                # (7.98)
    q_i = ph @ np.linalg.solve(Cm, ts_)                               # (7.99)
    lhs = m.log_evidence_direct(Phi, ts_, alpha, beta) - m.log_evidence_direct(Phi[:, others], ts_, alpha[others], beta)
    rhs = m.lam(alpha[i], q_i ** 2, s_i)                              # (7.96)-(7.97)
    C_full = np.eye(10) / beta + (Phi / alpha) @ Phi.T
    S_i, Q_i = ph @ np.linalg.solve(C_full, ph), ph @ np.linalg.solve(C_full, ts_)
    s2 = alpha[i] * S_i / (alpha[i] - S_i)                            # (7.105)
    q2 = alpha[i] * Q_i / (alpha[i] - S_i)                            # (7.104)
    worst = max(worst, abs(lhs - rhs), abs(s2 - s_i), abs(q2 - q_i))
check("L(alpha) = L(alpha_-i) + lambda(alpha_i) with s_i, q_i from (7.98)-(7.99) (7.96)-(7.97); and (7.104)-(7.105) give the same s_i, q_i", worst < 1e-8, "worst difference %.1e" % worst)
al = np.exp(np.linspace(-6, 6, 200001))
for q2, s in ((4.0, 1.0), (1.0, 2.0), (9.0, 2.0)):
    lm = m.lam(al, q2, s)
    am = al[np.argmax(lm)]
    dl = (m.lam(al + 1e-6 * al, q2, s) - m.lam(al - 1e-6 * al, q2, s)) / (2e-6 * al)
    d_formula = (al ** -1 * s ** 2 - (q2 - s)) / (2 * (al + s) ** 2)          # (7.100)
    if q2 > s:
        ok = abs(am - s ** 2 / (q2 - s)) / am < 1e-3
    else:
        ok = np.argmax(lm) == len(al) - 1 and np.all(np.diff(lm) > 0)
    check("lambda(alpha), q^2 = %g, s = %g: %s (7.101); derivative (7.100) matches finite differences" % (q2, s, "maximum at s^2/(q^2-s)" if q2 > s else "increasing, maximum at infinity"),
          ok and np.abs(dl - d_formula).max() < 1e-5, "argmax %.4g, expected %s" % (am, "%.4g" % (s ** 2 / (q2 - s)) if q2 > s else "inf"))
beta = 1 / 0.04
fast = m.rvm_fast_regression(Phi, ts_, beta)
df = np.diff(fast["trace"])
check("sequential sparse Bayesian learning: every add, re-estimate or delete leaves ln p(t|alpha) non-decreasing", len(df) > 3 and df.min() > -1e-9,
      "%d actions, min change %.1e, final ln p %.4f with %d basis functions" % (len(df), df.min(), fast["trace"][-1], len(fast["active"])))
act = fast["active"]
al_f = np.full(Phi.shape[1], np.inf)
al_f[act] = fast["alpha"]
cur = m.log_evidence_direct(Phi[:, act], ts_, fast["alpha"], beta)
best_gain = -np.inf
for i in range(Phi.shape[1]):
    for cand in np.exp(np.linspace(-8, 12, 161)):
        a2 = al_f.copy()
        a2[i] = cand
        ac2 = np.isfinite(a2)
        best_gain = max(best_gain, m.log_evidence_direct(Phi[:, ac2], ts_, a2[ac2], beta) - cur)
    a2 = al_f.copy()
    a2[i] = np.inf
    if np.isfinite(a2).any():
        ac2 = np.isfinite(a2)
        best_gain = max(best_gain, m.log_evidence_direct(Phi[:, ac2], ts_, a2[ac2], beta) - cur)
check("the fast algorithm stops at a coordinate-wise maximum: no change of a single alpha_i (any value, or infinity) raises the evidence, and few basis functions remain",
      len(act) <= 4 and best_gain < 1e-3, "best single-coordinate gain on a grid %.1e (grid resolution 0.125 in ln alpha), %d basis functions" % (best_gain, len(act)))
mu_p, var_p = m.rvm_predict(m.design_matrix(Xs_, kk, centres=Xs_), mk)
check("predictive variance (7.91) = 1/beta + phi' Sigma phi >= 1/beta", np.all(var_p >= 1 / mk["beta"] - 1e-12), "min %.4f, 1/beta %.4f" % (var_p.min(), 1 / mk["beta"]))

# ----------------------------------------------------------------------------- 7.2.3 classification
print("== RVM for classification (Section 7.2.3)")
Xc, tc_pm = m.ripley_like(60, seed=5)
tc = (tc_pm > 0).astype(float)
Phic = m.design_matrix(Xc, m.Kernel("rbf", 0.45))
alc = rng.uniform(0.2, 2, Phic.shape[1])
w0 = rng.normal(0, 0.3, Phic.shape[1])
g = m.logpost_grad(w0, Phic, tc, alc)
fdg = np.zeros_like(w0)
for i in range(len(w0)):
    e = np.zeros_like(w0)
    e[i] = 1e-6
    fdg[i] = (m._bern_logpost(w0 + e, Phic, tc, alc) - m._bern_logpost(w0 - e, Phic, tc, alc)) / 2e-6
check("gradient (7.110) of the log posterior (7.109) agrees with finite differences", np.abs(g - fdg).max() < 1e-6, "max |difference| %.1e" % np.abs(g - fdg).max())
Hn = np.zeros((len(w0), len(w0)))
for i in range(len(w0)):
    e = np.zeros_like(w0)
    e[i] = 1e-6
    Hn[:, i] = (m.logpost_grad(w0 + e, Phic, tc, alc) - m.logpost_grad(w0 - e, Phic, tc, alc)) / 2e-6
check("Hessian (7.111) = -(Phi'B Phi + A) agrees with finite differences of the gradient", np.abs(Hn - m.logpost_hess(w0, Phic, tc, alc)).max() < 1e-6,
      "max |difference| %.1e" % np.abs(Hn - m.logpost_hess(w0, Phic, tc, alc)).max())
ws, Sg = m.laplace_mode(Phic, tc, alc)
yy = common.sigmoid(Phic @ ws)
check("IRLS mode satisfies w* = A^-1 Phi'(t - y) (7.112)", np.abs(ws - (Phic.T @ (tc - yy)) / alc).max() < 1e-8, "max |difference| %.1e" % np.abs(ws - (Phic.T @ (tc - yy)) / alc).max())
check("Laplace covariance (7.113) = (Phi'B Phi + A)^-1", np.abs(Sg - np.linalg.inv(Phic.T @ (Phic * (yy * (1 - yy))[:, None]) + np.diag(alc))).max() < 1e-10)
rv = m.rvm_classification(Phic, tc)
dtr = np.diff(rv["trace"])
check("RVM classification: Laplace evidence (7.114) rises overall under the updates (7.116) and the solution is sparse",
      rv["trace"][-1] > rv["trace"][0] and len(rv["active"]) <= 0.2 * Phic.shape[1],
      "evidence %.3f -> %.3f, %d of %d basis functions survive, %d outer iterations" % (rv["trace"][0], rv["trace"][-1], len(rv["active"]), Phic.shape[1], rv["iters"]))
sv_ = m.nu_svc_fit(Xc, tc_pm, 0.6, m.Kernel("rbf", 0.45))
check("the RVM is sparser than the nu-SVM on the same data", len(rv["active"]) < sv_.sv.sum() / 3, "RVM %d vs SVM %d basis functions" % (len(rv["active"]), sv_.sv.sum()))
p_ = m.rvm_class_prob(Phic, rv, False)
check("RVM class probabilities lie in (0, 1) and classify better than chance", np.all((p_ > 0) & (p_ < 1)) and np.mean((p_ > 0.5) == (tc > 0.5)) > 0.7,
      "training accuracy %.3f" % np.mean((p_ > 0.5) == (tc > 0.5)))

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