# Suspicious Cakes: Official Solution

*English translation by SOTA – AI Community of the Russian original. Organisers who would like this translation removed can email sota.ai.community@gmail.com.*

All-Russian School Olympiad in Informatics 2025–2026, Final Stage<br>
“Artificial Intelligence” profile, Tour 2, Moscow, 25 March 2026<br>
Task E

*Translator's note: the official solution is the authors' Jupyter notebook exported to PDF. Code cells are marked “In [n]”, their results “Out [n]”, and printed output follows the cell. Lines that the PDF wrapped are joined again. Comments, docstrings and messages are translated; the code itself is unchanged. Figures are described in words; see the original PDF for the images.*

The solution consists of the following stages:

1. We load the pre-trained CNN and obtain embeddings for the images.
2. For each image, we compute the outlier metrics `robust_mahalanobis`, `global_knn`, `class_knn`.
3. We aggregate the metrics with the weights (2.0, 1.5, 0.5).
4. We sort the images by the value of the metric, select the top-K and save them to `submission_author.csv`.

In [1]:

```python
import json
import numpy as np
import pandas as pd
import torch
import torch.nn as nn
from torch.utils.data import DataLoader, TensorDataset

ARTIFACTS_DIR = "."
TEST_PACKAGE = f"{ARTIFACTS_DIR}/public_test_package.npz"
TEST_META = f"{ARTIFACTS_DIR}/public_test_meta.json"
WEIGHTS_PATH = f"{ARTIFACTS_DIR}/model_weights.pt"
SUBMISSION_PATH = "submission.csv"

DEVICE = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print("Device:", DEVICE)
```

```text
Device: cpu
```

The code of the network and of obtaining the image embeddings:

In [2]:

```python
class SmallCNN(nn.Module):
    def __init__(self, n_classes=10):
        super().__init__()
        self.features = nn.Sequential(
            nn.Conv2d(1, 16, kernel_size=3, padding=1),
            nn.ReLU(inplace=True),
            nn.MaxPool2d(2),
            nn.Conv2d(16, 32, kernel_size=3, padding=1),
            nn.ReLU(inplace=True),
            nn.MaxPool2d(2),
            nn.Conv2d(32, 64, kernel_size=3, padding=1),
            nn.ReLU(inplace=True),
        )
        self.fc1 = nn.Linear(64 * 8 * 8, 64)
        self.act = nn.ReLU(inplace=True)
        self.fc2 = nn.Linear(64, n_classes)

    def forward(self, x):
        z = self.features(x)
        z = z.flatten(1)
        h = self.act(self.fc1(z))
        logits = self.fc2(h)
        return logits, h

def infer_embeddings_and_preds(model, x_np, batch_size=256):
    ds = TensorDataset(torch.from_numpy(x_np).float())
    loader = DataLoader(ds, batch_size=batch_size, shuffle=False)
    all_emb, all_pred = [], []
    model.eval()
    with torch.no_grad():
        for (xb,) in loader:
            xb = xb.to(DEVICE)
            logits, emb = model(xb)
            all_emb.append(emb.cpu().numpy())
            all_pred.append(logits.argmax(dim=1).cpu().numpy())
    emb = np.concatenate(all_emb, axis=0).astype(np.float32)
    pred = np.concatenate(all_pred, axis=0).astype(np.int64)
    return emb, pred

state = torch.load(WEIGHTS_PATH, map_location=DEVICE)
n_classes = int(state["fc2.weight"].shape[0]) if "fc2.weight" in state else 10
model = SmallCNN(n_classes=n_classes)
model.load_state_dict(state)
model = model.to(DEVICE).eval()
print("Model loaded")
```

```text
Model loaded
```

In [3]:

```python
test = np.load(TEST_PACKAGE)
x_test = test["images"].astype(np.float32)

K = 1000

emb_test, pred_test = infer_embeddings_and_preds(model, x_test)

n_test = len(x_test)
print(f"n_test={n_test}, K={K}, emb_dim={emb_test.shape[1]}")
```

```text
n_test=10000, K=1000, emb_dim=64
```

### 1.12.9 The idea of the solution

We will look for outliers in the embedding space of the neural network.

We assume that, after passing through the trained model, objects of the same class form compact clusters in the feature space. Then outliers can be sought as points that are located "atypically" relative to the other points.

For this, several types of metrics are used:

1. **Global kNN score.**
   For each object, the mean distance to its nearest neighbours among all objects of the dataset is computed.
   If an object is isolated, its score will be large.
2. **Class kNN score.**
   Similar to the previous item, but the neighbours are sought only among objects of the same predicted class.
   This makes it possible to find points that look unusual specifically within their own class, even if globally they are not very far from the others.
3. **Robust Mahalanobis score within the class.**
   For each predicted class, the mean $\mu$ and the covariance matrix $\Sigma$ of the embeddings are estimated, after which, for each object, the Mahalanobis distance to the centre of its class is computed. The Mahalanobis distance is computed as follows:

   $$(x_i - \mu)^T \Sigma^{-1} (x_i - \mu)$$

   The idea of using the covariance is that Mahalanobis looks at a point relative to the shape of the class distribution. Sometimes an outlier point can be caught by the fact that it deviates not in length but in a "wrong direction" in the space.
   So that outliers do not distort the initial sample estimates of $\mu$ and $\Sigma$, we use a robust scheme: at each iteration we temporarily discard the farthest points and estimate the parameters only from the remaining "core" of the class.

Thus, the solution relies on the following hypothesis:

**outliers are objects that, in the feature space, are either far from the typical distribution of their class or fit poorly into the local structure of their neighbours.**

**Note**: without using the Mahalanobis distance, one can get up to 80% of the points for the task.

The metric functions are implemented below:

In [4]:

```python
import numpy as np

def score_robust_mahalanobis(emb, pred, lam=0.15, trim_frac=0.10, iters=2):
    """
    Computes a robust Mahalanobis score within each predicted class.

    Idea:
    - for each class, we estimate the "normal" distribution of the embeddings;
    - then, for each object, we compute the Mahalanobis distance to the centre of the class;
    - so that outliers do not spoil the estimates of the mean and the covariance,
      we discard the farthest points several times and recompute the statistics.

    Parameters:
    - emb: matrix of embeddings, shape [N, d]
    - pred: predicted classes, shape [N]
    - lam: covariance regularisation coefficient
    - trim_frac: fraction of the farthest points temporarily excluded during trimming
    - iters: number of iterations of the robust recomputation
    """
    d = emb.shape[1]  # dimension of the feature space
    out = np.zeros(len(emb), dtype=np.float64)  # final score for all objects
    eye = np.eye(d, dtype=np.float64)  # identity matrix for regularisation

    # Process each class separately
    for c in np.unique(pred):
        m = pred == c          # mask of the objects of class c
        h = emb[m].astype(np.float64)  # embeddings of this class only
        n = len(h)

        # If there are too few objects, a robust estimate is unreliable.
        # Then we simply compute the ordinary Mahalanobis distance.
        if n < 8:
            mu = h.mean(axis=0, keepdims=True)   # class centre
            cen = h - mu                         # centred embeddings

            # Covariance + diagonal regularisation for stability
            cov = (cen.T @ cen) / max(1, n - 1) + lam * eye

            # Pseudo-inverse instead of the ordinary inverse:
            # more stable if the covariance is ill-conditioned
            inv = np.linalg.pinv(cov)

            # Squared Mahalanobis distance for each point:
            # (x - mu)^T inv (x - mu)
            out[m] = np.einsum("bi,ij,bj->b", cen, inv, cen)
            continue

        # Initially, we assume that all points of the class are kept
        keep = np.ones(n, dtype=bool)

        # Initial rough estimate of the mean and of the inverse covariance
        mu = h.mean(axis=0, keepdims=True)
        inv = np.linalg.pinv(np.cov(h.T) + lam * eye)

        # Several iterations of robust trimming
        for _ in range(max(1, iters)):
            # Take only the current "reliable" points
            hh = h[keep]

            # Recompute the centre from the remaining points
            mu = hh.mean(axis=0, keepdims=True)

            # Centre the remaining points
            cen_hh = hh - mu

            # Recompute the covariance on the cleaned subset
            cov = (cen_hh.T @ cen_hh) / max(1, len(hh) - 1) + lam * eye
            inv = np.linalg.pinv(cov)

            # Compute the Mahalanobis distances, now for all points of the class
            cen_all = h - mu
            dist_all = np.einsum("bi,ij,bj->b", cen_all, inv, cen_all)

            # Threshold: keep the (1 - trim_frac) fraction of the nearest points
            thr = float(np.quantile(dist_all, 1.0 - trim_frac))
            keep = dist_all <= thr

            # Protection against too few points remaining:
            # in that case, relax the trimming
            if keep.sum() < max(5, int(0.5 * n)):
                keep = dist_all <= float(np.quantile(dist_all, 0.7))

        # After the final estimate, compute the final Mahalanobis score
        # for all points of this class
        cen = h - mu
        out[m] = np.einsum("bi,ij,bj->b", cen, inv, cen)

    return out


def score_global_knn(emb, k=10):
    """
    Computes the global kNN score:
    the mean distance to the k nearest neighbours over the whole dataset.

    If a point is isolated from the others,
    its score will be large.
    """
    n = emb.shape[0]

    # If there are too few points, there is no meaningful score
    if n <= 2:
        return np.zeros(n, dtype=np.float64)

    # We cannot take more neighbours than there are other points
    kk = min(k, n - 1)

    x = emb.astype(np.float64)

    # Squared norms of all points
    x2 = np.sum(x * x, axis=1, keepdims=True)

    # Matrix of squared pairwise Euclidean distances:
    # ||xi - xj||^2 = ||xi||^2 + ||xj||^2 - 2 <xi, xj>
    d2 = x2 + x2.T - 2.0 * (x @ x.T)

    # Remove the distance from each point to itself
    np.fill_diagonal(d2, np.inf)

    # Take the kk smallest distances in each row
    knn = np.partition(d2, kk - 1, axis=1)[:, :kk]

    # Return the mean ordinary distance to the nearest neighbours
    return np.mean(np.sqrt(np.maximum(knn, 0.0)), axis=1)


def score_class_knn(emb, pred, k=8):
    """
    Computes the class-wise kNN score:
    the mean distance to the k nearest neighbours within the predicted class.

    This helps to find objects that are atypical
    specifically for their own class.
    """
    out = np.zeros(len(emb), dtype=np.float64)

    # Compute the score separately within each class
    for c in np.unique(pred):
        m = pred == c
        h = emb[m].astype(np.float64)
        n = len(h)

        # If the class has too few points, the score is taken to be zero
        if n <= 2:
            out[m] = 0.0
            continue

        kk = min(k, n - 1)

        # Squared norms of the points within the class
        h2 = np.sum(h * h, axis=1, keepdims=True)

        # Matrix of squared pairwise distances within the class
        d2 = h2 + h2.T - 2.0 * (h @ h.T)

        # Exclude the distance to the point itself
        np.fill_diagonal(d2, np.inf)

        # Find the kk nearest neighbours
        knn = np.partition(d2, kk - 1, axis=1)[:, :kk]

        # Mean distance to the nearest neighbours within the class
        out[m] = np.mean(np.sqrt(np.maximum(knn, 0.0)), axis=1)

    return out

def rank_score(x):
    """
    Converts an array of values into ranks.

    The smallest element gets rank 0,
    the next one gets 1, and so on.
    This is convenient if several scores of different scales
    have to be combined later.
    """
    # Indices of the elements in ascending order of the values of x
    order = np.argsort(x)

    # Array for the ranks
    r = np.empty_like(order, dtype=np.float64)

    # For each element, write its position in the sorted order
    r[order] = np.arange(len(x), dtype=np.float64)
    return r
```

In [5]:

```python
# Build the outlier score on the test embeddings from an ensemble of the ranks
# of three metrics: robust Mahalanobis + global kNN + class kNN.

# Robust Mahalanobis distance within the predicted class:
# shows how atypical an object is relative to the distribution of its class.
s_robust = score_robust_mahalanobis(emb_test, pred_test, lam=0.15, trim_frac=0.10, iters=2)

# Global kNN score:
# the mean distance to the nearest neighbours in the whole set.
# Large values correspond to more isolated points.
s_gknn = score_global_knn(emb_test, k=10)

# Class-wise kNN score:
# the mean distance to the nearest neighbours only within the object's own predicted class.
# Helps to find points that look strange specifically within the class.
s_cknn = score_class_knn(emb_test, pred_test, k=8)

# Combine the three scores into one final suspiciousness indicator.
# Before combining, each score is converted into ranks,
# so that the different value scales do not interfere with each other.
# Here robust Mahalanobis has the largest weight,
# global kNN a medium one, and class kNN a smaller additional contribution.
final_score = (
    2.0 * rank_score(s_robust)
    + 1.5 * rank_score(s_gknn)
    + 0.5 * rank_score(s_cknn)
)

# Select the K objects with the largest final score
# as the most likely outliers.
topk = np.argsort(-final_score)[:K]

# Form the binary answer vector:
# 1 means the object is considered an outlier, 0 an ordinary object.
is_outlier = np.zeros(n_test, dtype=np.int64)
is_outlier[topk] = 1

# Assemble the file for submission:
# for each id, give the prediction is_outlier.
submission = pd.DataFrame({
    "id": np.arange(n_test, dtype=np.int64),
    "is_outlier": is_outlier,
})

# Save the submission to CSV without the DataFrame index.
submission.to_csv(SUBMISSION_PATH, index=False)

# Print service information:
# where the file was saved and how many objects are marked as outliers.
print("Saved:", SUBMISSION_PATH)
print("Predicted outliers:", int(is_outlier.sum()))

# Show the first rows of the submission table.
submission.head()
```

```text
Saved: submission.csv
Predicted outliers: 1000
```

Out [5]:

```text
   id  is_outlier
0   0           0
1   1           0
2   2           0
3   3           0
4   4           0
```

Let us check the solution: the cell below computes the score on the public and private parts of the dataset.

In [6]:

```python
#!/usr/bin/env python3
"""
Validate submission.csv against y_test.csv and compute hits@k metrics.

Expected columns:
  y_test.csv: id, is_outlier_true, Usage
  submission.csv: id, is_outlier
"""

from __future__ import annotations

import argparse
import json
from pathlib import Path
from typing import Dict, List

import pandas as pd

BASELINE_SOLUTION_SCORE_PUBLIC=113
BASELINE_SOLUTION_SCORE_PRIVATE=95
AUTHOR_SOLUTION_SCORE_PUBLIC=338
AUTHOR_SOLUTION_SCORE_PRIVATE=308


def _strip_jupyter_kernel_args(unknown: List[str]) -> List[str]:
    cleaned: List[str] = []
    i = 0
    while i < len(unknown):
        tok = unknown[i]
        if tok == "-f" and i + 1 < len(unknown) and unknown[i + 1].endswith(".json"):
            i += 2
            continue
        cleaned.append(tok)
        i += 1
    return cleaned


def _validate_binary_column(col: pd.Series, name: str) -> pd.Series:
    numeric = pd.to_numeric(col, errors="coerce")
    bad = numeric.isna() | (~numeric.isin([0, 1]))
    if bad.any():
        first_bad_idx = int(col.index[bad][0])
        raise ValueError(
            f"Column '{name}' must contain only 0/1. "
            f"First invalid row: {first_bad_idx}"
        )
    return numeric.astype(int)


def _clip_score_0_100(value: float) -> float:
    return max(0.0, min(100.0, float(value)))


def _compute_hits_metrics(
    df: pd.DataFrame, include_score_0_100: bool = True
) -> Dict[str, float | int | bool]:
    k_true = int(df["is_outlier_true"].sum())
    k_pred = int(df["is_outlier"].sum())
    hits_at_k = int(((df["is_outlier_true"] == 1) & (df["is_outlier"] == 1)).sum())
    if AUTHOR_SOLUTION_SCORE <= BASELINE_SOLUTION_SCORE:
        normalized_score_0_100 = 0.0
    else:
        normalized_score_0_100 = 100.0 * (
            (hits_at_k - BASELINE_SOLUTION_SCORE)
            / (AUTHOR_SOLUTION_SCORE - BASELINE_SOLUTION_SCORE)
        )
        normalized_score_0_100 = _clip_score_0_100(normalized_score_0_100)

    metrics: Dict[str, float | int | bool] = {
        "n_samples": int(len(df)),
        "k_true": k_true,
        "k_pred": k_pred,
        "k_match": bool(k_true == k_pred),
        "hits_at_k": hits_at_k,
    }
    if include_score_0_100:
        metrics["score_0_100"] = normalized_score_0_100
    return metrics


def validate_and_score_submission(ytest_path: Path, submission_path: Path) -> Dict[str, object]:
    try:
        ytest = pd.read_csv(ytest_path)
    except Exception as exc:  # pragma: no cover
        raise ValueError(f"Error reading ytest.csv: {exc}") from exc

    try:
        submission = pd.read_csv(submission_path)
    except Exception as exc:  # pragma: no cover
        raise ValueError(f"Error reading submission.csv: {exc}") from exc

    if ytest.empty:
        raise ValueError("ytest.csv is empty")
    if submission.empty:
        raise ValueError("submission.csv is empty")

    required_ytest_cols = {"id", "is_outlier_true", "Usage"}
    required_submission_cols = {"id", "is_outlier"}

    missing_ytest_cols = required_ytest_cols - set(ytest.columns)
    if missing_ytest_cols:
        raise ValueError(f"Columns missing from ytest.csv: {sorted(missing_ytest_cols)}")

    missing_submission_cols = required_submission_cols - set(submission.columns)
    if missing_submission_cols:
        raise ValueError(
            f"Columns missing from submission.csv: {sorted(missing_submission_cols)}"
        )

    ytest = ytest[["id", "is_outlier_true", "Usage"]].copy()
    submission = submission[["id", "is_outlier"]].copy()

    ytest["id"] = ytest["id"].astype(str).str.strip()
    ytest["Usage"] = ytest["Usage"].astype(str).str.strip().str.title()

    submission["id"] = submission["id"].astype(str).str.strip()

    if (ytest["id"] == "").any():
        raise ValueError("ytest.csv contains empty ids")
    if (submission["id"] == "").any():
        raise ValueError("submission.csv contains empty ids")

    invalid_usage = sorted(set(ytest["Usage"]) - {"Public", "Private"})
    if invalid_usage:
        raise ValueError(f"Invalid Usage values in ytest.csv: {invalid_usage}")

    if ytest["id"].duplicated().any():
        duplicates = (
            ytest.loc[ytest["id"].duplicated(keep=False), "id"]
            .unique()
            .tolist()
        )
        suffix = " ..." if len(duplicates) > 10 else ""
        raise ValueError(
            f"ytest.csv contains duplicate ids: {duplicates[:10]}{suffix} "
            f"({len(duplicates)} in total)"
        )

    if submission["id"].duplicated().any():
        duplicates = (
            submission.loc[submission["id"].duplicated(keep=False), "id"]
            .unique()
            .tolist()
        )
        suffix = " ..." if len(duplicates) > 10 else ""
        raise ValueError(
            f"submission.csv contains duplicate ids: {duplicates[:10]}{suffix} "
            f"({len(duplicates)} in total)"
        )

    y_ids = set(ytest["id"])
    s_ids = set(submission["id"])

    missing_ids = y_ids - s_ids
    if missing_ids:
        miss = sorted(list(missing_ids))
        suffix = " ..." if len(miss) > 10 else ""
        raise ValueError(
            f"ids missing from submission.csv: {miss[:10]}{suffix} "
            f"({len(miss)} in total)"
        )

    extra_ids = s_ids - y_ids
    if extra_ids:
        extra = sorted(list(extra_ids))
        suffix = " ..." if len(extra) > 10 else ""
        raise ValueError(
            f"submission.csv contains extra ids: {extra[:10]}{suffix} "
            f"({len(extra)} in total)"
        )

    ytest["is_outlier_true"] = _validate_binary_column(ytest["is_outlier_true"], "is_outlier_true")
    submission["is_outlier"] = _validate_binary_column(submission["is_outlier"], "is_outlier")

    merged = ytest.merge(submission, on="id", how="left", validate="one_to_one")

    if merged["is_outlier"].isna().any():
        raise ValueError("NaN in the predictions after the merge")

    k_true_total = int(merged["is_outlier_true"].sum())
    k_pred_total = int(merged["is_outlier"].sum())
    if k_pred_total != k_true_total:
        raise ValueError(
            "submission.csv must contain exactly K ones in the is_outlier column, "
            f"where K={k_true_total}. Currently: {k_pred_total}"
        )

    overall_metrics = _compute_hits_metrics(merged, include_score_0_100=True)
    public_metrics = _compute_hits_metrics(
        merged.loc[merged["Usage"] == "Public"], include_score_0_100=False
    )
    private_metrics = _compute_hits_metrics(
        merged.loc[merged["Usage"] == "Private"], include_score_0_100=False
    )
    overall_metrics["score_0_100"] = _clip_score_0_100(overall_metrics["score_0_100"])

    return {
        **overall_metrics,
        "n_total": int(len(merged)),
        "n_public": int((merged["Usage"] == "Public").sum()),
        "n_private": int((merged["Usage"] == "Private").sum()),
        "public": public_metrics,
        "private": private_metrics,
    }


def main(argv: List[str] | None = None) -> None:
    parser = argparse.ArgumentParser(description="Check the submission and compute hits@k")
    parser.add_argument("--ytest", type=str, default="y_test.csv", help="Path to y_test.csv")
    parser.add_argument(
        "--submission",
        type=str,
        default="submission.csv",
        help="Path to submission.csv",
    )
    parser.add_argument(
        "--save-json",
        type=str,
        default=None,
        help="Optional path for saving the metrics as JSON",
    )
    args, unknown = parser.parse_known_args(argv)
    unknown = _strip_jupyter_kernel_args(unknown)
    if unknown:
        parser.error(f"unrecognized arguments: {' '.join(unknown)}")

    ytest_path = Path(args.ytest)
    submission_path = Path(args.submission)

    if not ytest_path.exists():
        raise FileNotFoundError(f"ytest file not found: {ytest_path}")
    if not submission_path.exists():
        raise FileNotFoundError(f"submission file not found: {submission_path}")

    metrics = validate_and_score_submission(ytest_path, submission_path)
    print(json.dumps(metrics, indent=2))

    if args.save_json is not None:
        out_path = Path(args.save_json)
        with out_path.open("w", encoding="utf-8") as file:
            json.dump(metrics, file, indent=2)
        print(f"Metrics saved to JSON: {out_path}")


if __name__ == "__main__":
    main()
```

*(Translator's note: the original shows no output for this cell. As printed, `_compute_hits_metrics` refers to `AUTHOR_SOLUTION_SCORE` and `BASELINE_SOLUTION_SCORE`, which the cell does not define; only the `_PUBLIC` and `_PRIVATE` variants are defined.)*