(base) liziyu@liziyudeMacBook-Pro ~ % ssh -R 7891:localhost:7890 zyli@lab
Welcome to Ubuntu 24.04.3 LTS (GNU/Linux 6.14.0-29-generic x86_64)

 * Documentation:  https://help.ubuntu.com
 * Management:     https://landscape.canonical.com
 * Support:        https://ubuntu.com/pro

Expanded Security Maintenance for Applications is not enabled.

124 updates can be applied immediately.
To see these additional updates run: apt list --upgradable

Enable ESM Apps to receive additional future security updates.
See https://ubuntu.com/esm or run: sudo pro status

*** System restart required ***
Last login: Fri Jan 23 15:16:55 2026 from 61.169.135.222
(base) zyli@lab:~$ 

# Load the FashionMNIST dataset from torchvision
train_data = torchvision.datasets.FashionMNIST(root='./data', train=True, download=True)

# Extract the image data and convert it to a numpy array of type float
images = train_data.data.numpy().astype(float)

# Extract the target labels as a numpy array
targets = train_data.targets.numpy()

# Create a dictionary mapping class indices to class names
class_dict = {i: class_name for i, class_name in enumerate(train_data.classes)}

# Map the target labels to their corresponding class names
labels = np.array([class_dict[t] for t in targets])

# Create a list of class names in order of their indices
class_names = [class_dict[i] for i in range(len(class_dict))]

# Get the total number of samples in the dataset
n = len(images)

# Ensure class_names is a list of class names (redundant but ensures consistency)
class_names = list(class_dict.values())

# Print dataset information for verification
print("Loaded FashionMNIST dataset with {} samples.".format(n))
print("Classes: {}".format(class_dict))
print("Image shape: {}".format(images[0].shape))  # Shape of a single image
print("Image dtype: {}".format(images[0].dtype))  # Data type of the image array
print("Image type: {}".format(type(images[0])))   # Type of the image object

def show_images(images, max_images=40, ncols=5, labels = None, reshape=False):
    """Visualize a subset of images from the dataset.
    Args:
        images (np.ndarray or list): Array of images to visualize [img,row,col].
        max_images (int): Maximum number of images to display.
        ncols (int): Number of columns in the grid.
        labels (np.ndarray, optional): Labels for the images, used for facet titles.
    Returns:
        plotly.graph_objects.Figure: A Plotly figure object containing the images.
    """
    if isinstance(images, list):
        images = np.stack(images)
    n = min(images.shape[0], max_images) # Number of images to show
    px_height = 220 # Height of each image in pixels
    if reshape:
        images = images.reshape(images.shape[0], 28, 28)
    fig = px.imshow(images[:n, :, :], color_continuous_scale='gray_r', 
                    facet_col = 0, facet_col_wrap=ncols,
                    height = px_height * int(np.ceil(n/ncols)))
    fig.update_layout(coloraxis_showscale=False)
    fig.update_xaxes(showticklabels=False, showgrid=False)
    fig.update_yaxes(showticklabels=False, showgrid=False)
    if labels is not None:
        # Extract the facet number and replace with the label.
        fig.for_each_annotation(lambda a: a.update(text=labels[int(a.text.split("=")[-1])]))
    return fig

# flatten features into 1D arrays
X_train = np.stack(train_df['image'].values)
y_train = train_df['label'].values
X_test = np.stack(test_df['image'].values)
y_test = test_df['label'].values

print(f"X_train shape: {X_train.shape}\t y_train shape: {y_train.shape}")
print(f"X_test shape: {X_test.shape}\t y_test shape: {y_test.shape}")

# TODO: Train the model using the scaled traning data and plott the loss curve (remeber to normalize your data!)
# NOTE: Your model must be named `model`
def flatten(images):
    return images.reshape(images.shape[0], -1)
image_scalar=StandardScaler()
image_scalar.fit(flatten(X_train))
X_train_sc=image_scalar.transform(flatten(X_train))
X_test_sc=image_scalar.transform(flatten(X_test))

if load_saved_models and os.path.exists('classifier.joblib'):

    model = joblib.load('classifier.joblib')

else:

    model=MLPClassifier(
    hidden_layer_sizes=(100, 50),
    max_iter=200, tol=1e-4, random_state=SEED)

    model.fit(X_train_sc,y_train)
if save_models:
    joblib.dump(model, 'classifier.joblib')

loss_df=pd.DataFrame({'epoch':np.arange(1,len(model.loss_curve_)+1),'loss':model.loss_curve_})
loss_df.plot(x='epoch', y='loss', title="Training Error")

# TODO: Calculate train and test accuracy per class 
# TODO: Use class_accuracy to create a grouped bar chart of class accuracy for train and test

train_class_acc=train_df.groupby('label')['correct'].mean().reset_index()
train_class_acc['split']='train'

test_class_acc=test_df.groupby('label')['correct'].mean().reset_index()
test_class_acc['split']='test'

all_df=pd.concat([train_class_acc,test_class_acc])

class_accuracy=all_df[['split','label','correct']]
# print(class_accuracy)

class_accuracy_pivot=class_accuracy.pivot(index='label',columns='split',values='correct')
print(class_accuracy_pivot)

class_accuracy_pivot.plot(kind='bar',text_auto=True)

# Initialize confusion matrix with zeros
conf_matrix = np.zeros((len(class_names), len(class_names)), dtype=int)
class_to_idx = {class_name: idx for idx, class_name in enumerate(class_names)}

# Fill the confusion matrix by counting predictions and plot it as a heatmap
y_true=test_df['label'].values
y_pred=test_df['predicted_label'].values

true_indices = np.array([class_to_idx[label] for label in y_true])
pred_indices = np.array([class_to_idx[label] for label in y_pred])

np.add.at(conf_matrix, (true_indices, pred_indices), 1)

fig=px.imshow(conf_matrix,
                labels=dict(x="Predicted Label", y="True Label", color="Count"),
                x=class_names,
                y=class_names,
                )

fig.show()

# Calculate accuracy from confusion matrix
accuracy_from_matrix = np.trace(conf_matrix)/np.sum(conf_matrix)
print(f"\nAccuracy calculated from confusion matrix: {accuracy_from_matrix:.3f}")

# Calculate per-class metrics from confusion matrix
per_class_metrics = []
print("\nPer-class metrics from confusion matrix:")
for i, class_name in enumerate(class_names):
    true_positives = conf_matrix[i][i]
    false_positives = np.sum(conf_matrix[:,i])-conf_matrix[i][i]
    false_negatives = np.sum(conf_matrix[i,:])-conf_matrix[i][i]
    precision = true_positives / (true_positives + false_positives) if (true_positives + false_positives) > 0 else 0.0
    recall = true_positives / (true_positives + false_negatives) if (true_positives + false_negatives) > 0 else 0.0
    per_class_metrics.append({
        'class': class_name,
        'precision': precision,
        'recall': recall
    })
    
pd.DataFrame(per_class_metrics)

def create_vertical_flip_matrix(height=28, width=28):
    """
    Returns a (height*width, height*width) matrix that vertically flips an image
    when applied to its flattened vector. Values are 0 or 1.
    """
    N = height * width  # Total number of pixels in the image
    T = np.zeros((N, N), dtype=int)  # Initialize the transformation matrix with zeros
    for i in range(height):  # Loop over each row
        for j in range(width):  # Loop over each column
            orig_idx = i * width + j  # Compute the flattened index for the original pixel
            flipped_i = height - 1 - i  # Compute the row index after vertical flip
            flipped_idx = flipped_i * width + j  # Compute the flattened index for the flipped pixel
            # Set the corresponding entry in the transformation matrix to 1
            # This means the pixel at (i, j) moves to (flipped_i, j)
            T[flipped_idx, orig_idx] = 1
    return T

def vertical_flip(image):
    T_flip = create_vertical_flip_matrix()
    return apply_transformation(image, T_flip)

def create_horizontal_flip_matrix(height=28, width=28):
    """
    Returns a (height*width, height*width) matrix that horizontally flips an image
    when applied to its flattened vector. Values are 0 or 1.
    """
    N = height * width
    T = np.zeros((N, N), dtype=int)
    for i in range(height):
        for j in range(width):
            orig_idx = i * width + j
            flipped_j = width - 1 - j
            flipped_idx = i * width + flipped_j
            T[flipped_idx, orig_idx] = 1
    return T

    
def horizontal_flip(image):
    T_flip = create_horizontal_flip_matrix()
    return apply_transformation(image, T_flip)

flipped_image = horizontal_flip(test_image)

show_images(np.stack([test_image, flipped_image]), labels=['Original', 'Horizontal Flipped'], reshape=True)

def create_shift_matrix(dx, dy, height=28, width=28):
    """
    Create a transformation matrix for shifting an image by dx pixels horizontally and dy pixels vertically.

    Args:
        dx (int): Number of pixels to shift horizontally.
        dy (int): Number of pixels to shift vertically.
        height (int): Height of the image.
        width (int): Width of the image.

    Returns:
        np.ndarray: A (height*width, height*width) transformation matrix for shifting.
    """
    N = height * width
    T = np.zeros((N, N))
    yy, xx = np.meshgrid(np.arange(height), np.arange(width), indexing='ij')
    yy_flat = yy.ravel()
    xx_flat = xx.ravel()

    y_shift = yy_flat + dy
    x_shift = xx_flat + dx

    valid = (
        (y_shift >= 0) & (y_shift < height) &
        (x_shift >= 0) & (x_shift < width)
    )

    src = yy_flat[valid] * width + xx_flat[valid]
    dst = y_shift[valid] * width + x_shift[valid]

    T[dst, src] = 1.0

    return T


def shift_image(image, dx, dy):
    """
    Shift an image by dx pixels horizontally and dy pixels vertically.

    Args:
        image (np.ndarray): Flattened image array of shape (height*width,).
        dx (int): Number of pixels to shift horizontally.
        dy (int): Number of pixels to shift vertically.

    Returns:
        np.ndarray: Shifted image as a flattened array.
    """
    T = create_shift_matrix(dx, dy)
    return apply_transformation(image, T)

shifted_right_image = shift_image(test_image, 5, 0)
shifted_left_image = shift_image(test_image, -5, 0)
shifted_up_image = shift_image(test_image, 0, -5)
shifted_down_image = shift_image(test_image, 0, 5)

all_images = np.stack([test_image, shifted_up_image, shifted_down_image, shifted_right_image, shifted_left_image])
plot_labels = ['Original', 'Shifted Up', 'Shifted Down', 'Shifted Right', 'Shifted Left']
show_images(all_images, labels=plot_labels, reshape=True)

def create_blur_matrix(kernel_size=3, height=28, width=28):
    """
    Create a transformation matrix T that applies a uniform mean blur using a centered, odd-sized square sliding window.

    For each output pixel (i, j):
      1) Place a `kernel_size × kernel_size` window centered at (i, j).
      2) If the window is outside the image, then it will have fewer neighbors (only average the pixels that exist)

    Args:
        kernel_size (int): Size of the square kernel (must be odd).
        height (int): Height of the image.
        width (int): Width of the image.

    Returns:
        np.ndarray: A (height*width, height*width) transformation matrix for blurring.
    """
    N = height * width
    T = np.zeros((N, N))
    pad = kernel_size // 2
    kernel_area = kernel_size * kernel_size
    # Each pixel contributes equally to a kernel's neighborhood
    for i in range (height):
        for j in range (width):
            col_start=max(0,j-pad)
            row_start=max(0,i-pad)
            col_end=min(width-1,j+pad)
            row_end=min(height-1,i+pad)

            valid_count=(row_end-row_start+1)*(col_end-col_start+1)
            arg_value=1.0/valid_count

            dst_idx=i*width+j
            for r in range(row_start,row_end+1):
                for c in range(col_start,col_end+1):
                    src_idx=r*width+c
                    T[dst_idx,src_idx]=arg_value
    return T

def blur_image(image, kernel_size=3):
    """
    Apply a blur transformation to a flattened image array or a batch of flattened images.

    Args:
        image (np.ndarray): Flattened image array of shape (height*width,) or batch of images (N, height*width).
        kernel_size (int): Size of the square kernel to use for blurring.

    Returns:
        np.ndarray: Blurred image(s) as a flattened array or batch of arrays.
    """
    T = create_blur_matrix(kernel_size)
    return apply_transformation(image, T)

blurred_1x1 = blur_image(test_image, kernel_size=1)
blurred_3x3 = blur_image(test_image, kernel_size=3)
blurred_5x5 = blur_image(test_image, kernel_size=5)

blurred_images = [test_image, blurred_1x1, blurred_3x3, blurred_5x5]
blurred_labels = ['Original', 'Blur 1x1', 'Blur 3x3', 'Blur 5x5']

show_images(blurred_images, labels=blurred_labels, reshape=True)

def create_rotation_matrix(theta, height=28, width=28):
    """
    Create a transformation matrix for rotating an image by theta degrees.

    Args:
        theta (float): Angle of rotation in degrees.
        height (int): Height of the image.
        width (int): Width of the image.

    Returns:
        np.ndarray: A (height*width, height*width) transformation matrix for rotating.
    """
    # Convert theta from degrees to radians
    theta = np.deg2rad(theta)
    N = height * width
    T = np.zeros((N, N))
    center_i = (height - 1) / 2.0
    center_j = (width - 1) / 2.0

    cos_theta = np.cos(theta)
    sin_theta = np.sin(theta)
    
    for src_i in range(height):
        for src_j in range(width):

            src_i_centered = src_i - center_i
            src_j_centered = src_j - center_j
            
            # 旋转（逆时针旋转theta）
            dst_i_centered = src_i_centered * cos_theta - src_j_centered * sin_theta
            dst_j_centered = src_i_centered * sin_theta + src_j_centered * cos_theta
            
            # 平移回去
            dst_i = dst_i_centered + center_i
            dst_j = dst_j_centered + center_j
            
            # 取整并检查边界
            dst_i_int = int(np.round(dst_i))
            dst_j_int = int(np.round(dst_j))
            
            if 0 <= dst_i_int < height and 0 <= dst_j_int < width:
                src_idx = src_i * width + src_j
                dst_idx = dst_i_int * width + dst_j_int
                T[dst_idx, src_idx] = 1.0

    return T


def rotate_image(image, theta):
    """
    Apply a rotation transformation to a flattened image array or a batch of flattened images.

    Args:
        image (np.ndarray): Flattened image array of shape (height*width,) or batch of images (N, height*width).
        theta (float): Angle of rotation in degrees.

    Returns:
        np.ndarray: Rotated image(s) as a flattened array or batch of arrays.
    """
    T = create_rotation_matrix(theta)
    return apply_transformation(image, T)

# rotate with matrix
rotated_45 = rotate_image(test_image, 45) 
rotated_90 = rotate_image(test_image, 90)
rotated_200 = rotate_image(test_image, 200)
rotated_270 = rotate_image(test_image, 270)

# visualize original and 4 augmentations in plotly image grid
all_images = np.stack([test_image, rotated_45, rotated_90, rotated_200, rotated_270])
plot_labels = ['Original', 'Rotated (45°)', 'Rotated (90°)', 'Rotated (200°)', 'Rotated (270°)']
show_images(all_images, labels=plot_labels, reshape=True)

def create_bilinear_rotation_matrix(theta, height=28, width=28):
    """
    Create a (height*width, height*width) matrix that applies bilinear interpolation
    for rotating a flattened image by theta degrees.
    Each row of the matrix gives the weights for the input pixels that contribute to each output pixel.

    Args:
        theta (float): Angle of rotation in degrees.
        height (int): Height of the image.
        width (int): Width of the image.

    Returns:
        np.ndarray: A (height*width, height*width) transformation matrix for rotating.
    """
    theta = np.deg2rad(theta)
    N = height * width
    T = np.zeros((N, N))
    center_i = (height - 1) / 2.0
    center_j = (width - 1) / 2.0
    cos_theta = np.cos(-theta)
    sin_theta = np.sin(-theta)

    for i in range(height):  # Loop over rows of the output image
        for j in range(width):  # Loop over columns of the output image
            output_idx = i * width + j

            # Translate output pixel to be relative to the center
            centered_j = j - center_j
            centered_i = i - center_i

            # Apply the inverse rotation to find the original pixel coordinates
            original_j_rotated = centered_j * cos_theta - centered_i * sin_theta
            original_i_rotated = centered_j * sin_theta + centered_i * cos_theta

            # Translate back to original image coordinates
            original_j = original_j_rotated + center_j
            original_i = original_i_rotated + center_i

            # Bilinear interpolation
            # Find the 4 nearest neighbor pixels
            j_floor = int(np.floor(original_j))
            i_floor = int(np.floor(original_i))
            j_ceil = int(np.ceil(original_j))
            i_ceil = int(np.ceil(original_i))

            # Calculate the fractional distances
            dj = original_j - j_floor
            di = original_i - i_floor

            # Define the weights
            w_tl = (1 - dj) * (1 - di)  # Top-left
            w_tr = dj * (1 - di)      # Top-right
            w_bl = (1 - dj) * di      # Bottom-left
            w_br = dj * di          # Bottom-right

            # Get the indices of the 4 neighbors, handling boundaries by setting weight to 0
            neighbors = []
            weights = []

            # Top-left neighbor
            if 0 <= i_floor < height and 0 <= j_floor < width:
                neighbors.append(i_floor * width + j_floor)
                weights.append(w_tl)

            # Top-right neighbor
            if 0 <= i_floor < height and 0 <= j_ceil < width:
                neighbors.append(i_floor * width + j_ceil)
                weights.append(w_tr)

            # Bottom-left neighbor
            if 0 <= i_ceil < height and 0 <= j_floor < width:
                neighbors.append(i_ceil * width + j_floor)
                weights.append(w_bl)

            # Bottom-right neighbor
            if 0 <= i_ceil < height and 0 <= j_ceil < width:
                neighbors.append(i_ceil * width + j_ceil)
                weights.append(w_br)

            # Assign the weighted sum to the output pixel in the transformation matrix
            for neighbor_idx, weight in zip(neighbors, weights):
                T[output_idx, neighbor_idx] = weight
    return T

# Top-left neighbor
if 0 <= i_floor < height and 0 <= j_floor < width:
    neighbors.append(i_floor * width + j_floor)
    weights.append(w_tl)

# Top-right neighbor
if 0 <= i_floor < height and 0 <= j_ceil < width:
    neighbors.append(i_floor * width + j_ceil)
    weights.append(w_tr)

# Bottom-left neighbor
if 0 <= i_ceil < height and 0 <= j_floor < width:
    neighbors.append(i_ceil * width + j_floor)
    weights.append(w_bl)

# Bottom-right neighbor
if 0 <= i_ceil < height and 0 <= j_ceil < width:
    neighbors.append(i_ceil * width + j_ceil)
    weights.append(w_br)

# Assign the weighted sum to the output pixel in the transformation matrix
for neighbor_idx, weight in zip(neighbors, weights):
    T[output_idx, neighbor_idx] = weight

def rotate_image_bilinear(image, theta):
    """
    Rotate an image using bilinear interpolation.

    Args:
        image (np.ndarray): Flattened image array of shape (height*width,) or batch of images (N, height*width).
        theta (float): Angle of rotation in degrees.

    Returns:
        np.ndarray: Rotated image as a flattened array.
    """
    T = create_bilinear_rotation_matrix(theta)
    return apply_transformation(image, T)
    
# rotate with matrix
rotated = rotate_image(test_image, 45)
rotated_interpolated = rotate_image_bilinear(test_image, 45)

all_images = np.stack([test_image, rotated, rotated_interpolated])
plot_labels = ['Original',  'Rotated 45°', 'Rotated 45° (Bilinear)']
show_images(all_images, labels=plot_labels, reshape=True)

# Test augmentation functions on a few examples
test_images = np.stack(test_df['image'])
test_labels = test_df['label']

shift_inputs = [(5, 0), (-5, 0), (0, 5), (0, -5)]
rotate_inputs = [45, 90, 200]
blur_inputs = [3, 5]
rotate_blur_inputs = [(45, 3), (90, 5)]
shift_rotate_blur_inputs = [((5, 0), 45, 3), ((-5, 0), 90, 5)]

augmented_data = []
# Randomly sample 100 datapoints from test_images
sample_idx = np.random.choice(len(test_images), 100, replace=False)
test_images_sample = test_images[sample_idx]
test_labels_sample = np.array(test_labels)[sample_idx]

# TODO: Apply the augmentation functions we just created (shift, blur, rotate w/ bilinear, rotate then blur, shift then rotate then blur) to every image from test_images_sample
# use the inputs defined above to apply the augmentations
# Save the augmented images in a new DataFrame aug_df
def add_augmented_sample(img, label, orig_idx, aug_name, aug_type, transform_fn):
    augmented_data.append({
        'original_idx': orig_idx,
        'augmentation': aug_name,
        'image': transform_fn(img).copy(),
        'label': label,
        'type': aug_type
    })


for orig_idx, img, label in zip(sample_idx, test_images_sample, test_labels_sample):
    for dx, dy in shift_inputs:
        add_augmented_sample(
            img, label, orig_idx,f'shift_{dx}_{dy}', 'shift',
            lambda x, dx=dx, dy=dy: shift_image(x, dx, dy)
        )
    for theta in rotate_inputs:
        add_augmented_sample(
            img, label, orig_idx,f'rotate_bilinear_{theta}', 'rotate',
            lambda x, theta=theta: rotate_image_bilinear(x, theta)
        )
    for k in blur_inputs:
        add_augmented_sample(
            img, label, orig_idx,f'blur_{k}x{k}', 'blur',
            lambda x, k=k: blur_image(x, kernel_size=k)
        )
    for theta, k in rotate_blur_inputs:
        add_augmented_sample(
            img, label, orig_idx,f'rotate_{theta}_blur_{k}', 'rotate_blur',
            lambda x, theta=theta, k=k: rotate_then_blur(x, theta, k)
        )
    for (shift_pair, theta, k) in shift_rotate_blur_inputs:
        dx, dy = shift_pair
        add_augmented_sample(
            img, label, orig_idx,f'shift_{dx}_{dy}_rotate_{theta}_blur_{k}', 'shift_rotate_blur',
            lambda x, dx=dx, dy=dy, theta=theta, k=k: shift_then_rotate_then_blur(x, dx, dy, theta, k)
        )


aug_df = pd.DataFrame(augmented_data)

# TODO: Select an image and visualize it with all the augmentations applied to it
example_image = test_images_sample[0]

example_variants = [('original', example_image)]
for dx, dy in shift_inputs:
    example_variants.append((f'shift_{dx}_{dy}', shift_image(example_image, dx, dy)))
for theta in rotate_inputs:
    example_variants.append((f'rot_bilin_{theta}', rotate_image_bilinear(example_image, theta)))
for k in blur_inputs:
    example_variants.append((f'blur_{k}x{k}', blur_image(example_image, kernel_size=k)))
for theta, k in rotate_blur_inputs:
    example_variants.append((f'rot_{theta}_blur_{k}', rotate_then_blur(example_image, theta, k)))
for (shift_pair, theta, k) in shift_rotate_blur_inputs:
    dx, dy = shift_pair
    example_variants.append((f'shift_{dx}_{dy}_rot_{theta}_blur_{k}',
                             shift_then_rotate_then_blur(example_image, dx, dy, theta, k)))

example_imgs = np.stack([img for _, img in example_variants])
example_labels = [name for name, _ in example_variants]

fig = show_images(example_imgs, ncols=4, labels=example_labels, reshape=True)
fig.show()

from sklearn.metrics import accuracy_score


def evaluate_augmented_data(aug_df, clf):
    records = []

    for aug_name, group in aug_df.groupby("augmentation"):
        # 将图像堆叠成 3D/4D 数组，再 reshape 成模型需要的形状
        images = np.stack(group["image"].to_numpy())
        n_samples = images.shape[0]
        X = images.reshape(n_samples, -1)

        # 若模型需要标准化，这里先拟合再变换（或对训练集拟合后在此仅 transform）
        scaler = StandardScaler()
        X_scaled = scaler.fit_transform(X)

        y_true = group["label"].to_numpy()
        y_pred = clf.predict(X_scaled)
        acc = accuracy_score(y_true, y_pred)

        records.append({
            "augmentation": aug_name,
            "accuracy": acc,
            "type": group["type"].iloc[0]
        })

    return pd.DataFrame(records)


aug_performance = evaluate_augmented_data(aug_df, model)
aug_performance = aug_performance.sort_values("accuracy", ascending=False).reset_index(drop=True)
aug_performance
# Visualize performance: sort by accuracy, color by augmentation type (blur, rotate, shift, none)

fig = px.bar(
    aug_performance,
    x="augmentation",
    y="accuracy",
    color="type",
    title="Classifier Accuracy on Augmented Data"
)
fig.update_layout(xaxis_title="Augmentation", yaxis_title="Accuracy")
fig.show()
<img src="/images/189_hw1/189_hw1_20.png" alt="label_plot" />