<jupyter_start><jupyter_text>Build CNN  1 Import data<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import random
from scipy.io import loadmat
from scipy.optimize import minimize
from sklearn.preprocessing import OneHotEncoder
data = loadmat("ex4data1.mat")
data
X = data["X"]
y = data["y"]
X.shape, y.shape
encoder = OneHotEncoder(sparse=False, categories="auto")
y_onehot = encoder.fit_transform(y)
y_onehot.shape
import random
#将400 转为 20*20图像
def toImg(row):
    width, height = 20, 20
    square = row.reshape(width, height)
    return square.T

#展示图像
def displayImg(X, indices_to_display = None):
    width, height = 20, 20
    nrows, ncols = 10, 10
    if indices_to_display == None:
        indices_to_display = random.sample(range(X.shape[0]), nrows * ncols)
    big_picture = np.zeros((height * nrows, width * ncols))
    irow , icol = 0, 0
    for idx in indices_to_display:
        if icol == ncols:
            irow += 1
            icol = 0
        iimg = toImg(X[idx])
        big_picture[(irow*height):(irow*height+iimg.shape[0]),
                   (icol*width):(icol*width+iimg.shape[1])] = iimg
        icol += 1
    fig = plt.figure(figsize=(6, 6))
    plt.imshow(big_picture, interpolation='nearest')
displayImg(X)<jupyter_output><empty_output><jupyter_text>2 Convolution<jupyter_code>def convolution(Xi, kernal): #(20, 20) (5, 5)
    xrow, xcol = Xi.shape
    krow, kcol = kernal.shape
    rrow, rcol = xrow - krow + 1, xcol - kcol + 1
    res = np.zeros((rrow, rcol))    
    for i in range(rrow):
        for j in range(rcol):
            res[i][j] = np.sum(Xi[i:i+krow, j:j+kcol] * kernal)
    return res # (16, 16)<jupyter_output><empty_output><jupyter_text>3 Pooling<jupyter_code>def maxPooling(conv):
    # use 2*2 pooling
    row, col = conv.shape
    res = np.zeros((row // 2, col // 2))
    pooling_grad = np.zeros(conv.shape)
    for i in range(0, row, 2):
        for j in range(0, col, 2):
            m = conv[i, j]
            mr, mc = i, j
            res[i//2,j//2] = np.max(conv[i:i+2,j:j+2])
            idx = int(np.argmax(conv[i:i+2,j:j+2]))
            pooling_grad[i+idx//2,j+idx%2] = 1
    return res, pooling_grad
theta = np.random.rand(5,5)
theta.shape
Xi = X[0].reshape((20,20))
%%time
for i in range(100):
    res, grad = maxPooling(convolution(Xi, theta))<jupyter_output><empty_output><jupyter_text>4 Sigmod<jupyter_code>def sigmod(z):
    return 1 / (1 + np.exp(-z))<jupyter_output><empty_output><jupyter_text>5 Initializate thetas<jupyter_code>input_size = 400
inrow, incol = 20, 20
num_units = 5
krow, kcol = 5, 5
prow, pcol = 2, 2
conrow, concol = inrow - krow + 1, incol - kcol + 1
hrow, hcol = (inrow - krow + 1) // 2, (incol - kcol + 1) // 2
output_size = 10

weights_size = (krow * kcol + 1 + # convolution layer
                hrow * hcol * output_size) * num_units + output_size # output layer

params = (np.random.random(size=weights_size) - 0.5) * 0.25
params.shape

# inrow, incol, num_units, krow, kcol, cnrow, cncol, hrow, hcol, output_size
def encode(theta1, theta2, b1, b2):
    return np.concatenate((theta1.flatten(), theta2.flatten(), b1.flatten(), b2.flatten()))
def decode(params, krow, kcol, hrow, hcol, num_units, output_size):
    theta1 = params[: krow * kcol * num_units].reshape((num_units, krow, kcol)) # (5, 5, 5)
    theta2 = params[krow * kcol * num_units: 
                   krow * kcol * num_units + 
                    hrow * hcol * output_size * num_units].reshape((num_units, output_size, hrow, hcol)) # (5, 10, 8,8)
    b1 = params[-num_units - output_size: - output_size].reshape((num_units, 1))
    b2 = params[-output_size:].reshape((output_size, 1))
    return theta1, theta2, b1, b2
t1, t2, b1, b2 = decode(params,krow, kcol, hrow, hcol, num_units, output_size)
t1.shape, t2.shape, b1.shape, b2.shape


len(encode(t1, t2, b1, b2))<jupyter_output><empty_output><jupyter_text>6 Forward propagate<jupyter_code>def forwardPropagate(Xi, theta1, theta2, b1, b2):
    a1 = Xi.reshape(inrow, incol) # (20, 20)
    z2 = np.zeros((num_units, conrow, concol)) # (5, 16, 16) 
    a2 = np.zeros((num_units, conrow, concol)) # (5, 16, 16)
    pooling_grad = np.zeros((num_units, conrow, concol)) # (5, 16, 16)
    a3 = np.zeros((num_units, hrow, hcol)) # (5, 8, 8) z3 = a3
    z4 = np.zeros((output_size, 1)) # (10, 1)
    a4 = np.zeros((output_size, 1)) # (10, 1)
    for i in range(num_units):
        z2[i] = convolution(a1, theta1[i]) + b1[i] # (16, 16)
    a2 = sigmod(z2)
    for i in range(num_units):
        a3[i], pooling_grad[i] = maxPooling(a2[i]) # (8, 8) (16, 16)
        for j in range(output_size):
            z4[j] += np.sum(a3[i] * theta2[i,j])
    z4 += b2
    a4 = sigmod(z4)
    return a1, z2, a2, pooling_grad, a3, z4, a4
a1, z2, a2, pooling_grad, a3, z4, a4  = forwardPropagate(X[0], t1, t2, b1, b2)
temp = a3[0].reshape(1, *a3[0].shape)
%%time 
for i in range(1000):
    a = forwardPropagate(X[0], t1, t2, b1, b2)<jupyter_output><empty_output><jupyter_text>7 Predict and compute accuracy<jupyter_code>def computeAccuracy(X, y, theta1, theta2, b1, b2):
    m = X.shape[0]
    correct = 0
    for i in range(m):
        a1, z2, a2, pooling_grad, a3, z4, a4 = forwardPropagate(X[i], theta1, theta2, b1, b2)
        #h = a3.flatten()
        correct += int(np.argmax(a4) + 1 == y[i])
    return "accuracy = {0}%".format(100 * correct / m)
computeAccuracy(X[:500], y, t1, t2, b1, b2)<jupyter_output><empty_output><jupyter_text>8 Cost function<jupyter_code>def cost(params, X, y, learning_rate=0.):
    m = X.shape[0]
    theta1, theta2, b1, b2 = decode(params, krow, kcol, hrow, hcol, num_units, output_size)
    J = 0.
    for i in range(m):
        a1, z2, a2, pooling_grad, a3, z4, a4 = forwardPropagate(X[i], theta1, theta2, b1, b2)
        tempy = y[i].reshape((output_size, 1)) # (10, 1)
        first_term = tempy * np.log(a4) # (10, 1)
        second_term = (1 - tempy) * np.log(1 - a4) # (10, 1)
        J += -np.sum(first_term + second_term)
    J /= m
    J += (float(learning_rate) / (2 * m)) * \
        (np.sum(theta1 **2) + np.sum(theta2 ** 2))
    return J
cost(params, X[:50,:], y_onehot)<jupyter_output><empty_output><jupyter_text>9 Sigmod gradient<jupyter_code>def sigmodGradient(z):
    t = sigmod(z)
    return t * (1 - t)<jupyter_output><empty_output><jupyter_text>10 Backpropagation<jupyter_code>def backPropagate(params, inrow, incol, num_units, krow, kcol, 
                  conrow, concol, hrow, hcol, output_size, X, y, learning_rate=0.):
    m = X.shape[0]
    theta1, theta2, b1, b2 = decode(params, krow, kcol, hrow, hcol, num_units, output_size)
    J = 0.
    delta1 = np.zeros(theta1.shape) # (5, 5, 5)
    delta2 = np.zeros(theta2.shape) # (5, 10, 8, 8)
    
    belta1 = np.zeros(b1.shape) # (5, 1)
    belta2 = np.zeros(b2.shape) # (10, 1)
    
    for i in range(m):
        a1, z2, a2, pooling_grad, a3, z4, a4 = forwardPropagate(X[i], theta1, theta2, b1, b2)
        tempy = y[i,:].reshape((output_size, 1)) # (10, 1)
        first_term = tempy * np.log(a4) # (10, 1)
        second_term = (1 - tempy) * np.log(1 - a4) # (10, 1)
        J += -np.sum(first_term + second_term)
        dt2 = a4 - tempy  # (10, 1)
        belta2 += dt2  # det b2
        for j in range(num_units):
            temp = np.zeros((output_size, hrow, hcol)) # (10, 8, 8)
            for k in range(output_size):
                temp[k] = a3[j].copy()
            delta2[j] += dt2.reshape((output_size, 1, 1)) * temp  # (10, 8, 8) det w2
                        
            dt1 = np.zeros((conrow, concol))
            for p in range(conrow):
                for q in range(concol):
                    if pooling_grad[j,p,q] == 1.:
                        val = np.zeros((output_size, 1))
                        p1 = p // 2
                        q1 = p // 2
                        for k in range(output_size):
                            val[k,0] = theta2[j,k,p1,q1]
                        
                        dt1[p, q] = np.sum(val * dt2) * z2[j,p,q] * (1 - z2[j,p,q])

            belta1[j] += np.sum(dt1) # det b1
            for p in range(krow):
                for q in range(kcol):
                    delta1[j, p, q] += np.sum(a1[p:p+conrow, q:q+concol] * dt1) # det w1
    J /= m
    delta1 /= m
    delta2 /= m
    belta1 /= m
    belta2 /= m
        
    J += (float(learning_rate) / (2 * m)) * (np.sum(theta1 ** 2) + np.sum(theta2 ** 2))
    delta1 += theta1 * learning_rate / m
    delta2 += theta2 * learning_rate / m
    return J, encode(delta1, delta2, b1, b2)
J, grad = backPropagate(params, inrow, incol, num_units, krow, kcol, 
                  conrow, concol, hrow, hcol, output_size, X[:50], y_onehot)
J<jupyter_output><empty_output><jupyter_text>11 Gradient checking<jupyter_code>def checkGradient(params, grad, krow, kcol, 
                  hrow, hcol, num_units, output_size, X, y, leraning_rate=0.):
    eps = 1e-4
    n = len(params)
    for i in range(10):
        x = random.randint(0, n - 1)
        epsvec = np.zeros(n)
        epsvec[x] = eps
        cost_high = cost(params + epsvec, X, y_onehot)
        cost_low = cost(params - epsvec, X, y_onehot)
        mygrad = (cost_high - cost_low) / (2 * eps)
        print(f"Element: {x}. Numerical grad = {mygrad} BP grad = {grad[x]}")

checkGradient(params, grad, krow, kcol, hrow, hcol, num_units, output_size, X[:50], y_onehot)
#len(grad)
learning_rate = 0.03<jupyter_output><empty_output><jupyter_text>12 Training<jupyter_code>%%time
fmin = minimize(fun=backPropagate, x0=params,
               args=(inrow, incol, num_units, krow, kcol, 
                  conrow, concol, hrow, hcol, output_size, X[:50], y_onehot[:50], learning_rate),
               method="TNC", jac=True, options={"maxiter": 50})
fmin<jupyter_output><empty_output><jupyter_text>13 Compute accuracy<jupyter_code>def ComputeAccuracy(params, X, y):
    t1, t2, b1, b2 = decode(params, krow, kcol, hrow, hcol, num_units, output_size)
    m = X.shape[0]
    correct = 0
    for i in range(m):
        a1, z2, a2, pooling_grad, a3, z4, a4  = forwardPropagate(X[0], t1, t2, b1, b2)
        correct += int(np.argmax(a4) + 1 == y[i])
    print("accuracy = {0}%".format(100 * correct / m))<jupyter_output><empty_output><jupyter_text>14 Mini data set<jupyter_code>tr_X = []
tr_y_onehot = []
tr_y = []

ts_X = []
ts_y_onehot = []
ts_y = []
m = X.shape[0]
n = m // 100
for i in range(0, m, m // 10):
    tr_X.append(X[i:i+n])
    tr_y.append(y[i:i+n])
    tr_y_onehot.append(y_onehot[i:i+n])
    
    ts_X.append(X[i+n:i+2*n])
    ts_y.append(y[i+n:i+2*n])
    ts_y_onehot.append(y_onehot[i+n:i+2*n])
    
def connect(a, b):
    return np.concatenate((a.ravel(), b.ravel()))
from functools import reduce

res = reduce(connect, tr_X)
tr_X = res.reshape((500, 400))

res = reduce(connect, ts_X)
ts_X = res.reshape((500, 400))
res = reduce(connect, tr_y)
tr_y = res.reshape((500,1))
res = reduce(connect, ts_y)
ts_y = res.reshape((500,1))
res = reduce(connect, tr_y_onehot)
tr_y_onehot = res.reshape((500,10))
res = reduce(connect, ts_y_onehot)
ts_y_onehot = res.reshape((500,10))
def connect(a, b):
    return np.concatenate((a.ravel(), b.ravel()))<jupyter_output><empty_output><jupyter_text>from functools import reduce<jupyter_code>res = reduce(connect, tr_X)
tr_X = res.reshape((500, 400))
res = reduce(connect, ts_X)
ts_X = res.reshape((500, 400))
res = reduce(connect, tr_y)
tr_y = res.reshape((500,1))
res = reduce(connect, ts_y)
ts_y = res.reshape((500,1))
res = reduce(connect, tr_y_onehot)
tr_y_onehot = res.reshape((500,10))
res = reduce(connect, ts_y_onehot)
ts_y_onehot = res.reshape((500,10))
%%time
fmin = minimize(fun=backPropagate, x0=params,
               args=(inrow, incol, num_units, krow, kcol, 
                  conrow, concol, hrow, hcol, output_size, tr_X, tr_y_onehot, 0.05),
               method="TNC", jac=True, options={"maxiter": 100})
fmin
ComputeAccuracy(fmin.x, tr_X, tr_y)
def predict(Xi, t1, t2, b1, b2):
    *a, h = forwardPropagate(Xi, t1, t2, b1, b2)
    return np.argmax(h) + 1

predict( tr_X[340], *decode(fmin.x, krow, kcol, hrow, hcol, num_units, output_size))
predict( tr_X[400], *decode(params, krow, kcol, hrow, hcol, num_units, output_size))
def compute(X, y, weights):
    m = X.shape[0]
    c = 0
    t1, t2, b1, b2 = decode(weights, krow, kcol, hrow, hcol, num_units, output_size)
    for i in range(m):
        res = predict( X[i], t1, t2, b1, b2)
        c += (res == y[i])[0]
    return c / m
compute(ts_X, ts_y, params)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15875299760191847, 'fraction_numerical': 0.04668265387689848, 'mean_word_length': 2.978060413354531, 'pattern_counts': {'":': 2, '<': 61, '<?xml version=': 0, '>': 61, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '533612', 'license': 'Apache-2.0', 'n_tokens_mistral': 5482, 'n_tokens_neox': 5138, 'n_words': 1168, 'path': 'StudyNotesOfML/4. Backforward /CNN.ipynb', 'repo_name': 'yang-233/ML', 'size': '12510'}