main.py

import numpy as np
L = 3
n = [2, 3, 3, 1]

# Randomly assign the weights and biases
W1 = np.random.randn(n[1], n[0])
W2 = np.random.randn(n[2], n[1])
W3 = np.random.randn(n[3], n[2])
b1 = np.random.randn(n[1], 1)
b2 = np.random.randn(n[2], 1)
b3 = np.random.randn(n[3], 1)


def prepare_data():
    X = np.array([
        [150, 70],  # 150 pounds, 70 inches tall.
        [254, 73],
        [312, 68],
        [120, 60],
        [154, 61],
        [212, 65],
        [216, 67],
        [145, 67],
        [184, 64],
        [130, 69]
    ])
    A0 = X.T

    y = np.array([
        0,
        1,
        1,
        0,
        0,
        1,
        1,
        0,
        1,
        0
    ])
    m = 10
    Y = y.reshape(n[3], m)
    return A0, Y, m


def cost(y_hat, y):
    """
    y_hat should be a n^L x m matrix
    y should be a n^L x m matrix
    """
    # 1. losses is a n^L x m
    losses = - ((y * np.log(y_hat)) + (1 - y)*np.log(1 - y_hat))

    m = y_hat.reshape(-1).shape[0]

    # 2. summing across axis = 1 means we sum across rows,
    #   making this a n^L x 1 matrix
    summed_losses = (1 / m) * np.sum(losses, axis=1)

    # 3. unnecessary, but useful if working with more than one node
    #   in output layer
    return np.sum(summed_losses)


def g(z):
    return 1 / (1 + np.exp(-1 * z))


def feed_forward(A0):
    # layer 1 calculations
    Z1 = W1 @ A0 + b1
    A1 = g(Z1)

    # layer 2 calculations
    Z2 = W2 @ A1 + b2
    A2 = g(Z2)

    # layer 3 calculations
    Z3 = W3 @ A2 + b3
    A3 = g(Z3)

    cache = {
        "A0": A0,
        "A1": A1,
        "A2": A2
    }

    return A3, cache


A0, Y, m = prepare_data()


def backprop_layer_3(y_hat, Y, m, A2, W3):
    A3 = y_hat

    # step 1. calculate dC/dZ3 using shorthand
    dC_dZ3 = (1/m) * (A3 - Y)
    assert dC_dZ3.shape == (n[3], m)

    # step 2. calculate dC/dW3 = dC/dZ3 * dZ3/dW3
    #   we matrix multiply dC/dZ3 with (dZ3/dW3)^T
    dZ3_dW3 = A2
    assert dZ3_dW3.shape == (n[2], m)

    dC_dW3 = dC_dZ3 @ dZ3_dW3.T
    assert dC_dW3.shape == (n[3], n[2])

    # step 3. calculate dC/db3 = np.sum(dC/dZ3, axis=1, keepdims=True)
    dC_db3 = np.sum(dC_dZ3, axis=1, keepdims=True)
    assert dC_db3.shape == (n[3], 1)

    # step 4. calculate propagator dC/dA2 = dC/dZ3 * dZ3/dA2
    dZ3_dA2 = W3
    dC_dA2 = W3.T @ dC_dZ3
    assert dC_dA2.shape == (n[2], m)

    return dC_dW3, dC_db3, dC_dA2


def backprop_layer_2(propagator_dC_dA2, A1, A2, W2):

    # step 1. calculate dC/dZ2 = dC/dA2 * dA2/dZ2

    # use sigmoid derivation to arrive at this answer:
    #   sigmoid'(z) = sigmoid(z) * (1 - sigmoid(z))
    #     and if a = sigmoid(z), then sigmoid'(z) = a * (1 - a)
    dA2_dZ2 = A2 * (1 - A2)
    dC_dZ2 = propagator_dC_dA2 * dA2_dZ2
    assert dC_dZ2.shape == (n[2], m)

    # step 2. calculate dC/dW2 = dC/dZ2 * dZ2/dW2
    dZ2_dW2 = A1
    assert dZ2_dW2.shape == (n[1], m)

    dC_dW2 = dC_dZ2 @ dZ2_dW2.T
    assert dC_dW2.shape == (n[2], n[1])

    # step 3. calculate dC/db2 = np.sum(dC/dZ2, axis=1, keepdims=True)
    dC_db2 = np.sum(dC_dW2, axis=1, keepdims=True)
    assert dC_db2.shape == (n[2], 1)

    # step 4. calculate propagator dC/dA1 = dC/dZ2 * dZ2/dA1
    dZ2_dA1 = W2
    dC_dA1 = W2.T @ dC_dZ2
    assert dC_dA1.shape == (n[2], m)

    return dC_dW2, dC_db2, dC_dA1


def backprop_layer_1(propagator_dC_dA1, A1, A0, W1):

    # step 1. calculate dC/dZ1 = dC/dA1 * dA1/dZ1

    # use sigmoid derivation to arrive at this answer:
    #   sigmoid'(z) = sigmoid(z) * (1 - sigmoid(z))
    #     and if a = sigmoid(z), then sigmoid'(z) = a * (1 - a)
    dA1_dZ1 = A1 * (1 - A1)
    dC_dZ1 = propagator_dC_dA1 * dA1_dZ1
    assert dC_dZ1.shape == (n[1], m)

    # step 2. calculate dC/dW1 = dC/dZ1 * dZ1/dW1
    dZ1_dW1 = A0
    assert dZ1_dW1.shape == (n[0], m)

    dC_dW1 = dC_dZ1 @ dZ1_dW1.T
    assert dC_dW1.shape == (n[1], n[0])

    # step 3. calculate dC/db1 = np.sum(dC/dZ1, axis=1, keepdims=True)
    dC_db1 = np.sum(dC_dW1, axis=1, keepdims=True)
    assert dC_db1.shape == (n[1], 1)

    return dC_dW1, dC_db1


def train():
    # must use global keyword in order to modify global variables
    global W3, W2, W1, b3, b2, b1

    epochs = 1000  # training for 1000 iterations
    alpha = 0.1  # set learning rate to 0.1
    costs = []  # list to store costs

    for e in range(epochs):
        # 1. FEED FORWARD
        y_hat, cache = feed_forward(A0)

        # 2. COST CALCULATION
        error = cost(y_hat, Y)
        costs.append(error)

        # 3. BACKPROP CALCULATIONS

        dC_dW3, dC_db3, dC_dA2 = backprop_layer_3(
            y_hat,
            Y,
            m,
            A2=cache["A2"],
            W3=W3
        )

        dC_dW2, dC_db2, dC_dA1 = backprop_layer_2(
            propagator_dC_dA2=dC_dA2,
            A1=cache["A1"],
            A2=cache["A2"],
            W2=W2
        )

        dC_dW1, dC_db1 = backprop_layer_1(
            propagator_dC_dA1=dC_dA1,
            A1=cache["A1"],
            A0=cache["A0"],
            W1=W1
        )

        # 4. UPDATE WEIGHTS
        W3 = W3 - (alpha * dC_dW3)
        W2 = W2 - (alpha * dC_dW2)
        W1 = W1 - (alpha * dC_dW1)

        b3 = b3 - (alpha * dC_db3)
        b2 = b2 - (alpha * dC_db2)
        b1 = b1 - (alpha * dC_db1)

        if e % 20 == 0:
            print(f"epoch {e}: cost = {error:4f}")

    return costs


costs = train()