<jupyter_start><jupyter_text>Let's look at some images<jupyter_code>X,y = read_fMNIST('training')
print('Training images')<jupyter_output>Training images<jupyter_text>Now we will build the code for our softmax model<jupyter_code>def softmax(scores):
    """
    takes the softmax along the second dimension of a matrix, returns class scores
    """
    # we will adjust the dynamic range by subtracting the max, to prevent potential underflow in exp
    exp_scores = np.exp(scores-np.max(scores,axis=1,keepdims=True))+1e-6
    return exp_scores/(np.sum(exp_scores,axis=1,keepdims=True))


def log_likelihood(X,w,y):
    """
    Inputs: 
      X: a Nxd matrix where each of the N rows is a datapont with d features
      Y: a Nx10 matrix, each row is a one-hot vector with the ith entry is 1 if that datapoint belongs to class i 
      w: is a dxC matrix containing your model parameters
    Outputs:
      LL: a scalar containing the AVERAGE log-likelihood of the dataset's labels y, given inputs X and model w
    """
    #compute the un-normalized 'scores' of each class for each datapoint
    scores = np.dot(X,w)
    
    #normalize the scores to get a distribution over classes for each datapoint
    predictions = softmax(scores)
    if np.min(predictions)<= 0:
        print('zeros')
        
    #Use the predicted distributions, and the true distributions to compute log-likelihood
    LL = np.sum(y*np.log(predictions))/y.shape[0]

    return LL

def objective(X,w,y):
    """
    Compute components of the optimization objective
    Output:
    logloss : value of log-loss (negative average log-likelihood)
    regularizer: value of the regularization term (Frobenius norm of w), NOT multiplied by lambda
    """
    regularizer = np.linalg.norm(w,'fro')
    logloss = - log_likelihood(X,w,y)
    return logloss, regularizer



def gradient(X,w,y):
    """
    Compute in grad the gradient of the log-loss of the model (w) given the features (X) and labels (y) w.r.t. w
    Also compute and return in l2grad the gradient of the regularizer lambda*norm(w)
    Note: we can't use the name 'lambda' since it's a keyword in Python
    Reminder: the full objective is -log p(y|X;w) + lmbda*norm2(w)
    """
    
    #Compute any extra variable needed to compute the gradient:
    scores = np.dot(X,w)
    predictions = softmax(scores)

    #Compute the gradient of the average log-loss
    grad = -np.dot(X.T,y-predictions)/y.shape[0]
    
    l2grad = .5*np.hstack((np.zeros((w.shape[0],1)),w[:,1:w.shape[1]]))
    
    return grad, l2grad

def infer(X,w):
    """
    Compute the class scores that model w gives for each class, given data X
    """
    scores = np.dot(X,w)
    return scores

def accuracy(X,w,y):
    """
    Compute accuracy (one minus average 0/1 loss) of model w relative to true labels y on data X
    """
    y_hat = np.argmax(infer(X,w),axis=1)
    y = np.argmax(y,axis=1)
    acc = np.sum(np.equal(y,y_hat).astype(np.float32))/y.shape[0]
    return acc<jupyter_output><empty_output><jupyter_text>Now we define the SGD (stochastic gradient descent) procedure for softmax, with optional regularization<jupyter_code>def minibatch_sgd(xTr, yTr, xVl, yVl, lmbda=0, lr = .01,maxepochs=10,bsize=32,objTol=.01,verbose=1):
    
    # initialization of w
    w = np.zeros((xTr.shape[1],yTr.shape[1]))
    
    old_obj = 1e9  
    old_w = w
    
    # we will organize the run in terms of epochs (one epcoh = one full pass over data)
    # to keep track of learning curves, allocate space
    trAcc = np.zeros(maxepochs)
    obj = np.zeros(maxepochs)
    
    for epoch in range(maxepochs):
        
        #shuffle the data
        index = np.random.permutation(xTr.shape[0])
        nBatches = np.floor(xTr.shape[0]/bsize).astype(np.integer) # ignore the remainder N-floor(N/bsize), for simplicity
    
        for b in range(nBatches):
            # fill in the batch
            iBatch = index[b*bsize:(b+1)*bsize]
            xBatch = xTr[iBatch]
            yBatch = yTr[iBatch]
            
            grad,l2grad = gradient(xBatch,w,yBatch)
    
            w = w-lr*(grad+lmbda*l2grad)
            
        
        # end of an epoch: test for convergence by looking at validation
        logloss, regularizer = objective(xTr,w,yTr)
        obj[epoch] = logloss+lmbda*regularizer
                
        # also record accuracy on training
        trAcc[epoch] = accuracy(xTr,w,yTr)

        obj_gain = (old_obj-obj[epoch])/np.abs(old_obj) # improvement in training objective this epich
        
        if verbose > 0.5:
            print('Epoch %d: obj=%.4f, gain %.4f  [train acc %.4f]'%(epoch,obj[epoch],obj_gain,trAcc[epoch]))
        
        
        if epoch == maxepochs-1:
            if verbose > 0:
                print('Reached max epochs, stopping')
            return w, obj, trAcc
        
        if obj_gain < objTol:  
            lr = lr/2
            if verbose > 0:
                print('Dropping learning rate to %.4f'%lr)
        
        # update bookkeeping before going to next epoch
        old_obj = obj[epoch]<jupyter_output><empty_output><jupyter_text>Before we can start experiments, we need one more piece: feature transform. In this case we are using the "raw" features (pixel values) plus the constant term, but with normalization that applies z-scoring to each pixel.<jupyter_code># if mu and sigma are provided, they are used to normalize each pixel
# otherwise, they are computed and returned
def preprocess_data(X,mu = None,sigma = None,visualize=False):

    if visualize:
        print('VISUALIZE BEFORE NORMALIZATION')
        #Visualize example before normalization
        show_fMNIST_example(X[300])
        show_fMNIST_example(X[500])
    
    if mu is None: # need to compute normalizing stats
        compute_stats = True
        mu = X.mean(axis=0, keepdims=True)
        sigma = X.std(axis=0, keepdims=True)
        sigma[np.equal(sigma,0)]=1.0 # avoid division by zero in case of degenerate features
    else:
        compute_stats = False
        
    #Normalize the data
    X = (X-mu)/sigma
    bias_feature = np.ones((X.shape[0],1))
    X = np.concatenate([np.ones((X.shape[0],1)),X],1)

    if visualize:
        print('VISUALIZE AFTER NORMALIZATION')
        #Visualize example after normalization
        show_fMNIST_example(X[300])
        show_fMNIST_example(X[500])
    
    if compute_stats:
        return X, mu, sigma
    else:
        return X



xTr, yTr = read_fMNIST('Training')
print(str(np.min(xTr[1]))+' ... '+str(np.max(xTr[1])))

xTr, mu, sigma = preprocess_data(xTr,visualize=True)<jupyter_output>0 ... 255
VISUALIZE BEFORE NORMALIZATION<jupyter_text>We now define the main routine for our experiments: load training/validation data sets, and tune the regularization parameter (and perhaps other parameters you want to tune)<jupyter_code># run minibatch SGD for every value of lambda and record training/val accuracy
def tune_regularization(xTr,yTr,xVl,yVl,lambdas,sgd_opt):
    
    obj=dict()
    trAcc=dict()
    model=dict()
    valAcc=dict()

    for lmbda in lambdas:
        model[lmbda], obj[lmbda],trAcc[lmbda] = minibatch_sgd(xTr,yTr,xVl,yVl,lmbda, sgd_opt['lr'],sgd_opt['maxepochs'],sgd_opt['bsize'],verbose=sgd_opt['verbose'])
        valAcc[lmbda] = accuracy(xVl,model[lmbda],yVl)
        if sgd_opt['verbose'] >= 0:
            print('-------------- lambda=%.5f, val Acc = %.4f  (tr Acc = %.4f)'%(lmbda,valAcc[lmbda],trAcc[lmbda][-1]))
        
    return model, obj, trAcc, valAcc
xTr, yTr = read_fMNIST('Training')
xVl, yVl = read_fMNIST('Validation')

xTr,mu,sigma = preprocess_data(xTr,visualize=False)
xVl = preprocess_data(xVl,mu,sigma,visualize=False)

lambdas = [0,.0001,.001,.01,.1,1.0,10.0]

sgd_opt = {'bsize':16, 'maxepochs':30, 'lr':.1,'verbose':0}


model, obj, trAcc, valAcc = tune_regularization(xTr,yTr,xVl,yVl,lambdas,sgd_opt)
best=0
bestAcc=0
for lmbda in lambdas:
    if valAcc[lmbda]>=bestAcc:
        best=lmbda
        bestAcc=valAcc[lmbda]
print('Best lambda: %.4f --> %.4f  /  %.4f'%(best,bestAcc,valAcc[0]))


#Visualize the weights of the best model for the Large dataset

w=model[best]

ind = 0
fig = plt.figure(figsize=(20,10))

for cn in class_names:
    plt.subplot(2,5,ind+1)
    plt.imshow(w[-28**2:, ind].reshape(28, 28))
    plt.title("Weights Viz: "+cn)
    plt.colorbar()
    ind += 1

plt.show()

from matplotlib.ticker import FuncFormatter, MaxNLocator

#Get the names for the tickmarks of the confusion matrix
def format_fn(tick_val, tick_pos):
    if int(tick_val) in range(10):
        return class_names[int(tick_val)]
    else:
        return ''

def build_confusion_matrix(y_hat,y):
    conf_mat = np.zeros((10,10))
    
    for i in range(10):
        for j in range(10):
            #Put the number of datapoints that had class i, and predicted class j in that entry of conf matrix
            conf_mat[i,j] = np.sum(np.equal(y,i).astype(np.float32)*np.equal(y_hat,j).astype(np.float32))
    
    for i in range(10):
        #normalize all the rows to sum to 1
        conf_mat[i,:] = conf_mat[i,:]/np.sum(conf_mat[i,:])
    
    #Make plot and add tick for each class
    fig = plt.figure(figsize=(10,10))
    ax = fig.add_subplot(111)
    plt.xticks(range(10))
    plt.yticks(range(10))

    #add labels and display matrix
    ax.xaxis.set_major_formatter(FuncFormatter(format_fn))
    ax.yaxis.set_major_formatter(FuncFormatter(format_fn))
    plt.ylabel('True Label')
    plt.xlabel('Predicted Label')
    ax.imshow(conf_mat)
    plt.show()
    
    return conf_mat

#Visualize Confusion Matrix for best model
yVl_pred = np.argmax(infer(xVl,model[best]),axis=1)
conv_mat = build_confusion_matrix(yVl_pred,np.argmax(yVl,axis=1))<jupyter_output><empty_output><jupyter_text>A few things to note here. The weights visualized above may look fairly random, but they encode evidence each pixel provides for the image beloning to a particular class. Note that this evidence may be positive or negative. In some of the classes one can see the prototypical shape of the object, or some aspects of it (e.g., in trousers there is space between the pant legs, absence of which would be a strong negative indicator for that class).The confusion matrix is close to diagonal (since the accuracy is high), and the higher off-diagonal entries (corresponding to classes most often confused) are not particularly suprising: for example, one can see that different types of footweak tend to be confused, or that shirt and t-shirt can be sometimes confused.Finally, the regularization has pretty minor effect, presumably since we have a relatively large amount of data for this task. Select the best model and use it to generate predictions for test data<jupyter_code>xTe = read_fMNIST('testing')
xTe = preprocess_data(xTe,mu,sigma,visualize=False)


#Decide which model you want to use for your submission
chosen_model = model[best]

#Make predictions and write them to a csv file
final_preds = np.argmax(infer(xTe,chosen_model),axis=1)
create_submission_file('./large_submission.csv', final_preds)<jupyter_output><empty_output><jupyter_text>Now we will repeat this experiment with the small training set<jupyter_code>xTr, yTr = read_fMNIST('Training', load_small=True)
xVl, yVl = read_fMNIST('Validation')


xTr,mu,sigma = preprocess_data(xTr,visualize=False)
xVl = preprocess_data(xVl,mu,sigma,visualize=False)

lambdas = [0,.001,.01,.1,1.0,10]

sgd_opt = {'bsize':3, 'maxepochs':30, 'lr':.1,'verbose':0}



model, obj, trAcc, valAcc = tune_regularization(xTr,yTr,xVl,yVl,lambdas,sgd_opt)
best=0
bestAcc=0
for lmbda in lambdas:
    if valAcc[lmbda]>=bestAcc:
        best=lmbda
        bestAcc=valAcc[lmbda]
print('Best lambda: %.4f --> %.4f  /  %.4f'%(best,bestAcc,valAcc[0]))

#Visualize the weights of the best model for the Small dataset

w=model[best]

ind = 0
fig = plt.figure(figsize=(20,10))

for cn in class_names:
    plt.subplot(2,5,ind+1)
    plt.imshow(w[-28**2:, ind].reshape(28, 28))
    plt.title("Weights Viz: "+cn)
    plt.colorbar()
    ind += 1

plt.show()
#Visualize Confusion Matrix for best model for small dataset
yVl_pred = np.argmax(infer(xVl,model[best]),axis=1)
conv_mat = build_confusion_matrix(yVl_pred,np.argmax(yVl,axis=1))<jupyter_output><empty_output><jupyter_text>The confusion matrix is obviously farther from diagonal. The model weights (parameters)  much more closely resemble the shapes of objects in each class, an in fact have "ghosting" effects from objects in other classes. This is primarily due to the much smaller data set, and the fact that effect of each individual example on the eventual model is much more profound.<jupyter_code>xTe = read_fMNIST('testing')
xTe = preprocess_data(xTe,mu,sigma,visualize=False)


#Decide which model you want to use for your submission
chosen_model = model[best]

#Make predictions and write them to a csv file
final_preds = np.argmax(infer(xTe,chosen_model),axis=1)
create_submission_file('./small_submission.csv', final_preds)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10996537129665256, 'fraction_numerical': 0.015698345517506734, 'mean_word_length': 3.523494604942569, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 39, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '592466', 'license': 'MIT', 'n_tokens_mistral': 4565, 'n_tokens_neox': 4389, 'n_words': 1259, 'path': 'Problem Set 2/.ipynb_checkpoints/PS2-GABR-IBRAHIM-checkpoint.ipynb', 'repo_name': 'igabr/TTIC_STATML', 'size': '12995'}