Excuse the messiness – not tidied up or optimised yet, but just to see how it’d work.
This gets 55.9% accuracy relatively quickly with the Kaggle Facial Expression dataset, with 7 possible expressions, after converting them to feature vectors first. However it fails to convert about 29% of the original dataset – presumably because it wasn’t trained on small, grayscale images. So you’d have to mark it down on the challenge. The best model on Kaggle got about 70%, with the runners up around 50%, but you might be able to get this higher with more training.
import torchvision.datasets as dsets
import torchvision.transforms as transforms
import csv
import os
from PIL import Image, ImageDraw
import face_recognition
import torch
from torch import nn
import numpy as np
from torch.utils.data import Dataset, DataLoader
import math


def str2act(s):
    if s is 'none':
        return None
    elif s is 'hardtanh':
        return nn.Hardtanh()
    elif s is 'sigmoid':
        return nn.Sigmoid()
    elif s is 'relu6':
        return nn.ReLU6()
    elif s is 'tanh':
        return nn.Tanh()
    elif s is 'tanhshrink':
        return nn.Tanhshrink()
    elif s is 'hardshrink':
        return nn.Hardshrink()
    elif s is 'leakyrelu':
        return nn.LeakyReLU()
    elif s is 'softshrink':
        return nn.Softshrink()
    elif s is 'softsign':
        return nn.Softsign()
    elif s is 'relu':
        return nn.ReLU()
    elif s is 'prelu':
        return nn.PReLU()
    elif s is 'softplus':
        return nn.Softplus()
    elif s is 'elu':
        return nn.ELU()
    elif s is 'selu':
        return nn.SELU()
    else:
        raise ValueError("[!] Invalid activation function.")


class MLP(nn.Module):
    def __init__(self, num_layers, in_dim, hidden_dim, out_dim, activation='relu'):
        super().__init__()
        self.num_layers = num_layers
        self.in_dim = in_dim
        self.hidden_dim = hidden_dim
        self.out_dim = out_dim
        self.activation = str2act(activation)

        nonlin = True
        if self.activation is None:
            nonlin = False

        layers = []
        for i in range(num_layers - 1):
            layers.extend(
                self._layer(
                    hidden_dim if i > 0 else in_dim,
                    hidden_dim,
                    nonlin,
                )
            )
        layers.extend(self._layer(hidden_dim, out_dim, False))

        self.model = nn.Sequential(*layers)

    def _layer(self, in_dim, out_dim, activation=True):
        if activation:
            return [
                nn.Linear(in_dim, out_dim),
                self.activation,
            ]
        else:
            return [
                nn.Linear(in_dim, out_dim),
            ]

    def forward(self, x):
        out = self.model(x.float())
        return out



def _load_data(path='fer2013.csv', expect_labels=True):

    assert path.endswith('.csv')

    # If a previous call to this method has already converted
    # the data to numpy format, load the numpy directly
    X_path = path[:-4] + '.X.npy'
    Y_path = path[:-4] + '.Y.npy'
    if os.path.exists(X_path):
        X = np.load(X_path)
        if expect_labels:
            y = np.load(Y_path)
        else:
            y = None
        return X, y

    csv_file = open(path, 'r')
    reader = csv.reader(csv_file)

    # Discard header
    row = next(reader)

    y_list = []
    X_list = []
    counter = 0
    skip_counter = 0

    for i, row in enumerate(reader):
        counter +=1
        y_str, X_row_str = (row[0], row[1])
        y = int(y_str)
        
        X_row_strs = X_row_str.split(' ')
        X_row = [float(x) for x in X_row_strs]
        X_row = np.reshape(X_row, (48,48))
        image = np.zeros((48, 48, 3), dtype=np.uint8)
        
        image[:,:,0] = X_row
        image[:,:,1] = X_row
        image[:,:,2] = X_row
        print(i)
        face_landmarks_list = face_recognition.face_landmarks(image)
        X_row = image

        landmarks_array = []
        for face_landmarks in face_landmarks_list:
            for facial_feature in face_landmarks.keys():
                for item in face_landmarks[facial_feature]:
                    landmarks_array.append(np.round(item[0] / 48, 5))
                    landmarks_array.append(np.round(item[1] / 48, 5))

        if face_landmarks_list:
            X_list.append(landmarks_array)
            y_list.append(y)
        else:
            skip_counter +=1
        
    X = np.asarray(X_list)
    y = np.asarray(y_list)
    np.save(X_path, X)
    np.save(Y_path, y)
    print(skip_counter,'missed, out of', counter, ' - total:', (counter - skip_counter) / counter, '%')

    return X, y


class PrepareData(Dataset):
    def __init__(self, x, y):
        self.x = torch.from_numpy(x) if not torch.is_tensor(x) else x
        self.y = torch.from_numpy(y) if not torch.is_tensor(y) else y

    def __len__(self):
        return len(self.x)

    def __getitem__(self, idx):
        return self.x[idx], self.y[idx]


X, y = _load_data()

NUM_TRAINING_IMAGES = int(len(X) * 0.9)
print(NUM_TRAINING_IMAGES)

trainer_loader = PrepareData(x=X[:NUM_TRAINING_IMAGES], y=y[:NUM_TRAINING_IMAGES])
test_loader = PrepareData(x=X[NUM_TRAINING_IMAGES:], y=y[NUM_TRAINING_IMAGES:])

torch.manual_seed(1)

# Hyper Parameters
EPOCH = 1000
BATCH_SIZE = 16
LR = 0.0001

# Data Loader for easy mini-batch return in training
train_loader = DataLoader(dataset=trainer_loader, batch_size=BATCH_SIZE, shuffle=True)
test_loader = DataLoader(test_loader, batch_size=1, shuffle=True)

model = MLP(num_layers=5, in_dim=144, hidden_dim=256, out_dim=7)
print(model)

optimizer = torch.optim.Adam(model.parameters(), lr=LR)
# optimizer = torch.optim.SGD(model.parameters(), lr=LR, momentum=0.7)
loss_func = nn.CrossEntropyLoss()
average_loss = 2
correct_ratio = 0


for epoch in range(EPOCH):
    for step, (b_x, b_y) in enumerate(train_loader):
        # b_x = b_x.view(-1, 28, 28)        # reshape x to (batch, time_step, input_size)
        output = model(b_x)
        loss = loss_func(output, b_y)
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()

        if step % 10 == 0:
            average_loss = average_loss * 0.95 + loss.item() * 0.05
            print('Epoch:', epoch, '- step:', step, '- loss:', np.round(loss.item(),5), 
                '- average loss:', np.round(average_loss, 5), '- correct:', correct_ratio, '%')
        
    if epoch % 5 == 0:
        correct = 0
        incorrect = 0

        for data, target in test_loader:
            predicted_answer = np.argmax(model.forward(data).detach().numpy())
            right_answer = target[0]
            if int(right_answer) == int(predicted_answer):
                correct+=1
            else:
                incorrect+=1
        print('Test Correct: {}/{}'.format(correct, correct+incorrect))
        correct_ratio = np.round(100 * (correct / (correct+incorrect)), 1)