# Assuming to have: trainloader, testloader, and model with model.backbone() and model.classifer() methods

# Extract the training set features. Dimensions 2D matrix: training_size x num_features
features = torch.vstack([model.backbone(images) for images, label in trainloader])
norm = torch.linalg.norm(features, dim=-1, keepdim=True)
features = features / norm

# Compute the training data empirical correlation matrix inverse
x_t_x_inv = torch.linalg.inv(features.T @ features)

# Calculate the regret: Large regret means out-of-distribution sample
for images, labels in testloader:
    features = model.backbone(images)

    # Normalize
    norm = torch.linalg.norm(features, dim=-1, keepdim=True)
    features = features / norm
    
    # Get the probabilities of the normalized features
    probs = torch.softmax(model.classifier(features), dim=-1)

    # Calc projection
    x_proj = features @ x_t_x_inv @ features.T
    xt_g = x_proj / (1 + x_proj)

    # Equation 20
    n_classes = probs.shape[-1]
    nf = torch.sum(probs / (probs + (1 - probs) * (probs ** x_t_g)), dim=-1)
    regrets = torch.log(nf) / torch.log(torch.tensor(n_classes))