#!/usr/bin/env python3
"""Een klein neuraal netwerk dat leert of een punt binnen een cirkel ligt.

Het rekenvoorbeeld uit het artikel "Deep learning" op
https://fermat360.nl/artikelen/deep-learning. De figuren en de getallen in dat
artikel komen uit dit script.

Het netwerk heeft twee invoerwaarden (de coordinaten van een punt), een
verborgen laag van acht neuronen en een uitvoerneuron: de kans dat het punt
binnen de cirkel ligt. Het script traint hetzelfde netwerk drie keer, met
dezelfde begingewichten: zonder activatiefunctie, met tanh en met ReLU. En het
controleert of backpropagation de juiste helling geeft.

Alleen de standaardbibliotheek van Python 3 is nodig:

    python3 neuraal_netwerk.py          # uitkomsten en de controle
    python3 neuraal_netwerk.py --json   # alles, ook de figuren, als JSON
"""

import json
import math
import random
import sys

PUNTEN = 200  # trainingsvoorbeelden, willekeurig in het vierkant [-1, 1] x [-1, 1]
NIEUW = 1000  # punten die het netwerk tijdens het trainen nooit ziet, om het na te meten
STRAAL = 0.8  # binnen deze cirkel: klasse 1, erbuiten: klasse 0 (ongeveer half om half)
VERBORGEN = 8  # neuronen in de verborgen laag
RONDES = 1000  # trainingsrondes: elke ronde ziet het netwerk alle voorbeelden
STAPLENGTE = 1.0
ZAAD = 3  # vaste startwaarde: elke keer dezelfde punten en begingewichten
ROOSTER = (-1.2, 0.05, 49)  # voor de figuur: van, stap, aantal punten per richting


# Activatiefuncties, elk met hun afgeleide (nodig voor backpropagation).


def geen(z):
    return z


def geen_afgeleide(z):
    return 1.0


def tanh(z):
    return math.tanh(z)


def tanh_afgeleide(z):
    return 1 - math.tanh(z) ** 2


def relu(z):
    return max(0.0, z)


def relu_afgeleide(z):
    return 1.0 if z > 0 else 0.0


def sigmoid(z):
    """Maakt van elk getal een kans tussen 0 en 1 (zonder overloop bij grote z)."""
    if z >= 0:
        return 1 / (1 + math.exp(-z))
    e = math.exp(z)
    return e / (1 + e)


ACTIVATIES = {
    "geen": (geen, geen_afgeleide),
    "tanh": (tanh, tanh_afgeleide),
    "ReLU": (relu, relu_afgeleide),
}


def maak_gegevens(aantal, zaad):
    """Willekeurige punten, elk met het juiste antwoord: 1 binnen de cirkel, 0 erbuiten."""
    toeval = random.Random(zaad)
    gegevens = []
    for _ in range(aantal):
        x = (2 * toeval.random() - 1, 2 * toeval.random() - 1)
        y = 1 if x[0] ** 2 + x[1] ** 2 < STRAAL**2 else 0
        gegevens.append((x, y))
    return gegevens


def begin_netwerk():
    """Willekeurige begingewichten, geschaald zoals Glorot en Bengio (2010) voorstellen."""
    toeval = random.Random(ZAAD + 1)

    def gewicht(n_in, n_uit):
        grens = math.sqrt(6 / (n_in + n_uit))
        return (2 * toeval.random() - 1) * grens

    w1 = [[gewicht(2, VERBORGEN) for _ in range(2)] for _ in range(VERBORGEN)]
    b1 = [0.0] * VERBORGEN
    w2 = [gewicht(VERBORGEN, 1) for _ in range(VERBORGEN)]
    b2 = 0.0
    return w1, b1, w2, b2


def vooruit(netwerk, x, activatie):
    """De voorspelling: de kans dat punt x binnen de cirkel ligt."""
    w1, b1, w2, b2 = netwerk
    # Elk verborgen neuron: gewogen som van de invoer, plus een drempel, door de activatie.
    z1 = [w[0] * x[0] + w[1] * x[1] + b for w, b in zip(w1, b1)]
    h = [activatie(z) for z in z1]
    # Het uitvoerneuron: gewogen som van de verborgen laag, door de sigmoid.
    z2 = sum(w * hi for w, hi in zip(w2, h)) + b2
    return z1, h, sigmoid(z2)


def verlies(p, y):
    """Kruisentropie: hoe verrast het netwerk is door het juiste antwoord y."""
    p = min(max(p, 1e-12), 1 - 1e-12)
    return -(y * math.log(p) + (1 - y) * math.log(1 - p))


def terug(netwerk, x, y, activatie, afgeleide):
    """Backpropagation voor een voorbeeld: de helling van het verlies naar elk gewicht.

    Van achteren naar voren, met de kettingregel. Eerst hoe het verlies
    verandert met de uitvoer, dan met elk verborgen neuron, dan met de
    gewichten daarvoor. Elke tussenuitkomst wordt een keer berekend en
    hergebruikt.
    """
    w1, b1, w2, b2 = netwerk
    z1, h, p = vooruit(netwerk, x, activatie)
    d2 = p - y  # het verlies naar z2: sigmoid en kruisentropie samen geven p - y
    g_w2 = [d2 * hi for hi in h]
    g_b2 = d2
    d1 = [d2 * w * afgeleide(z) for w, z in zip(w2, z1)]  # terug door de verborgen laag
    g_w1 = [[d * x[0], d * x[1]] for d in d1]
    g_b1 = d1
    return (g_w1, g_b1, g_w2, g_b2), verlies(p, y)


def helling(netwerk, gegevens, activatie, afgeleide):
    """De helling van het gemiddelde verlies over alle voorbeelden, en dat verlies zelf."""
    n = len(gegevens)
    g_w1 = [[0.0, 0.0] for _ in range(VERBORGEN)]
    g_b1 = [0.0] * VERBORGEN
    g_w2 = [0.0] * VERBORGEN
    g_b2 = 0.0
    totaal = 0.0
    for x, y in gegevens:
        (a_w1, a_b1, a_w2, a_b2), v = terug(netwerk, x, y, activatie, afgeleide)
        for i in range(VERBORGEN):
            g_w1[i][0] += a_w1[i][0] / n
            g_w1[i][1] += a_w1[i][1] / n
            g_b1[i] += a_b1[i] / n
            g_w2[i] += a_w2[i] / n
        g_b2 += a_b2 / n
        totaal += v / n
    return (g_w1, g_b1, g_w2, g_b2), totaal


def train(naam, gegevens):
    """RONDES stappen gradientdaling; geeft het getrainde netwerk en het verloop van het verlies."""
    activatie, afgeleide = ACTIVATIES[naam]
    w1, b1, w2, b2 = begin_netwerk()
    verloop = []
    for _ in range(RONDES):
        (g_w1, g_b1, g_w2, g_b2), v = helling((w1, b1, w2, b2), gegevens, activatie, afgeleide)
        verloop.append(v)
        w1 = [[w - STAPLENGTE * g for w, g in zip(rij, grij)] for rij, grij in zip(w1, g_w1)]
        b1 = [b - STAPLENGTE * g for b, g in zip(b1, g_b1)]
        w2 = [w - STAPLENGTE * g for w, g in zip(w2, g_w2)]
        b2 = b2 - STAPLENGTE * g_b2
    return (w1, b1, w2, b2), verloop


def nauwkeurigheid(netwerk, gegevens, activatie):
    """Het aandeel voorbeelden waarbij het netwerk aan de goede kant van 50% zit."""
    goed = sum(1 for x, y in gegevens if (vooruit(netwerk, x, activatie)[2] > 0.5) == (y == 1))
    return goed / len(gegevens)


def controleer_backpropagation(gegevens, epsilon=1e-6):
    """Vergelijkt backpropagation met de trage methode: elk gewicht een beetje verschuiven.

    Geeft het grootste relatieve verschil tussen de twee, over alle gewichten
    van het tanh-netwerk in zijn beginstand.
    """
    activatie, afgeleide = ACTIVATIES["tanh"]
    netwerk = begin_netwerk()
    (g_w1, g_b1, g_w2, g_b2), _ = helling(netwerk, gegevens, activatie, afgeleide)

    def gemiddeld_verlies(net):
        return sum(verlies(vooruit(net, x, activatie)[2], y) for x, y in gegevens) / len(gegevens)

    def verschoven(pad, delta):
        w1, b1, w2, b2 = netwerk
        w1 = [rij[:] for rij in w1]
        b1, w2 = b1[:], w2[:]
        if pad[0] == "w1":
            w1[pad[1]][pad[2]] += delta
        elif pad[0] == "b1":
            b1[pad[1]] += delta
        elif pad[0] == "w2":
            w2[pad[1]] += delta
        else:
            b2 += delta
        return w1, b1, w2, b2

    paden = [("w1", i, j) for i in range(VERBORGEN) for j in range(2)]
    paden += [("b1", i) for i in range(VERBORGEN)] + [("w2", i) for i in range(VERBORGEN)] + [("b2",)]
    snel = [g_w1[p[1]][p[2]] for p in paden[: 2 * VERBORGEN]]
    snel += g_b1 + g_w2 + [g_b2]
    grootste = 0.0
    for pad, s in zip(paden, snel):
        traag = (gemiddeld_verlies(verschoven(pad, epsilon)) - gemiddeld_verlies(verschoven(pad, -epsilon))) / (2 * epsilon)
        grootste = max(grootste, abs(s - traag) / max(abs(s), abs(traag), 1e-12))
    return len(paden), grootste


def main():
    gegevens = maak_gegevens(PUNTEN, ZAAD)
    nieuw = maak_gegevens(NIEUW, ZAAD + 2)
    uitkomst = {}
    for naam in ACTIVATIES:
        netwerk, verloop = train(naam, gegevens)
        activatie = ACTIVATIES[naam][0]
        van, stap, aantal = ROOSTER
        rooster = [
            [round(vooruit(netwerk, (van + stap * k, -van - stap * r), activatie)[2], 4) for k in range(aantal)]
            for r in range(aantal)
        ]
        uitkomst[naam] = {
            "nauwkeurigheid": nauwkeurigheid(netwerk, gegevens, activatie),
            "nauwkeurigheid_nieuw": nauwkeurigheid(netwerk, nieuw, activatie),
            "verloop": [round(v, 6) for v in verloop],
            "rooster": rooster,
        }
    gewichten, verschil = controleer_backpropagation(gegevens)
    if "--json" in sys.argv:
        curves = {
            naam: [[z / 10, f(z / 10)] for z in range(-30, 31)]
            for naam, f in (("sigmoid", sigmoid), ("tanh", tanh), ("ReLU", relu))
        }
        json.dump(
            {
                "punten": PUNTEN,
                "nieuw": NIEUW,
                "straal": STRAAL,
                "verborgen": VERBORGEN,
                "rondes": RONDES,
                "staplengte": STAPLENGTE,
                "rooster": {"van": ROOSTER[0], "stap": ROOSTER[1], "aantal": ROOSTER[2]},
                "gegevens": gegevens,
                "netwerken": uitkomst,
                "controle": {"gewichten": gewichten, "grootste_verschil": verschil},
                "activaties": curves,
            },
            sys.stdout,
            separators=(",", ":"),
        )
        print()
        return
    binnen = sum(y for _, y in gegevens)
    print(f"{PUNTEN} trainingspunten, {binnen} binnen de cirkel; netwerk 2-{VERBORGEN}-1, {RONDES} rondes")
    print(f"'goed nieuw': op {NIEUW} punten die het netwerk tijdens het trainen niet zag")
    print(f"  {'activatie':<10}{'verlies begin':>14}{'verlies eind':>13}{'goed':>8}{'goed nieuw':>12}")
    for naam, u in uitkomst.items():
        print(
            f"  {naam:<10}{u['verloop'][0]:>14.3f}{u['verloop'][-1]:>13.3f}"
            f"{u['nauwkeurigheid']:>8.1%}{u['nauwkeurigheid_nieuw']:>12.1%}"
        )
    print()
    print(f"Controle backpropagation, {gewichten} gewichten:")
    print(f"  grootste relatieve verschil met verschuiven: {verschil:.1e}")


if __name__ == "__main__":
    main()