"""Numerical checks of the derivations in Bishop (2006), Chapter 5 (Neural Networks).

Each check prints PASS or FAIL with the numbers behind it; the last line counts them.
Run from this folder: python3 checks_ch05.py (well under a minute). The algorithms checked are the
ones in ch05_neural_networks.py, which also produce the figures: backpropagation, the Jacobian and
the Hessian (exact, by the R-operator, outer product, diagonal) are compared with finite differences
and with each other, early stopping with weight decay in a quadratic error, the mixture density
network derivatives with finite differences, and the Laplace evidence with its closed form.
"""
import math

import numpy as np
from scipy import optimize, stats

import common
import ch05_neural_networks 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 kinds_data(kind, N, K, rng):
    X = rng.normal(size=(N, 2))
    if kind == "sse":
        return X, rng.normal(size=(N, K))
    if kind == "bce":
        return X, (rng.random((N, 1)) < 0.5).astype(float)
    return X, np.eye(K)[rng.integers(0, K, N)]


rng = np.random.default_rng(0)
CASES = [("sse", 2), ("bce", 1), ("ce", 3)]

print("== Section 5.3: backpropagation, the canonical links and the Jacobian")
for kind, K in CASES:
    net = m.Net([2, 4, K], kind)
    w = rng.normal(size=net.W) * 0.7
    X, T = kinds_data(kind, 9, K, rng)
    g, gf = net.grad(w, X, T), net.grad_fd(w, X, T)
    check("backprop gradient = finite differences (5.69), two layers, %s" % kind, np.abs(g - gf).max() < 1e-7, "max |diff| = %.1e" % np.abs(g - gf).max())
net3 = m.Net([2, 5, 4, 1], "bce")
w3 = rng.normal(size=net3.W) * 0.5
X3, T3 = kinds_data("bce", 8, 1, rng)
d3 = np.abs(net3.grad(w3, X3, T3) - net3.grad_fd(w3, X3, T3)).max()
check("backprop gradient = finite differences, three layers (two hidden)", d3 < 1e-7, "max |diff| = %.1e" % d3)
for kind, K in CASES:
    net = m.Net([2, 4, K], kind)
    w = rng.normal(size=net.W) * 0.7
    X, T = kinds_data(kind, 6, K, rng)
    c = net.forward(w, X)
    dE = np.zeros_like(c["a"])
    for n in range(len(X)):
        for k in range(K):
            e = np.zeros_like(c["a"])
            e[n, k] = 1e-6
            dE[n, k] = (net.error_from_a(c["a"] + e, T) - net.error_from_a(c["a"] - e, T)) / 2e-6
    ok = np.abs(dE - (net.act(c["a"]) - T)).max() < 1e-7
    check("canonical link: dE_n/da_k = y_k - t_k (5.18, 5.66, 5.69), %s" % kind, ok, "max |diff| = %.1e" % np.abs(dE - (net.act(c["a"]) - T)).max())
for kind, K in CASES:
    net = m.Net([3, 4, K], kind)
    w = rng.normal(size=net.W) * 0.7
    X = rng.normal(size=(5, 3))
    J, Jf = net.jacobian(w, X), net.jacobian_fd(w, X)
    check("Jacobian by backpropagation (5.73)-(5.76) = finite differences (5.77), %s outputs" % kind, np.abs(J - Jf).max() < 1e-7, "max |diff| = %.1e" % np.abs(J - Jf).max())
Jd = m.Net([2, 5, 4, 2], "sse")
wd_ = rng.normal(size=Jd.W) * 0.5
Xd = rng.normal(size=(4, 2))
check("Jacobian, three layers", np.abs(Jd.jacobian(wd_, Xd) - Jd.jacobian_fd(wd_, Xd)).max() < 1e-7)
net = m.Net([1, 3, 1])
w = rng.normal(size=net.W)
w2 = w.copy()
W1, W2 = net.unpack(w2)
W1[1, :] *= -1
W2[0, 2] *= -1
xs = np.linspace(-2, 2, 20)[:, None]
check("weight-space symmetry: flipping the signs into and out of a tanh unit leaves the function unchanged (Section 5.1.1)", np.abs(net.output(w, xs) - net.output(net.pack([W1, W2]), xs)).max() < 1e-12)
Wp = net.unpack(w)
perm = [2, 0, 1]
W1p, W2p = Wp[0][perm], Wp[1].copy()
W2p[:, 1:] = W2p[:, 1:][:, perm]
check("weight-space symmetry: permuting the hidden units leaves the function unchanged (M! symmetries)", np.abs(net.output(w, xs) - net.output(net.pack([W1p, W2p]), xs)).max() < 1e-12)

print("== Section 5.4: the Hessian")
for kind, K in CASES:
    net = m.Net([2, 4, K], kind)
    w = rng.normal(size=net.W) * 0.7
    X, T = kinds_data(kind, 9, K, rng)
    H = net.hessian_exact(w, X, T)
    Hf = net.hessian_fd(w, X, T)
    Hr = net.hessian_rop(w, X, T)
    check("exact Hessian (5.93)-(5.95) = finite differences of the gradient (5.91), %s" % kind, np.abs(H - Hf).max() < 1e-6, "max |diff| = %.1e" % np.abs(H - Hf).max())
    check("exact Hessian = W R-operator products H e_i (5.96)-(5.108), %s" % kind, np.abs(H - Hr).max() < 1e-10, "max |diff| = %.1e" % np.abs(H - Hr).max())
    v = rng.normal(size=net.W)
    check("R-operator product H v = (exact H) v, %s" % kind, np.abs(net.hessian_vector(w, X, T, v) - H @ v).max() < 1e-9, "max |diff| = %.1e" % np.abs(net.hessian_vector(w, X, T, v) - H @ v).max())
    check("Hessian symmetric", np.abs(H - H.T).max() < 1e-12)
Hd = net3.hessian_rop(w3, X3, T3)
check("R-operator Hessian = finite differences of the gradient, three layers", np.abs(Hd - net3.hessian_fd(w3, X3, T3)).max() < 1e-6, "max |diff| = %.1e" % np.abs(Hd - net3.hessian_fd(w3, X3, T3)).max())
v = rng.normal(size=net3.W)
eps = 1e-6
fd = (net3.grad(w3 + eps * v, X3, T3) - net3.grad(w3 - eps * v, X3, T3)) / (2 * eps)
check("R{grad E} = v^T grad grad E: directional derivative of the gradient (5.96)-(5.97), three layers", np.abs(net3.hessian_vector(w3, X3, T3, v) - fd).max() < 1e-6)

# outer product approximation
for kind, K in CASES:
    net = m.Net([2, 4, K], kind)
    w = rng.normal(size=net.W) * 0.7
    X = rng.normal(size=(9, 2))
    T = net.output(w, X)                               # zero residuals y_n = t_n
    Ho, He = net.hessian_outer(w, X, T), net.hessian_exact(w, X, T)
    check("outer product approximation = exact Hessian when y_n = t_n (5.83)-(5.85), %s" % kind, np.abs(Ho - He).max() < 1e-10, "max |diff| = %.1e" % np.abs(Ho - He).max())
lin = m.Net([3, 2], "sse")                              # a linear model: no hidden layer
wl = rng.normal(size=lin.W)
Xl = rng.normal(size=(11, 3))
Tl = rng.normal(size=(11, 2))
Hl = lin.hessian_rop(wl, Xl, Tl)
check("linear model: outer product approximation equals the Hessian for ANY targets (grad grad y = 0 in (5.83))", np.abs(lin.hessian_outer(wl, Xl) - Hl).max() < 1e-10, "max |diff| = %.1e" % np.abs(lin.hessian_outer(wl, Xl) - Hl).max())
Xt = np.hstack([np.ones((11, 1)), Xl])
check("linear model: H = Phi^T Phi per output (the least-squares Hessian)", np.abs(Hl[:4, :4] - Xt.T @ Xt).max() < 1e-10)
net = m.Net([2, 4, 1], "sse")
w = rng.normal(size=net.W) * 0.7
X, T = kinds_data("sse", 40, 1, rng)
Ho, He = net.hessian_outer(w, X, T), net.hessian_exact(w, X, T)
check("outer product approximation differs from the Hessian for a poorly fitted network (second term of (5.83) not negligible)", np.abs(Ho - He).max() > 1e-2, "max |diff| = %.2f" % np.abs(Ho - He).max())
# diagonal
net = m.Net([2, 4, 1], "sse")
w = rng.normal(size=net.W) * 0.7
X, T = kinds_data("sse", 12, 1, rng)
check("diagonal elements (5.79)-(5.80) = diagonal of the exact Hessian", np.abs(net.hessian_diagonal(w, X, T) - np.diag(net.hessian_exact(w, X, T))).max() < 1e-10)
check("Becker-Le Cun diagonal (5.81) = exact diagonal for a single output with M = I", np.abs(net.hessian_diagonal(w, X, T, True) - np.diag(net.hessian_exact(w, X, T))).max() < 1e-10)
net = m.Net([2, 4, 3], "ce")
w = rng.normal(size=net.W) * 0.7
X, T = kinds_data("ce", 12, 3, rng)
dd = np.abs(net.hessian_diagonal(w, X, T, True) - np.diag(net.hessian_exact(w, X, T))).max()
check("Becker-Le Cun (5.81) drops the off-diagonal M_kk' and so differs from the exact diagonal for softmax outputs", dd > 1e-4, "max |diff| = %.3f" % dd)
# inverse Hessian
net = m.Net([2, 4, 1], "sse")
w = rng.normal(size=net.W) * 0.7
X = rng.normal(size=(15, 2))
B = net.grad_activation(w, X)[:, 0, :]
for a0 in (1e-4, 1e-1):
    Pinv = m.inverse_hessian_sequential(B, a0)
    ref = np.linalg.inv(B.T @ B + a0 * np.eye(net.W))
    check("sequential inverse Hessian (5.87)-(5.89) = inverse of (sum b b^T + alpha I), alpha = %g" % a0, np.abs(Pinv - ref).max() < 1e-6 * max(1, np.abs(ref).max()), "max |diff| = %.1e" % np.abs(Pinv - ref).max())

print("== Section 5.5: weight decay, early stopping, invariances")
lams = np.array([0.1, 0.5, 2.0, 10.0])
w_ml = np.array([1.0, -2.0, 0.5, 1.5])
eta, tau = 0.01, 300
w = np.zeros(4)
for _ in range(tau):                                      # gradient descent on E = 1/2 sum lam (w - w_ML)^2
    w = w - eta * lams * (w - w_ml)
es, wdk = m.early_stopping_vs_weight_decay(lams, w_ml, eta, tau, 1.0)
check("gradient descent from 0 on a quadratic: w_i = (1 - (1 - eta lam_i)^tau) w_ML,i (Section 5.5.2)", np.abs(w - es).max() < 1e-12)
check("weight decay solution (H + alpha I)^-1 H w_ML = lam/(lam + alpha) w_ML (3.28)", np.abs(wdk - np.linalg.solve(np.diag(lams) + np.eye(4), np.diag(lams) @ w_ml)).max() < 1e-12)
s = np.linspace(0, 40, 400001)
gap = np.abs(1 - np.exp(-s) - s / (1 + s)).max()
check("early stopping ~ weight decay with alpha = 1/(eta tau): the shrinkage factors 1 - exp(-s) and s/(1 + s) (s = eta tau lam) differ by at most %.3f" % gap, abs(gap - 0.2036) < 1e-3)
lam_grid = np.logspace(-3, 3, 50)
a1, a2 = m.early_stopping_vs_weight_decay(lam_grid, np.ones(50), 1e-5, 20000, 1.0 / (1e-5 * 20000))
check("early stopping (eta -> 0 at eta tau = 0.2) is within 0.21 of weight decay with alpha = 1/(eta tau) for every eigenvalue (the same maximum)", np.abs(a1 - a2).max() < 0.21, "max |diff| = %.3f" % np.abs(a1 - a2).max())
check("both shrink small-curvature directions toward 0 and keep large-curvature ones (eigenvalue >> alpha)", abs(a1[-1] - 1) < 1e-6 and abs(a2[-1] - 1) < 1e-2 and a1[0] < 1e-2 and a2[0] < 1e-2)
# weight decay in a network: (5.112) gradient and the equivalent MAP
net = m.Net([1, 3, 1])
w = rng.normal(size=net.W)
X, T = rng.normal(size=(7, 1)), rng.normal(size=(7, 1))
lam = 0.3
gfun = lambda w_: net.error(w_, X, T) + 0.5 * lam * w_ @ w_
gfd = np.array([(gfun(w + 1e-6 * e) - gfun(w - 1e-6 * e)) / 2e-6 for e in np.eye(net.W)])
check("weight decay (5.112): gradient of E + lambda/2 w^T w is grad E + lambda w", np.abs(gfd - (net.grad(w, X, T) + lam * w)).max() < 1e-6)
# tangent vector
rngd = np.random.default_rng(3)
Xm, ym = common.mnist(train=False, digits=(6,))
img = Xm[0].reshape(28, 28)
tau_v = m.rotate_tangent(img)
errs = []
for deg in (1.0, 2.0, 4.0):
    sm = m.ndimage.gaussian_filter(img, 0.6)
    lin_img = sm + math.radians(deg) * tau_v
    errs.append(np.abs(lin_img - m.rotate_image(sm, deg)).max())
check("tangent vector: x + epsilon tau approximates the rotated image with error O(epsilon^2) (5.126)", errs[1] > 3 * errs[0] and errs[2] > 3 * errs[1], "errors %.4f %.4f %.4f for 1, 2, 4 degrees" % tuple(errs))

print("== Section 5.6: mixture density networks")
K = 3
net = m.Net([1, 4, 3 * K], "mdn", K_mix=K)
w = rng.normal(size=net.W) * 0.6
x = rng.random((10, 1))
t = rng.random((10, 1))
c = net.forward(w, x)
pi, sg, mu = m.mdn_parameters(c["a"], K)
check("mixing coefficients are positive and sum to 1 (softmax (5.150))", np.abs(pi.sum(1) - 1).max() < 1e-12 and pi.min() > 0, "max |sum - 1| = %.1e" % np.abs(pi.sum(1) - 1).max())
tg = np.linspace(-30, 30, 200001)
dens = m.mdn_density(c["a"][:3], tg, K)
check("p(t|x) of (5.148) integrates to 1", np.abs(np.trapz(dens, tg, axis=1) - 1).max() < 1e-6)
dE = m.mdn_delta(c["a"], t[:, 0], K)
dEf = np.zeros_like(dE)
for n in range(10):
    for j in range(3 * K):
        e = np.zeros_like(c["a"])
        e[n, j] = 1e-6
        dEf[n, j] = (m.mdn_error(c["a"] + e, t[:, 0], K) - m.mdn_error(c["a"] - e, t[:, 0], K)) / 2e-6
check("MDN output derivatives (5.155), (5.156), (5.157 with the erratum) = finite differences", np.abs(dE - dEf).max() < 1e-7, "max |diff| = %.1e" % np.abs(dE - dEf).max())
check("MDN weight gradient by backpropagation = finite differences", np.abs(net.grad(w, x, t) - net.grad_fd(w, x, t)).max() < 1e-7, "max |diff| = %.1e" % np.abs(net.grad(w, x, t) - net.grad_fd(w, x, t)).max())
lp = m.mdn_log_components(c["a"], t[:, 0], K)
gam = np.exp(lp - common.logsumexp(lp, axis=1)[:, None])
check("posterior responsibilities gamma_k (5.154) sum to 1 over k", np.abs(gam.sum(1) - 1).max() < 1e-12)
sg_check = np.abs(m.mdn_error(c["a"], t[:, 0], K) - (-np.sum(np.log((pi * stats.norm.pdf(t, mu, sg)).sum(1))))) < 1e-9
check("MDN error (5.153) = - sum ln sum_k pi_k N(t | mu_k, sigma_k^2)", sg_check)
x1, t1 = m.inverse_problem_data(80)
netm = m.Net([1, 4, 3], "mdn", K_mix=1)
wm = m.init_weights(netm, m.JSRandom(1), 1.0, 0.5)
wm2 = m.adam_train(netm, wm, x1[:, None], t1[:, None], lr=0.02, steps=2000)
a2 = netm.forward(wm2, x1[:, None])["a"]
_, s1, mu1 = m.mdn_parameters(a2, 1)
res = t1 - mu1[:, 0]
check("K = 1 mixture density network = Gaussian regression with input-dependent variance: error falls below the least-squares value", netm.error(wm2, x1[:, None], t1[:, None]) < netm.error(wm, x1[:, None], t1[:, None]) and np.all(s1 > 0))

print("== Section 5.7: Bayesian neural networks")
# exact evidence for a linear-Gaussian model equals the Laplace formula (5.175), which is exact when H is constant
N = 15
xx = rng.random(N)
tt = np.sin(2 * np.pi * xx) + 0.3 * rng.standard_normal(N)
Phi = np.c_[np.ones(N), xx, xx ** 2]
alpha, beta = 2.0, 8.0
cov = np.eye(N) / beta + Phi @ Phi.T / alpha
exact = stats.multivariate_normal(np.zeros(N), cov).logpdf(tt)
linnet = m.Net([2, 1], "sse")
Xl2 = np.c_[xx, xx ** 2]
H = beta * 0 + linnet.hessian_rop(np.zeros(linnet.W), Xl2, tt[:, None])
A = beta * H + alpha * np.eye(3)
wmap = np.linalg.solve(A, beta * Phi.T @ tt)
Emap = beta * 0.5 * np.sum((Phi @ wmap - tt) ** 2) + 0.5 * alpha * wmap @ wmap
lap = m.laplace_log_evidence(Emap, H, alpha, beta, N, 3)
check("Laplace evidence (5.175) = exact marginal likelihood (3.86) for a model linear in its weights", abs(lap - exact) < 1e-9, "%.6f vs %.6f" % (lap, exact))
# gamma and the alpha, beta fixed points
net = m.Net([1, 3, 1])
xs_ = rng.random(20)[:, None]
ts_ = (np.sin(2 * np.pi * xs_) + 0.25 * rng.standard_normal((20, 1)))
ev = m.evidence_regression(net, xs_, ts_, m.init_weights(net, m.JSRandom(3), 1.0, 0.5), 1.0, 10.0, outer=40, history=True)
w, al, be, gam = ev["w"], ev["alpha"], ev["beta"], ev["gamma"]
ED = net.error(w, xs_, ts_)
check("evidence procedure converges: alpha = gamma / w^T w (5.178)", abs(al - gam / (w @ w)) < 1e-4 * al, "alpha = %.4f, gamma/w^Tw = %.4f" % (al, gam / (w @ w)))
check("evidence procedure converges: 1/beta = 2 E_D / (N - gamma) (5.180-5.182)", abs(1 / be - 2 * ED / (len(xs_) - gam)) < 1e-4 / be, "1/beta = %.5f" % (1 / be))
lam = ev["eig"]
check("gamma = sum lam_i/(alpha + lam_i) (5.179) lies in (0, W) and alpha + lam_i > 0 (the regularised error has a minimum there; beta H itself may have eigenvalues down to -alpha)", 0 < gam < net.W and lam.min() > -al, "gamma = %.2f of W = %d, min eig %.1e" % (gam, net.W, lam.min()))
# the identities behind (5.178) and (5.180), for a fixed w and a positive semi-definite H (the outer product approximation)
Hp = net.hessian_outer(w, xs_)
lam_p = np.linalg.eigvalsh(Hp)
f_alpha = lambda a: -m.laplace_log_evidence(be * ED + 0.5 * a * w @ w, Hp, a, be, len(xs_), net.W)
opt = optimize.minimize_scalar(lambda la: f_alpha(math.exp(la)), bounds=(-8, 8), method="bounded", options=dict(xatol=1e-10))
a_star = math.exp(opt.x)
g_star = np.sum(be * lam_p / (a_star + be * lam_p))
check("the evidence over alpha (w fixed, H = outer product) is maximised where alpha = gamma / w^T w: (5.178) is the stationarity condition", abs(a_star - g_star / (w @ w)) < 1e-4 * a_star, "argmax alpha = %.5f, gamma/w^Tw = %.5f" % (a_star, g_star / (w @ w)))
def f_beta(b):
    return -m.laplace_log_evidence(b * ED + 0.5 * al * w @ w, Hp, al, b, len(xs_), net.W)
optb = optimize.minimize_scalar(lambda lb: f_beta(math.exp(lb)), bounds=(-6, 8), method="bounded", options=dict(xatol=1e-10))
b_star = math.exp(optb.x)
g_b = np.sum(b_star * lam_p / (al + b_star * lam_p))
check("the evidence over beta (w fixed, H = outer product) is maximised where 1/beta = 2 E_D/(N - gamma): (5.180) is the stationarity condition", abs(1 / b_star - 2 * ED / (len(xs_) - g_b)) < 1e-4 / b_star, "argmax 1/beta = %.5f, 2E_D/(N-gamma) = %.5f" % (1 / b_star, 2 * ED / (len(xs_) - g_b)))
lev_hist = [h[3] for h in ev["history"]]
check("the log evidence rises along the alternating procedure (final value above the first)", lev_hist[-1] >= lev_hist[0] - 1e-9, "%.3f -> %.3f" % (lev_hist[0], lev_hist[-1]))
# predictive distribution
Xg = np.linspace(0, 1, 15)[:, None]
ym, var = m.predictive_regression(net, w, ev["A"], be, Xg)
G = net.grad_activation(w, Xg)[:, 0, :]
check("predictive variance (5.173) = 1/beta + g^T A^-1 g >= 1/beta, g = grad_w y", np.all(var >= 1 / be - 1e-12) and np.abs(var - 1 / be - np.einsum("nw,wv,nv->n", G, np.linalg.inv(ev["A"]), G)).max() < 1e-12)
fd_g = np.zeros((15, net.W))
for i in range(net.W):
    e = np.zeros(net.W)
    e[i] = 1e-6
    fd_g[:, i] = (net.output(w + e, Xg)[:, 0] - net.output(w - e, Xg)[:, 0]) / 2e-6
check("g = d y / d w from backpropagation = finite differences", np.abs(G - fd_g).max() < 1e-7)
# classification evidence and the kappa approximation
Xc, tc = m.classification_data()
netc = m.Net([2, 6, 1], "bce")
w0 = m.fit_best_of(netc, Xc, tc[:, None], np.random.default_rng(1), starts=4, sd1=1.0, maxiter=2000)
evc = m.evidence_classification(netc, Xc, tc[:, None], w0, alpha0=0.05)
alc = evc["alpha"]
check("classification evidence procedure converges: alpha = gamma / w^T w", abs(alc - evc["gamma"] / (evc["w"] @ evc["w"])) < 1e-3 * alc, "alpha = %.4f" % alc)
check("one-group gamma: sum_i (1 - alpha [A^-1]_ii) = sum lam_i/(alpha + lam_i)", abs(evc["gamma"] - np.sum(evc["eig"] / (evc["eig"] + alc))) < 1e-8)
pm, pp, s2 = m.predictive_classification(netc, evc["w"], evc["A"], Xc[:5])
mu_a = netc.forward(evc["w"], Xc[:5])["a"][:, 0]
exact_p = []
for m_, s_ in zip(mu_a, s2):
    a_ = np.linspace(m_ - 12 * math.sqrt(s_), m_ + 12 * math.sqrt(s_), 40001)
    exact_p.append(np.trapz(common.sigmoid(a_) * stats.norm.pdf(a_, m_, math.sqrt(s_)), a_))
exact_p = np.array(exact_p)
check("marginalised predictive sigma(kappa a_MAP) (5.190) is within 0.02 of the numerical integral of sigma(a) N(a | a_MAP, sigma_a^2) (4.151)", np.abs(pm - exact_p).max() < 0.02, "max |diff| = %.4f" % np.abs(pm - exact_p).max())
check("marginalising pulls probabilities toward 0.5 and leaves the p = 0.5 contour unchanged (|a| shrinks by kappa < 1)", np.all(np.abs(pm - 0.5) <= np.abs(pp - 0.5) + 1e-12) and np.all((pm > 0.5) == (pp > 0.5)))

print("== The demos' random numbers")
js = m.JSRandom(7)
x1, t1 = m.data_sine(4, 0.3, js)
net = m.Net([1, 3, 1])
w0 = m.init_weights(net, m.JSRandom(11), 1.0)
print("   data_sine(4, 0.3, JSRandom(7)) x =", np.round(x1[:, 0], 8).tolist(), "t =", np.round(t1[:, 0], 8).tolist())
print("   init_weights([1,3,1], JSRandom(11)) first 4 =", np.round(w0[:4], 8).tolist())
check("data_sine and init_weights are reproducible with JSRandom", np.array_equal(m.data_sine(4, 0.3, m.JSRandom(7))[0], x1) and np.array_equal(m.init_weights(net, m.JSRandom(11), 1.0), w0))
xa, ta = m.data_sine(10, 0.3, m.JSRandom(1))
neta = m.Net([1, 3, 1])
wa = m.init_weights(neta, m.JSRandom(2), 1.0)
wb = m.adam_train(neta, wa, xa, ta, lr=0.02, steps=300)
check("Adam training lowers the error (the demos' optimiser)", neta.error(wb, xa, ta) < neta.error(wa, xa, ta))

print("\n%d of %d checks passed" % (sum(RESULTS), len(RESULTS)))
raise SystemExit(0 if all(RESULTS) else 1)