<jupyter_start><jupyter_text>MA 544 Spring 2021: Programming Assignment - 4---<jupyter_code>import os
import numpy as np
import pandas as pd
from glob import glob
import matplotlib.pyplot as pt
from typing import Iterable, List<jupyter_output><empty_output><jupyter_text>Question 1---Find the nonnegative basis representation of images from one of the following databases. Show a comparison of some sample images with their reconstruction from the basis.>- The ORL database of faces at [ORL link](http://cam-orl.co.uk/facedatabase.html/:~:text=Our%20Database%20of%20Faces%2C%20%28formerly%20%27The%20ORL%20Database,Robotics%20Group%20of%20the%20Cambridge%20University%20Engineering%20Department.).>- The CBCL database of faces at [MIT link](http://cbcl.mit.edu/software-datasets/heisele/facerecognition-database.html).>- Yale faces B facial images at [UCSD](http://vision.ucsd.edu/~leekc/ExtYaleDatabase/ExtYaleB.html).<jupyter_code>def list_files(directory: str, search: Iterable) -> List[str]:
    """
    list files
    """
    if isinstance(search, str):
        search = [search]

    output: List[str] = []
    for term in search:
        for filepath in glob(os.path.join(directory, term)):
            output.append(filepath)

    return output
# to avoid having to write code to fully parse the data by subject I 
# will write it to handle a single subject, yaleB01
data_folder: str = os.path.join('data', 'CroppedYale', 'yaleB01')
class PGMReader():
    def __init__(self, filepath, mode):
        assert filepath.split('.')[-1] == 'pgm'
        self.filepath = filepath
        self.mode = mode
        
    def __enter__(self):
        # return the file object after opening and setting it up
        self.pgm_image = open(self.filepath, self.mode)
        return self.read_pgm(self.pgm_image)
    
    def __exit__(self, type, value, traceback):
        self.pgm_image.close()
    
    def read_pgm(self, pgmf):
        """Return a raster of integers from a PGM as a list of lists."""
        assert pgmf.readline() == b'P5\n'
        (width, height) = [int(i) for i in pgmf.readline().split()]
        depth = int(pgmf.readline())
        assert depth <= 255

        raster = []
        for y in range(height):
            row = []
            for y in range(width):
                row.append(ord(pgmf.read(1)))
            raster.append(row)
        return np.array(raster)
m, n = 192, 168
filepaths = list_files(data_folder, '*.pgm')
images = np.zeros((m*n, len(filepaths)))
for i, filepath in enumerate(filepaths):
    with PGMReader(filepath, 'rb') as image:
        images[:, i] = image.reshape((m*n,))
        
print(f"Shape of images: ({m}, {n})")
pt.imshow(images[:, 2].reshape((m, n)), cmap='gray')
pt.imshow(images[:100, :], cmap='gray')
class LSNMF():
    def __init__(self, M):
        self.m, self.n = M.shape
        self.M = M
    
    def fit(self, threshold = 5e3, max_iter = 200):
        # initialize a random start state for W and H
        W = np.random.rand(self.m, self.n)
        H = np.random.rand(self.n, self.n)
        
        # @ delineates matrix multiplication
        # np.multiply delineates elementwise mutliplication
        # np.divide delineates elementwise division
        for i in range(max_iter):
            norm = np.linalg.norm(self.M - W @ H)
            if i % 50 == 0 or i == max_iter:
                print(f'Norm for iteration {i}: {norm}')
            if norm < threshold:
                print("Reached convergence through threshold")
                break
                
            wktv = np.dot(W.T, self.M)
            wktwk = np.dot(W.T, W)
            wktwkhk = np.dot(wktwk, H)
            
            H_kp1 = np.multiply(H, np.divide(wktv, wktwkhk))
            
            vhkp1t = np.dot(self.M, H_kp1.T)
            wkhkp1 = np.dot(W, H_kp1)
            wkhkp1hkp1t = np.dot(wkhkp1, H_kp1.T)
            
            W_kp1 = np.multiply(W, np.divide(vhkp1t, wkhkp1hkp1t))
            
            H = H_kp1
            W = W_kp1
            
            if i == max_iter-1:
                print("Reached convergence through max_iter")
        
        return W, H
lsnmf = LSNMF(images)
W, H = lsnmf.fit()<jupyter_output>Norm for iteration 0: 123871.53692533601
Norm for iteration 50: 14007.521071852227
Norm for iteration 100: 10615.854694153719
Norm for iteration 150: 8693.851355028339
Reached convergence through max_iter<jupyter_text>Original Image<jupyter_code># original image 2
pt.imshow(images[:, 2].reshape(m, n), cmap='gray')<jupyter_output><empty_output><jupyter_text>NMF Reconstruction using the Lee-Seung Algorithm<jupyter_code>combinations = H[:, 2]
reconstruction = W @ combinations
pt.imshow(reconstruction.reshape(m, n), cmap='gray')<jupyter_output><empty_output><jupyter_text>Question 2---Set up a linear regression model for the miles per gallon on the data at automobile [UCI](https://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/).  Discard the categorical data.>1. Analyze the data to get relevant insight.>1. Get feature matrix X, and target variable y.>1. Split data into training and testing.>1. Normalize data using MinMaxScaler.>1. Creat a LinearRegression object for modeling.>1. Train the model with training data.>1. Look at R^2 score for the goodness of fit for the train and test data.>1. Present a graphical comparison of true and observed responses for the test data.>1. Improve the performance of your model on the test data.<jupyter_code>from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import train_test_split
data_folder: str = os.path.join('data')
data = pd.read_csv(os.path.join(data_folder,'auto-mpg.data'),
                                 delim_whitespace=True,
                                 header=None)
column_names = ["mpg", "cylinders", "displacement", "horsepower", "weight", "acceleration", "model_year", "origin", "car_name"]
data.columns = column_names
data
data.info()
data.describe()
def is_float(x):
    try:
        float(x)
    except ValueError:
        return False
    return True
data = data[data['horsepower'].apply(lambda x: is_float(x))]
names = data.pop('car_name')
print(data.head())
data = data.astype(np.float)
scaler = MinMaxScaler()
scaler.fit(data)
data = scaler.transform(data)
y = data[:, 0]
X = data[:, 1:]
class LinearRegression():
    def __init__(self, validation_split: float = 0.2, reg: str = 'l2', lam: float = 1e-2):
        assert 0. <= validation_split <= 1.
        assert reg in ['none', 'l2']
        if reg != 'none': 
            assert lam is not None
        self.reg = reg
        self.lam = lam
        self.validation_split = validation_split
        
    def fit(self, x_in, y_in):
        x_in = np.asarray(x_in)
        y_in = np.asarray(y_in)
        X, x_test, y, y_test = train_test_split(x_in, y_in, test_size=self.validation_split)
        
        X = np.hstack([X, np.ones((X.shape[0], 1))])
        if self.reg == 'none':
            self.coef_ = np.linalg.inv(X.T @ X) @ np.dot(X.T, y) 
        if self.reg == 'l2':
            self.coef_ = np.linalg.inv(X.T @ X + self.lam*np.eye(X.shape[1])) @ np.dot(X.T, y)

        self.score_(x_test, y_test)
    
    def predict(self, features):
        features = np.hstack([features, np.array([1])])
        return np.dot(self.coef_, features)
    
    def score_(self, x_test, y_test):
        y_pred = []
        for x in x_test:
            y_pred.append(self.predict(x))
            
        self.r2 = np.corrcoef(y_test, y_pred)[0, 1]**2
linreg = LinearRegression(validation_split=0.1, reg='none')
linreg.fit(X, y)<jupyter_output><empty_output><jupyter_text>Results<jupyter_code>print(f"Calculated Coefficients [C0, C1, ..., Cn] (C0 = intercept):\n {linreg.coef_}")
print(f"Calculated r2 on the test set: {linreg.r2}")<jupyter_output>Calculated r2 on the test set: 0.8279471018412106<jupyter_text>Below, our predictions v true values are shown<jupyter_code>predictions = []
for x in X:
    predictions.append(linreg.predict(np.hstack(x)))
pt.plot(range(len(X)), predictions, 'go')
pt.plot(range(len(X)), y, 'ro')

pt.legend(["predictions", "true values"])
pt.show()<jupyter_output><empty_output><jupyter_text>Custom linear regression vs sklearn benchmark<jupyter_code>from sklearn.linear_model import Ridge
model = Ridge(alpha=1, fit_intercept=True)
model.fit(X, y)
print(f"Calculated Coefficients:\n {model.coef_}")
bench = []
for row in X:
    bench.append(model.predict(row.reshape(1, -1)))
pt.plot(range(len(X)), predictions, 'go')
pt.plot(range(len(X)), bench, 'bo')
pt.legend(["predictions", "sklearn predictions"])
pt.show()<jupyter_output><empty_output><jupyter_text>Improve through adding regularization<jupyter_code>linreg2 = LinearRegression(validation_split=0.1, reg='l2', lam=1.)
linreg2.fit(X, y)
print(f"Calculated Coefficients [C0, C1, ..., Cn] (C0 = intercept):\n {linreg2.coef_}")
print(f"Calculated r2 on the test set: {linreg2.r2}")<jupyter_output>Calculated r2 on the test set: 0.8300385504557488<jupyter_text>Below, our predictions v true values are shown<jupyter_code>predictions2 = []
for x in X:
    predictions2.append(linreg2.predict(np.hstack(x)))
pt.plot(range(len(X)), predictions2, 'go')
pt.plot(range(len(X)), y, 'ro')

pt.legend(["improved predictions", "true values"])
pt.show()<jupyter_output><empty_output><jupyter_text>Custom linear regression vs sklearn benchmark<jupyter_code>y_pred2 = []
pt.plot(range(len(X)), predictions2, 'go')
pt.plot(range(len(X)), bench, 'bo')
pt.legend(["improved predictions", "sklearn predictions"])
pt.show()<jupyter_output><empty_output><jupyter_text>Question 3---Setup a logistic regression model on the data at [adultUCI](http://archive.ics.uci.edu/ml/datasets/Adult).  Discuss the performance of your model using appropriate statistics. Use dummy variables to handle categorical variables. >1. Prepare the data. Create dummy variables for categorical variables. [See this](https://queirozf.com/entries/one-hot-encoding-a-feature-on-a-pandas-dataframe-an-example)>1. Analyze the data to get relevant insight.>2. Get feature matrix X, and target variable y (>50k or <50k)>3. Split data into training and testing>4. Normalize data using MinMaxScaler>5. Creat a LogisticRegression object for modeling>6. Train the model with training data>7. Compare the precision, recall, and F1-score on the train and test data.>1. Improve the performance of your model on the test data.  As an interesting thought experiment, a logistic regression can be thought of as a dense network with num_features input nodes, num_classes output nodes, and a sigmoid activation. I will build this network using the Tensorflow API<jupyter_code>data = pd.read_csv(os.path.join(data_folder, 'adult.data'), header=None)
orig_cols = ['age', 'workclass', 'fnlwgt', 'education', 'education_num', 'marital_status', 'occupation', 'relationship', 'race', 'sex', 'capital_gain', 'capital_loss', 'hours_per_week', 'native_country', 'target']
data.columns = orig_cols
data.head()
import tensorflow as tf
from tensorflow.data import Dataset
from tensorflow.keras.layers import Dense, Input
from tensorflow.keras.models import Model
data = pd.concat([data, pd.get_dummies(data["workclass"], prefix="wc")], axis=1)
data = pd.concat([data, pd.get_dummies(data["education"], prefix="edu")], axis=1)
data = pd.concat([data, pd.get_dummies(data["marital_status"], prefix="ms")], axis=1)
data = pd.concat([data, pd.get_dummies(data["occupation"], prefix="occ")], axis=1)
data = pd.concat([data, pd.get_dummies(data["relationship"], prefix="rel")], axis=1)
data = pd.concat([data, pd.get_dummies(data["race"], prefix="rac")], axis=1)
data = pd.concat([data, pd.get_dummies(data["sex"], prefix="s")], axis=1)
data = pd.concat([data, pd.get_dummies(data["native_country"], prefix="s")], axis=1)

data.target = data.target.astype('category')
data.target = data.target.cat.codes

for colname in ["workclass", "education", "marital_status", "occupation", "relationship", "race", "sex", "native_country"]:
    data.drop([colname], axis=1, inplace=True)
data.head(1)
y = data.pop('target').values
X = data.values
scaler = MinMaxScaler()
scaler.fit(X)
X = scaler.transform(X)
x_train, x_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
class LogisticRegression(tf.keras.Model):
    def __init__(self, num_features: int = None, reg = None):
        assert num_features is not None and num_features >= 1 and isinstance(num_features, int)
        assert len(reg) == 2 and 0. <= reg[0] and 0. <= reg[1]
        
        super(LogisticRegression, self).__init__()
        self.inp = Input((num_features,))
        self.outputs = Dense(1, activation='sigmoid', kernel_regularizer=tf.keras.regularizers.L1L2(l1=reg[0], l2=reg[1]))
    
    def call(self, features):
        return self.outputs(features)
logistic = LogisticRegression(num_features=X.shape[1], reg = [0.001, 0.001])
logistic.compile(optimizer='rmsprop',
                 loss='binary_crossentropy',
                 metrics=[tf.keras.metrics.Precision(), 
                          tf.keras.metrics.Recall(), 
                          tf.keras.metrics.TruePositives(), 
                          tf.keras.metrics.TrueNegatives(), 
                          tf.keras.metrics.FalsePositives(), 
                          tf.keras.metrics.FalseNegatives(), 
                          tf.keras.metrics.AUC()])

print(f'random output for initialization = {logistic(X[0, :].reshape(1, -1))}')
logistic.summary()
history = logistic.fit(x_train, y_train, batch_size=32, epochs=5)
history.history.keys()
precision = history.history['precision'][-1]
recall = history.history['recall'][-1]
tp = history.history['true_positives'][-1]
tn = history.history['true_negatives'][-1]
fp = history.history['false_positives'][-1]
fn = history.history['false_negatives'][-1]
auc = history.history['auc'][-1]
f1 = 2*precision*recall / (precision+recall)<jupyter_output><empty_output><jupyter_text>Training stats<jupyter_code>print(f'precision: {precision}')
print(f'recall: {recall}')
print(f'f1 score: {f1}')
print(f'true pos: {tp}')
print(f'true neg: {tn}')
print(f'false pos: {fp}')
print(f'false neg: {fn}')
print(f'auc: {auc}')<jupyter_output>precision: 0.7135922312736511
recall: 0.4946322739124298
f1 score: 0.5842717913534946
true pos: 3087.0
true neg: 18568.0
false pos: 1239.0
false neg: 3154.0
auc: 0.8746693134307861<jupyter_text>Test stats<jupyter_code>results = logistic.evaluate(x_test, y_test, batch_size=32)
testf1 = 2 * results[1] * results[2] / (results[1] + results[2])
print(f'test f1 score: {testf1}')<jupyter_output>test f1 score: 0.572188456556238<jupyter_text>Improve through tuning elasticnet regularization<jupyter_code>logistic_reg = LogisticRegression(num_features=X.shape[1], reg = [0.0001, 0.0001])
logistic_reg.compile(optimizer='rmsprop',
                 loss='binary_crossentropy',
                 metrics=[tf.keras.metrics.Precision(), 
                          tf.keras.metrics.Recall(), 
                          tf.keras.metrics.TruePositives(), 
                          tf.keras.metrics.TrueNegatives(), 
                          tf.keras.metrics.FalsePositives(), 
                          tf.keras.metrics.FalseNegatives(), 
                          tf.keras.metrics.AUC()])

print(f'random output for initialization = {logistic_reg(X[0, :].reshape(1, -1))}')
logistic_reg.summary()
history_reg = logistic_reg.fit(x_train, y_train, epochs=5, batch_size=32)
history_reg.history.keys()
precision = history_reg.history['precision_1'][-1]
recall = history_reg.history['recall_1'][-1]
tp = history_reg.history['true_positives_1'][-1]
tn = history_reg.history['true_negatives_1'][-1]
fp = history_reg.history['false_positives_1'][-1]
fn = history_reg.history['false_negatives_1'][-1]
auc = history_reg.history['auc_1'][-1]<jupyter_output><empty_output><jupyter_text>Improved training stats<jupyter_code>print(f'precision: {precision}')
print(f'recall: {recall}')
print(f'f1 score: {f1}')
print(f'true pos: {tp}')
print(f'true neg: {tn}')
print(f'false pos: {fp}')
print(f'false neg: {fn}')
print(f'auc: {auc}')<jupyter_output>precision: 0.7128393054008484
recall: 0.5302035212516785
f1 score: 0.5842717913534946
true pos: 3309.0
true neg: 18474.0
false pos: 1333.0
false neg: 2932.0
auc: 0.8847647309303284<jupyter_text>Improved test stats<jupyter_code>results = logistic_reg.evaluate(x_test, y_test, batch_size=32)
testf1 = 2 * results[1] * results[2] / (results[1] + results[2])
print(f'test f1 score: {testf1}')<jupyter_output>test f1 score: 0.6130075325051406
{'fraction_non_alphanumeric': 0.13990451441348886, 'fraction_numerical': 0.03589774581495135, 'mean_word_length': 3.7226027397260273, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 86, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '549075', 'license': 'MIT', 'n_tokens_mistral': 6109, 'n_tokens_neox': 5649, 'n_words': 1307, 'path': 'assignments/PA-4/MA544_S21_PA4.ipynb', 'repo_name': 'jschmidtnj/ma544-final-project', 'size': '16547'}