{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "37bb766a",
      "metadata": {},
      "source": [
        "# Trabalho final: Modelos de recomendação e análise de espaços latentes\n",
        "\n",
        "**Disciplina:** *Aprendizado não-supervisionado* - *MC847*\\\n",
        "**Instituição:** *Universidade Estadual de Campinas*\\\n",
        "**Professor(a):** *Marcos Medeiros Raimundo*\\\n",
        "**Semestre:** *2s2026*\n",
        "\n",
        "| Nome | RA |\n",
        "|---|---|\n",
        "| *Integrante 1* | *RA 1* |\n",
        "| *Integrante 2* | *RA 2* |\n",
        "| *Integrante 3 (opcional)* | *RA 3* |\n",
        "\n",
        "Preencha a tabela acima com os nomes e RAs dos integrantes da equipe (2 ou 3 alunos, impreterivelmente).\n",
        "\n",
        "## Introdução\n",
        "\n",
        "Em sistemas de recomendação, a principal tarefa é aprender a preferência de um usuário a partir de um histórico parcial de interações.\n",
        "Isso pode ser visto como um problema de *completar informação*: dado o histórico parcial de um usuário, quais itens provavelmente faltam?\n",
        "Neste trabalho vocês vão prever notas não observadas com modelos que fazem suposições diferentes sobre essa estrutura e, depois, investigar o que o espaço aprendido contém e quais suposições o produziram.\n",
        "\n",
        "O trabalho é dividido em três etapas, sendo as duas primeiras comuns a todos os grupos e a terceira a ser definida por cada grupo.\n",
        "A primeira etapa consiste em definir uma baseline que será usada como referência para as demais etapas.\n",
        "A segunda consiste em realizar e analisar experimentos para entender como diferentes suposições sobre a estrutura latente dos dados afetam o desempenho do modelo.\n",
        "A última etapa é a contribuição de cada grupo para o problema, que pode ser uma proposta de melhoria do modelo ou uma análise mais profunda de algum aspecto do problema.\n",
        "\n",
        "As duas primeiras etapas têm códigos de referência, que podem ser usados como base para a implementação dos modelos.\n",
        "Vocês podem usar bibliotecas externas com implementações de modelos, mas devem deixar claro quais bibliotecas foram usadas e como elas foram usadas.\n",
        "Por exemplo, vocês podem reescrever o modelo de fatoração de matrizes em PyTorch ou Jax para melhorar a performance, desde que justifiquem a escolha e expliquem como a implementação foi mantida fiel ao modelo original.\n",
        "\n",
        "## Dataset\n",
        "\n",
        "Utilizaremos o conjunto de dados MovieLens 100K [(Harper e Konstan, 2015)](https://www.scirp.org/reference/referencespapers?referenceid=3586299): 100 mil notas inteiras de 1 a 5, dadas por 943 usuários a 1.682 filmes, com pelo menos 20 notas por usuário.\n",
        "A matriz tem densidade de cerca de 6%, cada usuário tem idade (em anos), gênero e ocupação (21 categorias) e cada filme tem título, ano e 19 indicadores de gênero, e um filme pode ter vários gêneros.\n",
        "\n",
        "Devido à natureza do problema, a ausência de notas não pode ser tratada como aleatória: Os motivos pelos quais um usuário não avaliou um item podem ser variados, incluindo falta de interesse, falta de exposição ao item ou simplesmente esquecimento.\n",
        "Isso significa que os dados de treinamento são enviesados e não representam a distribuição real de preferências dos usuários."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "93f516da",
      "metadata": {},
      "outputs": [],
      "source": [
        "import random\n",
        "from pathlib import Path\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import seaborn as sns\n",
        "import torch\n",
        "from numpy.typing import NDArray\n",
        "\n",
        "sns.set_theme(style=\"whitegrid\")\n",
        "\n",
        "SEED = 42\n",
        "DEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
        "MOVIELENS_URL = \"https://files.grouplens.org/datasets/movielens/ml-100k.zip\"\n",
        "DATA_DIR = \"./data/\"\n",
        "\n",
        "def set_seed(seed: int) -> None:\n",
        "    \"\"\"Set the random seed for reproducibility.\"\"\"\n",
        "    random.seed(seed)\n",
        "    np.random.seed(seed)\n",
        "    torch.manual_seed(seed)\n",
        "    if torch.cuda.is_available():\n",
        "        torch.cuda.manual_seed_all(seed)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "cc9f7984",
      "metadata": {},
      "outputs": [],
      "source": [
        "from pandas.core.frame import DataFrame\n",
        "\n",
        "GENRES = [\n",
        "    \"unknown\", \"Action\", \"Adventure\", \"Animation\", \"Children's\", \"Comedy\",\n",
        "    \"Crime\", \"Documentary\", \"Drama\", \"Fantasy\", \"Film-Noir\", \"Horror\",\n",
        "    \"Musical\", \"Mystery\", \"Romance\", \"Sci-Fi\", \"Thriller\", \"War\", \"Western\"\n",
        "]\n",
        "\n",
        "\n",
        "def load_movielens() -> tuple[DataFrame, DataFrame, DataFrame]:\n",
        "    \"\"\"Download and extract the MovieLens dataset.\"\"\"\n",
        "    import zipfile\n",
        "\n",
        "    import requests\n",
        "\n",
        "    data_dir = Path(DATA_DIR)\n",
        "    data_dir.mkdir(parents=True, exist_ok=True)\n",
        "\n",
        "    zip_path = data_dir / \"ml-100k.zip\"\n",
        "    target_dir = data_dir / \"ml-100k\"\n",
        "\n",
        "    if not target_dir.exists():\n",
        "        if not zip_path.exists():\n",
        "            print(\"Downloading MovieLens dataset...\")\n",
        "            response = requests.get(MOVIELENS_URL)\n",
        "\n",
        "            with open(zip_path, \"wb\") as f:\n",
        "                f.write(response.content)\n",
        "\n",
        "            print(\"Download complete.\")\n",
        "\n",
        "        print(\"Extracting MovieLens dataset...\")\n",
        "        with zipfile.ZipFile(zip_path) as zip_ref:\n",
        "            zip_ref.extractall(data_dir)\n",
        "\n",
        "        print(f\"Extracted to {data_dir}.\")\n",
        "\n",
        "    ratings = pd.read_csv(\n",
        "            target_dir / \"u.data\",\n",
        "            sep=\"\\t\",\n",
        "            names=[\"user_id\", \"item_id\", \"rating\", \"timestamp\"],\n",
        "            dtype={\"user_id\": np.int32, \"item_id\": np.int32, \"rating\": np.float32, \"timestamp\": np.int64},\n",
        "        )\n",
        "\n",
        "    users = pd.read_csv(\n",
        "            target_dir / \"u.user\",\n",
        "            sep=\"|\",\n",
        "            names=[\"user_id\", \"age\", \"gender\", \"occupation\", \"zip_code\"],\n",
        "            dtype={\"user_id\": np.int32, \"age\": np.int32, \"gender\": str, \"occupation\": str, \"zip_code\": str},\n",
        "        )\n",
        "\n",
        "    item_cols = [\"item_id\", \"title\", \"release_date\", \"video_release_date\", \"imdb_url\", *GENRES]\n",
        "\n",
        "    movies = pd.read_csv(\n",
        "        target_dir / \"u.item\",\n",
        "        sep=\"|\",\n",
        "        encoding=\"latin-1\",\n",
        "        names=item_cols,\n",
        "        usecols=[\"item_id\", \"title\", \"release_date\", *GENRES],\n",
        "        dtype={\"item_id\": np.int32, \"title\": str, \"release_date\": str},\n",
        "    )\n",
        "\n",
        "    movies[\"year\"] = pd.to_numeric(movies[\"title\"].str.extract(r\"\\((\\d{4})\\)\")[0], errors=\"coerce\")\n",
        "\n",
        "    movies[\"genres\"] = movies[GENRES].apply(\n",
        "        lambda row: \"|\".join([g for g in GENRES if row[g] == 1]), axis=1\n",
        "    )\n",
        "\n",
        "    # Normalização para que os IDs de usuário e item comecem em 0\n",
        "    ratings[\"user_id\"] -= 1\n",
        "    ratings[\"item_id\"] -= 1\n",
        "    users[\"user_id\"] -= 1\n",
        "    movies[\"item_id\"] -= 1\n",
        "\n",
        "    return ratings, users, movies\n",
        "\n",
        "\n",
        "def ratings_to_matrix(ratings: pd.DataFrame, fill_value: float = float(\"nan\")) -> NDArray[np.float64]:\n",
        "    \"\"\"Convert ratings DataFrame to a user-item matrix.\"\"\"\n",
        "    num_users = ratings[\"user_id\"].max() + 1\n",
        "    num_items = ratings[\"item_id\"].max() + 1\n",
        "\n",
        "    user_item_matrix = np.full((num_users, num_items), fill_value, dtype=np.float64)\n",
        "    user_item_matrix[ratings[\"user_id\"], ratings[\"item_id\"]] = ratings[\"rating\"].values.astype(np.float64)\n",
        "\n",
        "    return user_item_matrix\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "dff62670",
      "metadata": {},
      "outputs": [],
      "source": [
        "ratings, users, movies = load_movielens()\n",
        "density = len(ratings) / (len(users) * len(movies))\n",
        "print(f\"notas={len(ratings):,}; usuários={len(users):,}; filmes={len(movies):,}; densidade={density:.3f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "dcd6e449",
      "metadata": {},
      "source": [
        "## Formulação do problema\n",
        "\n",
        "Uma suposição comum em sistemas de recomendação é que a preferência de um usuário por um item pode ser representada como uma função sobre um espaço latente compartilhado entre usuários e itens.\n",
        "Ou seja, considerando que cada usuário e cada item podem ser representados como vetores em um espaço latente, a nota de um usuário para um item pode ser modelada como uma função que mede a compatibilidade entre os vetores do usuário e do item.\n",
        "\n",
        "Utilizaremos uma formulação geral para este problema, que pode acomodar diferentes modelos e suposições sobre como os vetores latentes são gerados e como a função de compatibilidade é definida.\n",
        "Considere um conjunto de usuários $U$ e um conjunto de itens $I$, e seja $R$ a matriz de notas, a nota do usuário $u$ para o item $i$, $R_{ui}$ é modelada como:\n",
        "\n",
        "$$\n",
        "R_{ui} = \\mu + b_u + b_i + f(\\mathbf{z}_u, \\mathbf{w}_i) + \\epsilon_{ui}\n",
        "$$\n",
        "onde $\\mu$ é a média global das notas, $b_u$ e $b_i$ são os *biases* do usuário e do item, respectivamente, $\\mathbf{z}_u$ e $\\mathbf{w}_i$ são os vetores latentes do usuário e do item, $f$ é a função de compatibilidade entre os vetores latentes e $\\epsilon_{ui}$ é o erro de modelagem para aquele par usuário-item.\n",
        "\n",
        "Um desafio importante é que a matriz de notas $R$ é muito esparsa, ou seja, grande parte dos pares usuário-item não possuem notas observadas. Isso significa que aprender a estrutura latente compartilhada pode ser difícil sem fazer suposições adicionais sobre como os vetores latentes são gerados e como a função de compatibilidade é definida.\n",
        "\n",
        "Vamos explorar como diferentes suposições afetam o desempenho do modelo e o que podemos aprender sobre a estrutura latente dos dados através da inspeção do espaço latente aprendido.\n",
        "\n",
        "> **Importante:** Os modelos têm uma abordagem transdutiva, ou seja, eles aprendem os vetores latentes apenas para os usuários e itens presentes no conjunto de treinamento.\n",
        "As features de um usuário não são as notas que ele deu, mas sim apenas um identificador único, e o mesmo vale para os itens.\n",
        "Isso permite que o modelo possua parâmetros específicos para cada usuário, mas também significa que ele não pode generalizar para novos usuários ou itens sem re-treinamento.\n",
        "\n",
        "## Diretrizes para o trabalho\n",
        "\n",
        "Para cada experimento, sigam as seguintes diretrizes:\n",
        "\n",
        "> 1. Escrevam a hipótese e a previsão **antes** de rodar cada experimento.\n",
        "> 2. Resultados negativos ou inconclusivos valem, se vierem com **análise crítica**.\n",
        "> 3. Todos os modelos, inclusive os mais simples, passam pelo mesmo protocolo (a seguir).\n",
        "> 4. Separem o que um resultado diz sobre a previsão do que ele diz sobre a representação latente.\n",
        "> 5. Citem fontes sempre que necessário, seja para motivar uma escolha, ou para justificar uma decisão.\n",
        "\n",
        "A nota não depende da métrica obtida pelo modelo, mas da qualidade das hipóteses, aplicação correta do protocolo e da análise das evidências a favor, contra ou insuficientes.\n",
        "\n",
        "### Protocolo experimental\n",
        "\n",
        "Os experimentos devem ser controlados, com uma variável de interesse por vez e com o mesmo protocolo:\n",
        "\n",
        "**Split dos dados.** O conjunto de dados deve ser dividido em treino e teste, com 80% das notas para treino e 20% para teste.\n",
        "A divisão deve ser feita de forma estratificada por usuário, garantindo que cada usuário tenha pelo menos uma nota em cada conjunto, e todo filme deve ter pelo menos uma nota no conjunto de treino.\n",
        "\n",
        "**Métricas.** RMSE (principal) e MAE, com previsões truncadas em $[1,5]$, e NLL quando o modelo define uma verossimilhança comparável.\n",
        "\n",
        "**Seleção de hiperparâmetros.** O conjunto de teste é usado uma única vez, para reportar o desempenho final.\n",
        "Os hiperparâmetros (dimensão latente, regularização, variância do ruído, número de épocas e assim por diante) são escolhidos por *cross-validation* com $K$ folds dentro do conjunto de treino, respeitando as mesmas restrições da divisão de treino e teste.\n",
        "Depois da escolha, o modelo é treinado com todo o treino e avaliado no teste.\n",
        "\n",
        "**Incerteza.** Compare resultados com **intervalos de confiança** sobre as métricas.\n",
        "Vocês podem usar *cross-validation*, reportando a média e o desvio padrão das métricas entre os folds, que refletem a variação causada pelos dados de treino e pela inicialização, ou *bootstrapping*, amostrando com reposição o conjunto de teste e reportando a média e o intervalo de percentis de 95%, que reflete a variação da amostra de teste."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "414fea8b",
      "metadata": {},
      "outputs": [],
      "source": [
        "def split_ratings(\n",
        "    ratings: pd.DataFrame,\n",
        "    test_frac: float = 0.2,\n",
        "    seed: int = 0,\n",
        ") -> tuple[pd.DataFrame, pd.DataFrame]:\n",
        "    df = (\n",
        "        ratings.sample(frac=1, random_state=seed)\n",
        "               .sort_values(\"user_id\", kind=\"stable\")\n",
        "               .reset_index(drop=True)\n",
        "    )\n",
        "\n",
        "    grouped = df.groupby(\"user_id\")\n",
        "    pos = (\n",
        "        grouped.cumcount() /\n",
        "        grouped[\"user_id\"].transform(\"count\")\n",
        "    ).to_numpy()\n",
        "\n",
        "    splits = np.zeros(len(df), dtype=np.int8)\n",
        "    splits[pos >= (1 - test_frac)] = 1\n",
        "\n",
        "    # Trocar a primeira aparição de cada filme para o conjunto de treino\n",
        "    first_appearance = df.groupby(\"item_id\").head(1).index\n",
        "    splits[first_appearance] = 0\n",
        "\n",
        "    train_df = df.iloc[splits == 0]\n",
        "    test_df = df.iloc[splits == 1]\n",
        "\n",
        "    return train_df, test_df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "06179a16",
      "metadata": {},
      "outputs": [],
      "source": [
        "from typing import Protocol\n",
        "\n",
        "from sklearn.decomposition import PCA\n",
        "from sklearn.metrics import (\n",
        "    mean_absolute_error,\n",
        "    root_mean_squared_error,\n",
        ")\n",
        "\n",
        "# Funções auxiliares\n",
        "\n",
        "def clip_ratings(prediction: NDArray[np.float64]) -> NDArray[np.float64]:\n",
        "    return np.clip(np.asarray(prediction, dtype=np.float64), 1.0, 5.0)\n",
        "\n",
        "\n",
        "class Model(Protocol):\n",
        "    def fit(self, train_df: DataFrame) -> None: ...\n",
        "\n",
        "    def predict(\n",
        "        self, user_id: NDArray[np.intp], item_id: NDArray[np.intp],\n",
        "    ) -> NDArray[np.float64]: ...\n",
        "\n",
        "\n",
        "class ProbabilisticModel(Model, Protocol):\n",
        "    def predict_log_proba(\n",
        "        self,\n",
        "        user_id: NDArray[np.intp],\n",
        "        item_id: NDArray[np.intp],\n",
        "        ratings: NDArray[np.float64],\n",
        "    ) -> NDArray[np.float64]: ...\n",
        "\n",
        "\n",
        "def evaluate_model(\n",
        "    model: Model | ProbabilisticModel, **dfs: DataFrame\n",
        ") -> dict[str, float]:\n",
        "    metrics: dict[str, float] = {}\n",
        "    \n",
        "    for name, df in dfs.items():\n",
        "        users = df[\"user_id\"].to_numpy(dtype=np.intp)\n",
        "        items = df[\"item_id\"].to_numpy(dtype=np.intp)\n",
        "        ratings = df[\"rating\"].to_numpy(dtype=np.float64)\n",
        "\n",
        "        preds = model.predict(users, items)\n",
        "        preds_clipped = clip_ratings(preds)\n",
        "\n",
        "        metrics |= {\n",
        "            f\"{name}_rmse\": root_mean_squared_error(ratings, preds_clipped),\n",
        "            f\"{name}_mae\": mean_absolute_error(ratings, preds_clipped),\n",
        "            # Vocês podem adicionar mais métricas aqui\n",
        "        }\n",
        "\n",
        "        if hasattr(model, \"predict_log_proba\"):\n",
        "            log_proba = model.predict_log_proba(users, items, ratings)\n",
        "            metrics[f\"{name}_nll\"] = float(-np.mean(log_proba))\n",
        "\n",
        "    return metrics\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "ca2fd483",
      "metadata": {},
      "outputs": [],
      "source": [
        "train_df, test_df = split_ratings(ratings, seed=SEED)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "ad8eb3f2",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Explorem os dados, plotem gráficos, calculem estatísticas, para entender melhor o problema,\n",
        "# os dados disponíveis, a distribuição das notas, a distribuição de usuários e itens, e assim por diante."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5d432a85",
      "metadata": {},
      "source": [
        "## Etapa 1: Baselines\n",
        "\n",
        "O primeiro passo para entender o problema é estabelecer uma baseline, um modelo simples que sirva como referência da performance que podemos esperar de um modelo mais sofisticado.\n",
        "Podemos simplificar a formulação geral do problema, considerando apenas a média global das notas e os *biases* do usuário e do item, sem considerar os vetores latentes:\n",
        "\n",
        "$$R_{ui} = \\mu + b_u + b_i + \\epsilon_{ui}$$\n",
        "\n",
        "Como o modelo remove a interação entre usuários e itens, ele não é capaz de capturar a estrutura latente compartilhada entre eles, mas ainda assim pode fornecer uma boa referência para avaliar o desempenho de modelos mais complexos.\n",
        "[Koren, Y. et. al.](https://ieeexplore.ieee.org/document/5197422) discute este exato modelo como uma aproximação de primeira ordem para o problema, que pode ser melhorada com a adição de termos de segunda ordem, como os vetores latentes introduzidos na formulação geral.\n",
        "Recomenda-se a leitura do artigo para entender melhor a motivação e interpretação do problema e as soluções propostas.\n",
        "\n",
        "Devido à linealidade do modelo,podemos aprender os *biases* do usuário e do item de forma eficiente, mesmo com a matriz de notas esparsa, através de regressão linear regularizada.\n",
        "\n",
        "> **Tarefa 1.1**\n",
        "> Treine o modelo de baseline e avalie seu desempenho, reportando RMSE e MAE no conjunto de teste.\n",
        "> Lembre-se de usar o mesmo protocolo experimental para todos os modelos para encontrar o melhor conjunto de hiperparâmetros.\n",
        "> Após treinar o modelo, analise os *biases* aprendidos e discuta o que eles podem nos dizer sobre a preferência dos usuários e a popularidade dos itens."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "19686fbf",
      "metadata": {},
      "outputs": [],
      "source": [
        "class BaselineModel:\n",
        "    def __init__(\n",
        "        self,\n",
        "        user_reg: float = 10.0,\n",
        "        item_reg: float = 10.0,\n",
        "        n_epochs: int = 20,\n",
        "    ) -> None:\n",
        "        self.user_reg = user_reg\n",
        "        self.item_reg = item_reg\n",
        "        self.n_epochs = n_epochs\n",
        "\n",
        "    def fit(self, train_df: pd.DataFrame) -> None:\n",
        "        ratings = train_df[\"rating\"].to_numpy(dtype=np.float64)\n",
        "        users = train_df[\"user_id\"].to_numpy(dtype=np.intp)\n",
        "        items = train_df[\"item_id\"].to_numpy(dtype=np.intp)\n",
        "\n",
        "        self.global_mean = float(ratings.mean())\n",
        "        self.user_bias = np.zeros(int(users.max()) + 1, dtype=np.float64)\n",
        "        self.item_bias = np.zeros(int(items.max()) + 1, dtype=np.float64)\n",
        "\n",
        "        n_users = len(self.user_bias)\n",
        "        n_items = len(self.item_bias)\n",
        "\n",
        "        user_counts = np.bincount(users, minlength=n_users)\n",
        "        item_counts = np.bincount(items, minlength=n_items)\n",
        "\n",
        "        for _ in range(self.n_epochs):\n",
        "            # Os vieses de usuário e item são atualizados alternadamente, mantendo o outro fixo.\n",
        "            # cada um resolve \\min_b (sum_{(u,i) in R} (r_{ui} - mu - b_u - b_i)^2 + lambda * b^2)\n",
        "            # Que é equivalente à agrupar por usuário e item, e então resolver a equação normal para cada um.\n",
        "\n",
        "            user_residuals = ratings - self.global_mean - self.item_bias[items]\n",
        "            self.user_bias = np.bincount(\n",
        "                users, weights=user_residuals, minlength=n_users\n",
        "            ) / (user_counts + self.user_reg)\n",
        "\n",
        "            item_residuals = ratings - self.global_mean - self.user_bias[users]\n",
        "            self.item_bias = np.bincount(\n",
        "                items, weights=item_residuals, minlength=n_items\n",
        "            ) / (item_counts + self.item_reg)\n",
        "\n",
        "    def predict(\n",
        "        self,\n",
        "        user_id: NDArray[np.intp],\n",
        "        item_id: NDArray[np.intp],\n",
        "    ) -> NDArray[np.float64]:\n",
        "        return (\n",
        "            self.global_mean\n",
        "            + self.user_bias[user_id]\n",
        "            + self.item_bias[item_id]\n",
        "        )"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "d002ed38",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Faça o treinamento e a avaliação do modelo baseline aqui"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "540c85bc",
      "metadata": {},
      "source": [
        "**Discussão:**\n",
        "*Escreva suas conclusões aqui*"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0619f34c",
      "metadata": {},
      "source": [
        "---\n",
        "### Introduzindo vetores latentes\n",
        "\n",
        "A baseline linear captura diferenças sistemáticas entre usuários e itens, mas não explica por que determinados grupos de usuários atribuem notas semelhantes a determinados grupos de itens.\n",
        "Para investigar essa estrutura, podemos olhar para a matriz de notas como um todo e perguntar se sua variabilidade pode ser explicada por um número reduzido de dimensões.\n",
        "Essa perspectiva permite explorar interações globais nos dados e construir uma representação mais compacta das avaliações, codificando usuários e itens em um espaço latente compartilhado.\n",
        "\n",
        "Uma maneira de investigar essa hipótese é completar a matriz de notas e aplicar Análise de Componentes Principais (PCA, *Principal Component Analysis*).\n",
        "Para isso, substituímos as notas ausentes pela média do respectivo item, calculada apenas com as avaliações observadas no conjunto de treinamento.\n",
        "Em seguida, o PCA identifica as principais direções de variação da matriz completada, permitindo representar cada usuário e cada item como vetores em um espaço de baixa dimensionalidade.\n",
        "\n",
        "Seja $\\tilde{R}$ a matriz de notas em que cada nota ausente foi substituída por um valor de imputação, por exemplo, a média do item.\n",
        "O PCA com $k$ componentes aproxima $\\tilde{R}$ por uma matriz de posto $k$:\n",
        "\n",
        "$$\\hat{R}_{ui} = \\bar{R}_i + \\mathbf{z}_u^\\top \\mathbf{w}_i$$\n",
        "\n",
        "onde $\\bar{R}_i$ é a média da coluna $i$ de $\\tilde{R}$, $\\mathbf{z}_u$ e $\\mathbf{w}_i$ são as coordenadas do usuário $u$ e do item $i$ nos $k$ primeiros componentes principais.\n",
        "Sob a perspectiva do PCA, um usuário é codificado como um vetor de pesos que indica a importância de cada componente principal para aquele usuário, que então é decodificado em uma nota para cada item, através da combinação linear com a transformação dos itens no mesmo espaço latente.\n",
        "\n",
        "> **Tarefa 1.2**\n",
        "Treine o modelo de PCA com imputação de notas ausentes, utilizando a média de cada item para preencher os valores faltantes. Utilize o protocolo experimental para selecionar o número de componentes principais $k$ e então reporte RMSE e MAE no conjunto de teste para o valor de $k$ selecionado.\n",
        ">\n",
        "> Analise e discuta também como o desempenho varia com o número de componentes, discutindo se o aumento da dimensionalidade melhora as previsões, produz ganhos decrescentes ou degrada o desempenho.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "c1a67624",
      "metadata": {},
      "outputs": [],
      "source": [
        "class ImputedPCA:\n",
        "    def __init__(self, n_components: int = 10) -> None:\n",
        "        self.n_components = n_components\n",
        "        self.user_embeddings: NDArray[np.float64]\n",
        "        self.item_embeddings: NDArray[np.float64]\n",
        "        self.item_means: NDArray[np.float64]\n",
        "        self.pca = PCA(n_components=n_components)\n",
        "\n",
        "    def fit(self, train_df: pd.DataFrame) -> None:\n",
        "        X = ratings_to_matrix(train_df)\n",
        "\n",
        "        # Imputação de valores ausentes com a média do item\n",
        "        observed = ~np.isnan(X)\n",
        "        global_mean = float(X[observed].mean())\n",
        "        counts = observed.sum(axis=0)\n",
        "        sums = np.where(observed, X, 0.0).sum(axis=0)\n",
        "\n",
        "        self.item_means = np.divide(\n",
        "            sums,\n",
        "            counts,\n",
        "            out=np.full(X.shape[1], global_mean),\n",
        "            where=counts > 0,\n",
        "        )\n",
        "\n",
        "        rows, cols = np.where(~observed)\n",
        "        X[rows, cols] = self.item_means[cols]\n",
        "\n",
        "        user_scores = self.pca.fit_transform(X)\n",
        "        self.user_embeddings = user_scores\n",
        "        self.item_embeddings = self.pca.components_.T\n",
        "\n",
        "    def predict(\n",
        "        self,\n",
        "        user_id: NDArray[np.intp],\n",
        "        item_id: NDArray[np.intp],\n",
        "    ) -> NDArray[np.float64]:\n",
        "        return (\n",
        "            self.pca.mean_[item_id]\n",
        "            + np.sum(\n",
        "                self.user_embeddings[user_id]\n",
        "                * self.item_embeddings[item_id],\n",
        "                axis=1,\n",
        "            )\n",
        "        )"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "9737954c",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Faça o treinamento e a avaliação do PCA imputado aqui"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5de7cdb0",
      "metadata": {},
      "source": [
        "**Discussão:**\n",
        "*Escreva suas conclusões aqui*"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "144ab559",
      "metadata": {},
      "source": [
        "---\n",
        "\n",
        "### Fatoração de matrizes (*Matrix Factorization*)\n",
        "\n",
        "A principal limitação do modelo anterior é que ele depende de notas imputadas, que não correspondem a avaliações reais e podem introduzir vieses no aprendizado dos vetores latentes.\n",
        "Para superar essa limitação, devemos aprender os vetores latentes diretamente a partir das notas observadas, sem precisar completar previamente a matriz.\n",
        "\n",
        "A fatoração de matrizes parte da hipótese de que a matriz de notas pode ser aproximada pelo produto de duas matrizes de menor dimensão, ou seja\n",
        "\n",
        "$$\\hat R = Z W^\\top,$$\n",
        "em que $Z\\in\\mathbb R^{|U|\\times k}$ e $W\\in\\mathbb R^{|I|\\times k}$ contêm os vetores latentes dos usuários e dos itens, respectivamente.\n",
        "Dessa forma, a nota prevista para um par usuário–item é dada pelo produto interno entre suas representações:\n",
        "\n",
        "$$\\hat R_{ui}=\\mathbf z_u^\\top\\mathbf w_i.$$\n",
        "\n",
        "Entretanto, como vimos na baseline, parte da variação das notas pode ser explicada por diferenças sistemáticas entre usuários e itens, independentemente de suas interações (por exemplo, alguns filmes são mais populares que outros).\n",
        "Para considerar esses efeitos, incorporamos a média global e os *biases* de usuário e item ao modelo, obtendo a formulação geral com $f(\\mathbf{z}_u, \\mathbf{w}_i) = \\mathbf z_u^\\top\\mathbf w_i$:\n",
        "\n",
        "$$\\hat R_{ui}=\\mu+b_u+b_i+\\mathbf z_u^\\top\\mathbf w_i.$$\n",
        "\n",
        "Para aprender os parâmetros, minimizamos o erro quadrático entre as notas observadas e as respectivas previsões, considerando apenas os pares usuário–item avaliados no conjunto de treinamento.\n",
        "Introduzimos também regularização $L_2$, que penaliza valores elevados dos *biases* e das normas dos vetores latentes:\n",
        "\n",
        "\\begin{aligned}\n",
        "\\min_{\\mu,b,Z_U,Z_I}\\quad\n",
        "&\\sum_{(u,i)\\in\\Omega}\n",
        "\\left(r_{ui}-\\mu-b_u-b_i-\\mathbf z_u^\\top\\mathbf w_i\\right)^2\\\\\n",
        "&+\\lambda_b\\left(\\sum_u b_u^2+\\sum_i b_i^2\\right)\\\\\n",
        "&+\\lambda_z\\left(\\sum_u\\|\\mathbf z_u\\|_2^2+\n",
        "\\sum_i\\|\\mathbf w_i\\|_2^2\\right),\n",
        "\\end{aligned}\n",
        "\n",
        "em que $\\Omega$ é o conjunto de avaliações observadas no treinamento, enquanto $\\lambda_b$ e $\\lambda_z$ controlam a intensidade da regularização dos *biases* e dos fatores latentes, respectivamente.\n",
        "Como o problema é não-convexo, sem solução fechada, usamos descida no gradiente para atualizar os parâmetros, iterativamente com base em pequenas amostras do conjunto de treinamento.\n",
        "\n",
        "> **Tarefa 1.3**\n",
        "Treine um modelo de fatoração de matrizes utilizando regularização $L_2$.\n",
        "Selecione o número de dimensões latentes $k$ e os coeficientes de regularização com base no RMSE.\n",
        "> \n",
        "> Reporte RMSE e MAE no conjunto de teste para a configuração selecionada.\n",
        "Analise como o desempenho varia com $k$ e com a regularização.\n",
        "Por fim, compare a fatoração com *biases* e sem *biases*, mantendo o mesmo protocolo experimental, para investigar a contribuição dos efeitos aditivos e das interações latentes para a qualidade das previsões.\n",
        ">\n",
        "> Siga o mesmo protocolo experimental utilizado para os modelos anteriores.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b400448a",
      "metadata": {},
      "outputs": [],
      "source": [
        "class MatrixFactorization:\n",
        "    user_embeddings: NDArray[np.float64]\n",
        "    item_embeddings: NDArray[np.float64]\n",
        "    global_mean: float\n",
        "    user_bias: NDArray[np.float64]\n",
        "    item_bias: NDArray[np.float64]\n",
        "\n",
        "    def __init__(\n",
        "        self,\n",
        "        n_factors: int = 10,\n",
        "        n_epochs: int = 20,\n",
        "        learning_rate: float = 0.01,\n",
        "        reg: float = 0.1,\n",
        "        bias_reg: float | None = None,\n",
        "        use_biases: bool = True,\n",
        "    ) -> None:\n",
        "        self.n_factors = n_factors\n",
        "        self.n_epochs = n_epochs\n",
        "        self.learning_rate = learning_rate\n",
        "        self.reg = reg\n",
        "        self.bias_reg = reg if bias_reg is None else bias_reg\n",
        "        self.use_biases = use_biases\n",
        "\n",
        "    def fit(self, train_df: pd.DataFrame) -> None:\n",
        "        ratings = train_df[\"rating\"].to_numpy(dtype=np.float64)\n",
        "        users = train_df[\"user_id\"].to_numpy(dtype=np.intp)\n",
        "        items = train_df[\"item_id\"].to_numpy(dtype=np.intp)\n",
        "\n",
        "        self.global_mean = float(ratings.mean())\n",
        "        n_users = int(users.max()) + 1\n",
        "        n_items = int(items.max()) + 1\n",
        "\n",
        "        self.user_embeddings = np.random.normal(scale=0.1, size=(n_users, self.n_factors))\n",
        "        self.item_embeddings = np.random.normal(scale=0.1, size=(n_items, self.n_factors))\n",
        "        self.user_bias = np.zeros(n_users, dtype=np.float64)\n",
        "        self.item_bias = np.zeros(n_items, dtype=np.float64)\n",
        "\n",
        "        for _ in range(self.n_epochs):\n",
        "            for n in np.random.permutation(len(ratings)):\n",
        "                u, i, r = users[n], items[n], ratings[n]\n",
        "                pred = (\n",
        "                    self.global_mean\n",
        "                    + self.user_bias[u]\n",
        "                    + self.item_bias[i]\n",
        "                    + np.dot(self.user_embeddings[u], self.item_embeddings[i])\n",
        "                )\n",
        "                err = r - pred\n",
        "\n",
        "                # Atualização dos vieses\n",
        "                if self.use_biases:\n",
        "                    self.user_bias[u] += self.learning_rate * (err - self.bias_reg * self.user_bias[u])\n",
        "                    self.item_bias[i] += self.learning_rate * (err - self.bias_reg * self.item_bias[i])\n",
        "\n",
        "                user_emb_old = self.user_embeddings[u].copy()\n",
        "                item_emb_old = self.item_embeddings[i].copy()\n",
        "\n",
        "                self.user_embeddings[u] += self.learning_rate * (err * item_emb_old - self.reg * user_emb_old)\n",
        "                self.item_embeddings[i] += self.learning_rate * (err * user_emb_old - self.reg * item_emb_old)\n",
        "\n",
        "    def predict(\n",
        "        self,\n",
        "        user_id: NDArray[np.intp],\n",
        "        item_id: NDArray[np.intp],\n",
        "    ) -> NDArray[np.float64]:\n",
        "        return (\n",
        "            self.global_mean\n",
        "            + self.user_bias[user_id]\n",
        "            + self.item_bias[item_id]\n",
        "            + np.sum(\n",
        "                self.user_embeddings[user_id] * self.item_embeddings[item_id],\n",
        "                axis=1,\n",
        "            )\n",
        "        )\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "10026bbc",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Faça o treinamento e a avaliação do modelo de fatoração de matrizes aqui"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6b780ee8",
      "metadata": {},
      "source": [
        "**Discussão:**\n",
        "*Escreva suas conclusões aqui*"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b3e62cee",
      "metadata": {},
      "source": [
        "---\n",
        "## Etapa 2: Suposições sobre a estrutura latente\n",
        "\n",
        "Na etapa anterior, investigamos modelos que estimam notas a partir de vieses aditivos e representações latentes de usuários e itens.\n",
        "A fatoração de matrizes procura representações que aproximem as avaliações observadas, mas não define, por si só, uma distribuição probabilística completa para os dados.\n",
        "\n",
        "Nesta etapa, exploraremos uma abordagem probabilística baseada no PCA Probabilístico (*Probabilistic Principal Component Analysis*, PPCA), adaptada ao problema de recomendação.\n",
        "Essa formulação descreve as avaliações por meio de variáveis latentes e ruído gaussiano, permitindo explicitar a variância do ruído, estimar a verossimilhança dos dados e investigar como as suposições probabilísticas afetam as previsões.\n",
        "\n",
        "Ao contrário das implementações clássicas de PPCA e VAE, em que o encoder é treinado para mapear as observações para a distribuição posterior, aqui, como os usuários do conjunto de dados são conhecidos durante o treinamento, podemos estimar um embedding para cada um deles e armazenar essas representações após o ajuste do modelo.\n",
        "Isso é equivalente a um *encoder* que mapeia cada usuário para uma distribuição posterior com o máximo de liberdade, sem lidar com features de usuário, vieses indutivos ou generalização para novos usuários.\n",
        "\n",
        "\n",
        "### PCA Probabilístico\n",
        "\n",
        "O PCA clássico identifica um subespaço linear de dimensão reduzida que representa a estrutura de variância dos dados. O PPCA estende essa ideia ao assumir que as observações são geradas por transformações lineares de variáveis latentes gaussianas, acrescidas de ruído gaussiano.\n",
        "\n",
        "No problema de recomendação, podemos representar a matriz de avaliações por $R\\in\\mathbb R^{n\\times m}$, em que $n$ é o número de usuários e $m$ é o número de filmes. Cada usuário corresponde a uma observação, e cada filme corresponde a uma dimensão observada. Como os usuários do conjunto de dados são conhecidos durante o treinamento, podemos estimar um embedding para cada um deles e armazenar essas representações após o ajuste do modelo.\n",
        "\n",
        "Seja $\\mathbf z_u\\in\\mathbb R^k$ a variável latente do usuário $u$, e seja $\\mathbf w_i\\in\\mathbb R^k$ o vetor de cargas associado ao filme $i$. O modelo assume que\n",
        "\n",
        "$$\n",
        "\\mathbf z_u\\sim\\mathcal N(\\mathbf 0,I_k)\n",
        "$$\n",
        "\n",
        "e que a avaliação do usuário $u$ para o filme $i$ segue\n",
        "\n",
        "$$\n",
        "R_{ui}\\mid\\mathbf z_u\n",
        "\\sim\n",
        "\\mathcal N\\left(\n",
        "\\mu_i+\\mathbf z_u^\\top\\mathbf w_i,\n",
        "\\sigma^2\n",
        "\\right),\n",
        "$$\n",
        "\n",
        "em que $\\mu_i$ é a média associada ao filme $i$, $k$ é a dimensão latente e $\\sigma^2$ é a variância do ruído, considerada fixa durante o treinamento.\n",
        "A média condicional da avaliação é, portanto,\n",
        "\n",
        "$$\n",
        "\\mathbb E[R_{ui}\\mid\\mathbf z_u]\n",
        "=\n",
        "\\mu_i+\\mathbf z_u^\\top\\mathbf w_i.\n",
        "$$\n",
        "\n",
        "Essa formulação representa as avaliações por meio de uma estrutura linear de baixa dimensão, sem introduzir vieses aditivos específicos de usuários.\n",
        "A diferença fundamental em relação à fatoração de matrizes determinística está na distribuição probabilística atribuída às avaliações e às variáveis latentes.\n",
        "\n",
        "O treinamento pode ser realizado pelo algoritmo EM (*Expectation-Maximization*).\n",
        "Seja $\\Omega_u$ o conjunto de filmes avaliados pelo usuário $u$ e $\\Omega_i$ o conjunto de usuários que avaliaram o filme $i$.\n",
        "Na etapa E, como o modelo é linear e gaussiano, a posterior de cada usuário, condicionada às avaliações observadas e aos parâmetros atuais, é gaussiana, $p(\\mathbf z_u\\mid R_u)=\\mathcal N(\\mathbf m_u, S_u)$, com\n",
        "\n",
        "$$\n",
        "S_u=\\left( I_k + \\frac{1}{\\sigma^2} \\sum_{i\\in\\Omega_u} \\mathbf w_i\\mathbf w_i^\\top\\right)^{-1},\n",
        "\\qquad\n",
        "\\mathbf m_u=\\frac{1}{\\sigma^2} \\, S_u\\sum_{i\\in\\Omega_u}\\mathbf w_i\\,(r_{ui}-\\mu_i).\n",
        "$$\n",
        "\n",
        "Na etapa M, atualizamos as médias dos filmes e os vetores de cargas para maximizar a verossimilhança esperada.\n",
        "Definindo $\\boldsymbol\\theta_i=(\\mu_i,\\mathbf w_i)$ e $\\boldsymbol\\phi_u=(1,\\mathbf z_u)$, cada filme resolve um sistema de mínimos quadrados em que as posteriores entram pelos momentos $\\mathbb E[\\boldsymbol\\phi_u]=(1,\\mathbf m_u)$ e $\\mathbb E[\\mathbf z_u\\mathbf z_u^\\top]=S_u+\\mathbf m_u\\mathbf m_u^\\top$:\n",
        "\n",
        "$$\n",
        "\\boldsymbol\\theta_i=\\left(\\sum_{u\\in\\Omega_i}\\mathbb E[\\boldsymbol\\phi_u\\boldsymbol\\phi_u^\\top]\\right)^{-1}\\sum_{u\\in\\Omega_i}r_{ui}\\,\\mathbb E[\\boldsymbol\\phi_u].\n",
        "$$\n",
        "\n",
        "Como a matriz de avaliações é esparsa, cada usuário contribui para a estimação apenas por meio dos filmes que avaliou. As avaliações ausentes são tratadas como dados não observados, e não como notas iguais a zero. Após o treinamento, as médias dos embeddings posteriores podem ser armazenadas como representações dos usuários, permitindo realizar previsões diretamente para os pares usuário–filme.\n",
        "\n",
        "Para uma avaliação observada $r_{ui}$, a densidade condicional é\n",
        "\n",
        "$$\n",
        "p(r_{ui}\\mid\\mathbf z_u)\n",
        "=\n",
        "\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\n",
        "\\exp\\left(\n",
        "-\\frac{\n",
        "(r_{ui}-\\mu_i-\\mathbf w_i^\\top\\mathbf z_u)^2\n",
        "}{2\\sigma^2}\n",
        "\\right).\n",
        "$$\n",
        "\n",
        "A log-verossimilhança conjunta das avaliações observadas, condicionada às variáveis latentes, é obtida somando as contribuições dos pares usuário–filme observados. O prior gaussiano das variáveis latentes também participa do modelo probabilístico, favorecendo representações compatíveis com a distribuição $\\mathcal N(\\mathbf 0,I_k)$.\n",
        "\n",
        "Para prever a avaliação de um par usuário–filme, utilizamos a média condicional calculada com o embedding posterior estimado:\n",
        "\n",
        "$$\n",
        "\\hat R_{ui}\n",
        "=\n",
        "\\mu_i+\\mathbf w_i^\\top\\mathbb E[\\mathbf z_u\\mid R_u],\n",
        "$$\n",
        "\n",
        "em que $R_u$ representa as avaliações observadas do usuário $u$. Essa previsão corresponde à média da distribuição preditiva quando os parâmetros do modelo são considerados fixos e a incerteza sobre a variável latente é integrada.\n",
        "\n",
        "A variância $\\sigma^2$ controla o ruído residual assumido pelo modelo. Valores menores atribuem maior concentração às avaliações em torno da média condicional; valores maiores representam maior dispersão. Como a variância é fixada durante o treinamento nesta implementação, ela deve ser selecionada como hiperparâmetro.\n",
        "\n",
        "> **Tarefa 2.1**\n",
        "Treine o modelo PPCA para a matriz esparsa de avaliações e selecione o número de dimensões latentes $k$ e a variância do ruído $\\sigma^2$ usando o protocolo experimental.\n",
        "Compare o RMSE com os modelos da etapa anterior e investigue se a configuração que minimiza o RMSE também produz a melhor log-verossimilhança negativa média (*Negative Log-Likelihood*, NLL).\n",
        ">\n",
        "> Discuta como a variância do ruído afeta a distribuição das avaliações, como o prior gaussiano influencia as representações latentes e quais são as implicações de estimar embeddings específicos para os usuários conhecidos durante o treinamento."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "188f15e4",
      "metadata": {},
      "outputs": [],
      "source": [
        "class PPCA:\n",
        "    user_embeddings: NDArray[np.float64]\n",
        "    item_embeddings: NDArray[np.float64]\n",
        "    item_means: NDArray[np.float64]\n",
        "    user_covariances: NDArray[np.float64]\n",
        "\n",
        "    def __init__(\n",
        "        self, n_factors: int = 10, n_epochs: int = 100, sigma2: float = 1.0\n",
        "    ) -> None:\n",
        "        self.n_factors = n_factors\n",
        "        self.n_epochs = n_epochs\n",
        "        self.sigma2 = sigma2\n",
        "\n",
        "    def fit(self, train_df: pd.DataFrame) -> None:\n",
        "        X = ratings_to_matrix(train_df)\n",
        "        observed = ~np.isnan(X)\n",
        "        values = np.where(observed, X, 0.0)\n",
        "        weights = observed.astype(np.float64)\n",
        "\n",
        "        counts = observed.sum(axis=0)\n",
        "        global_mean = values.sum() / observed.sum()\n",
        "        self.item_means = np.divide(\n",
        "            values.sum(axis=0),\n",
        "            counts,\n",
        "            out=np.full(X.shape[1], global_mean),\n",
        "            where=counts > 0,\n",
        "        )\n",
        "\n",
        "        n_users, n_items = X.shape\n",
        "        k = self.n_factors\n",
        "        sigma2 = self.sigma2\n",
        "        eye = np.eye(k)\n",
        "\n",
        "        self.item_embeddings = np.random.normal(\n",
        "            scale=0.1, size=(n_items, k)\n",
        "        )\n",
        "\n",
        "        item_outer = (\n",
        "            self.item_embeddings[:, :, None]\n",
        "            * self.item_embeddings[:, None, :]\n",
        "        )\n",
        "\n",
        "        for _ in range(self.n_epochs):\n",
        "            centered = np.where(observed, X - self.item_means, 0.0)\n",
        "\n",
        "            precision = eye + np.einsum(\"ui,ijk->ujk\", weights, item_outer) / sigma2\n",
        "\n",
        "            self.user_covariances = np.linalg.inv(precision)\n",
        "\n",
        "            rhs = centered @ self.item_embeddings / sigma2\n",
        "            self.user_embeddings = (\n",
        "                self.user_covariances @ rhs[:, :, None]\n",
        "            )[:, :, 0]\n",
        "\n",
        "            user_means = np.column_stack(\n",
        "                (np.ones(n_users), self.user_embeddings)\n",
        "            )\n",
        "\n",
        "            user_second = np.zeros((n_users, k + 1, k + 1))\n",
        "            user_second[:, 0, 0] = 1.0\n",
        "            user_second[:, 0, 1:] = self.user_embeddings\n",
        "            user_second[:, 1:, 0] = self.user_embeddings\n",
        "            user_second[:, 1:, 1:] = (\n",
        "                self.user_covariances\n",
        "                + self.user_embeddings[:, :, None]\n",
        "                * self.user_embeddings[:, None, :]\n",
        "            )\n",
        "\n",
        "            A = np.einsum(\"ui,uab->iab\", weights, user_second)\n",
        "\n",
        "            b = values.T @ user_means\n",
        "            theta = np.linalg.solve(A, b[:, :, None])[:, :, 0]\n",
        "\n",
        "            self.item_means = theta[:, 0]\n",
        "            self.item_embeddings = theta[:, 1:]\n",
        "\n",
        "            item_outer = (\n",
        "                self.item_embeddings[:, :, None]\n",
        "                * self.item_embeddings[:, None, :]\n",
        "            )\n",
        "\n",
        "        centered = np.where(observed, X - self.item_means, 0.0)\n",
        "\n",
        "        precision = eye + np.einsum(\"ui,ijk->ujk\", weights, item_outer) / sigma2\n",
        "\n",
        "        self.user_covariances = np.linalg.inv(precision)\n",
        "\n",
        "        rhs = centered @ self.item_embeddings / sigma2\n",
        "        self.user_embeddings = (\n",
        "            self.user_covariances @ rhs[:, :, None]\n",
        "        )[:, :, 0]\n",
        "\n",
        "    def predict(\n",
        "        self,\n",
        "        user_id: NDArray[np.intp],\n",
        "        item_id: NDArray[np.intp],\n",
        "    ) -> NDArray[np.float64]:\n",
        "        return self.item_means[item_id] + np.sum(\n",
        "            self.user_embeddings[user_id] * self.item_embeddings[item_id],\n",
        "            axis=1,\n",
        "        )\n",
        "\n",
        "    def predict_log_proba(\n",
        "        self,\n",
        "        user_id: NDArray[np.intp],\n",
        "        item_id: NDArray[np.intp],\n",
        "        ratings: NDArray[np.float64],\n",
        "    ) -> NDArray[np.float64]:\n",
        "        user_embeddings = self.user_embeddings[user_id]\n",
        "        item_embeddings = self.item_embeddings[item_id]\n",
        "        covariances = self.user_covariances[user_id]\n",
        "\n",
        "        means = self.item_means[item_id] + np.sum(\n",
        "            user_embeddings * item_embeddings, axis=1\n",
        "        )\n",
        "\n",
        "        variances = self.sigma2 + np.einsum(\n",
        "            \"ni,nij,nj->n\",\n",
        "            item_embeddings,\n",
        "            covariances,\n",
        "            item_embeddings,\n",
        "        )\n",
        "\n",
        "        return -0.5 * (\n",
        "            np.log(2 * np.pi * variances)\n",
        "            + (ratings - means) ** 2 / variances\n",
        "        )"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e0a454f5",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Faça o treinamento e a avaliação do PPCA aqui"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d6ff52d3",
      "metadata": {},
      "source": [
        "**Discussão:**\n",
        "*Escreva suas conclusões aqui*"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2b85525f",
      "metadata": {},
      "source": [
        "### Autoencoder variacional (VAE)\n",
        "\n",
        "O PPCA modela as avaliações por meio de uma transformação linear de variáveis latentes gaussianas. O VAE estende essa abordagem ao permitir um decodificador não linear e aprender uma distribuição latente específica para cada usuário.\n",
        "\n",
        "Seja $k$ a dimensão latente.\n",
        "O codificador é representado por uma matriz de embeddings de usuários, que fornece a média latente $\\mathbf m_u$ e uma matriz de parâmetros que fornece a variância diagonal $\\mathbf s_u^2$ de cada usuário.\n",
        "Assim, a distribuição latente aproximada é\n",
        "\n",
        "$$\n",
        "\\mathbf z_u\n",
        "\\sim\n",
        "\\mathcal N\\left(\n",
        "\\mathbf m_u,\n",
        "\\operatorname{diag}(\\mathbf s_u^2)\n",
        "\\right),\n",
        "$$\n",
        "que pode ser amostrada por meio do truque de reparametrização:\n",
        "\n",
        "$$\n",
        "\\mathbf z_u\n",
        "=\n",
        "\\mathbf m_u+\\mathbf s_u\\odot\\boldsymbol\\epsilon,\n",
        "\\qquad\n",
        "\\boldsymbol\\epsilon\\sim\\mathcal N(\\mathbf 0,I_k).\n",
        "$$\n",
        "\n",
        "O decodificador transforma essa amostra em previsões para todos os filmes através de uma função $g_\\theta:\\mathbb R^k\\to\\mathbb R^m$ parametrizada por $\\theta$:\n",
        "\n",
        "$$\n",
        "\\hat{\\mathbf R}_u = g_\\theta(\\mathbf z_u).\n",
        "$$\n",
        "\n",
        "Assumimos o seguinte modelo generativo:\n",
        "\n",
        "$$\n",
        "R_{ui}\\mid\\mathbf z_u\n",
        "\\sim\n",
        "\\mathcal N\\left(g_\\theta(\\mathbf z_u)_i,\\sigma^2\\right),\n",
        "$$\n",
        "\n",
        "em que $\\sigma^2$ é a variância do ruído das avaliações.\n",
        "Com um decodificador linear, $g_\\theta(\\mathbf z)=\\boldsymbol\\mu+W\\mathbf z$, a média condicional das avaliações é uma função linear da variável latente, como no PPCA.\n",
        "Um decodificador não linear permite representar relações mais complexas entre as dimensões latentes e as avaliações.\n",
        "\n",
        "O treinamento maximiza o limite inferior da log-verossimilhança (*Evidence Lower Bound*, ELBO), somado sobre os usuários:\n",
        "\n",
        "$$\n",
        "\\mathcal L_u\n",
        "=\n",
        "\\mathbb E_{q_\\phi(\\mathbf z_u)}\n",
        "\\left[\n",
        "-\\frac{1}{2\\sigma^2}\n",
        "\\sum_{i\\in\\Omega_u}\n",
        "(r_{ui}-g_\\theta(\\mathbf z_u)_i)^2\n",
        "\\right] -\\beta\\,\\mathrm{KL}\\left(\n",
        "q_\\phi(\\mathbf z_u)\n",
        "\\,\\|\\,\\mathcal N(\\mathbf 0,I_k)\n",
        "\\right),\n",
        "$$\n",
        "em que $\\Omega_u$ é o conjunto de filmes avaliados pelo usuário.\n",
        "Foram omitidas constantes que não dependem dos parâmetros do modelo.\n",
        "Com $\\beta=1$, obtemos o ELBO padrão, a menos dessas constantes, variar $\\beta$ altera o peso relativo da regularização latente.\n",
        "\n",
        "O primeiro termo incentiva o modelo a atribuir alta probabilidade às notas observadas, o segundo aproxima a distribuição latente de cada usuário do prior gaussiano padrão, regularizando as representações aprendidas.\n",
        "A esperança é estimada por amostragem da distribuição latente, permitindo propagar gradientes pelos parâmetros do codificador e do decodificador.\n",
        "\n",
        "Na predição, podemos usar a média latente $\\mathbf m_u$ como embedding do usuário e calcular\n",
        "\n",
        "$$\n",
        "\\hat R_{ui} = g_\\theta(\\mathbf m_u)_i.\n",
        "$$\n",
        "\n",
        "Essa escolha fornece uma previsão pontual determinística. A distribuição preditiva completa, por sua vez, também considera a variabilidade latente e o ruído das avaliações.\n",
        "\n",
        "Diferentemente do PPCA, que calcula a posterior de cada usuário analiticamente a partir das avaliações observadas, o VAE aprende por otimização uma distribuição latente específica para cada usuário, representada por uma média e uma variância. Essa formulação é transdutiva: as representações são aprendidas para os usuários presentes no treinamento e não se estendem automaticamente a usuários novos.\n",
        "\n",
        "> **Tarefa 2.2**\n",
        "Treine um VAE com decodificador linear e compare seu desempenho com o PPCA usando o mesmo número de dimensões latentes $k$ e a mesma variância do ruído $\\sigma^2$. Compare RMSE, MAE e NLL no conjunto de teste.\n",
        ">\n",
        "> Em seguida, substitua o decodificador linear por um não linear e investigue se a flexibilidade adicional melhora as previsões. Escolha um hiperparâmetros para variar isoladamente, $\\beta$, $\\sigma^2$ ou $k$, e aplique o protocolo experimental.\n",
        "Verifique a convergência pela curva de treinamento armazenada em `history` antes de comparar os modelos.\n",
        "\n",
        "Atenção, com um decoder não-linear, os itens "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "185c40d2",
      "metadata": {},
      "outputs": [],
      "source": [
        "import torch\n",
        "from torch import nn\n",
        "\n",
        "\n",
        "class VAE:\n",
        "    user_embeddings: NDArray[np.float64]\n",
        "    user_variances: NDArray[np.float64]\n",
        "    item_embeddings: NDArray[np.float64]\n",
        "    item_means: NDArray[np.float64]\n",
        "\n",
        "    def __init__(\n",
        "        self,\n",
        "        n_factors: int = 10,\n",
        "        n_epochs: int = 300,\n",
        "        sigma2: float = 1.0,\n",
        "        beta: float = 1.0,\n",
        "        hidden_dim: int | None = None,\n",
        "        learning_rate: float = 1e-2,\n",
        "    ) -> None:\n",
        "        self.n_factors = n_factors\n",
        "        self.n_epochs = n_epochs\n",
        "        self.sigma2 = sigma2\n",
        "        self.beta = beta\n",
        "        self.hidden_dim = hidden_dim\n",
        "        self.learning_rate = learning_rate\n",
        "        self.history: list[float] = []\n",
        "\n",
        "    def fit(self, train_df: pd.DataFrame) -> None:\n",
        "        X = ratings_to_matrix(train_df)\n",
        "        observed = ~np.isnan(X)\n",
        "        values = np.where(observed, X, 0.0)\n",
        "\n",
        "        X = torch.tensor(values, dtype=torch.float32)\n",
        "        mask = torch.tensor(observed, dtype=torch.float32)\n",
        "\n",
        "        n_users, n_items = X.shape\n",
        "        k = self.n_factors\n",
        "\n",
        "        self.user_means = nn.Parameter(torch.randn(n_users, k) * 0.01)\n",
        "        self.user_logvars = nn.Parameter(torch.zeros(n_users, k))\n",
        "\n",
        "        if self.hidden_dim is None:\n",
        "            self.decoder = nn.Linear(k, n_items)\n",
        "\n",
        "        else:\n",
        "            self.decoder = nn.Sequential(\n",
        "                nn.Linear(k, self.hidden_dim),\n",
        "                nn.Tanh(),\n",
        "                nn.Linear(self.hidden_dim, n_items),\n",
        "            )\n",
        "\n",
        "        # O viés do decodificador começa nas médias dos filmes, e o treino modela o desvio em relação a elas\n",
        "        counts = np.maximum(observed.sum(axis=0), 1)\n",
        "        item_means = np.where(observed.any(axis=0), values.sum(axis=0) / counts, values.sum() / observed.sum())\n",
        "        output_layer = self.decoder if self.hidden_dim is None else self.decoder[-1]\n",
        "        with torch.no_grad():\n",
        "            output_layer.bias.copy_(torch.tensor(item_means, dtype=torch.float32))\n",
        "\n",
        "        parameters = [\n",
        "            self.user_means,\n",
        "            self.user_logvars,\n",
        "            *self.decoder.parameters(),\n",
        "        ]\n",
        "        optimizer = torch.optim.Adam(parameters, lr=self.learning_rate)\n",
        "\n",
        "        self.history = []\n",
        "\n",
        "        for _ in range(self.n_epochs):\n",
        "            optimizer.zero_grad()\n",
        "\n",
        "            std = torch.exp(0.5 * self.user_logvars)\n",
        "            epsilon = torch.randn_like(std)\n",
        "            z = self.user_means + std * epsilon\n",
        "\n",
        "            predictions = self.decoder(z)\n",
        "\n",
        "            log_likelihood = -0.5 * (\n",
        "                mask\n",
        "                * (\n",
        "                    np.log(2 * np.pi * self.sigma2)\n",
        "                    + (X - predictions) ** 2 / self.sigma2\n",
        "                )\n",
        "            ).sum()\n",
        "\n",
        "            kl = -0.5 * (\n",
        "                1\n",
        "                + self.user_logvars\n",
        "                - self.user_means**2\n",
        "                - self.user_logvars.exp()\n",
        "            ).sum()\n",
        "\n",
        "            elbo = log_likelihood - self.beta * kl\n",
        "            loss = -elbo / n_users\n",
        "\n",
        "            loss.backward()\n",
        "            optimizer.step()\n",
        "\n",
        "            self.history.append((-loss.detach()).item())\n",
        "\n",
        "        with torch.no_grad():\n",
        "            self.user_embeddings = (\n",
        "                self.user_means.detach().cpu().numpy().astype(np.float64)\n",
        "            )\n",
        "            self.user_variances = (\n",
        "                self.user_logvars.exp()\n",
        "                .detach()\n",
        "                .cpu()\n",
        "                .numpy()\n",
        "                .astype(np.float64)\n",
        "            )\n",
        "\n",
        "            if self.hidden_dim is None:\n",
        "                self.item_embeddings = (\n",
        "                    self.decoder.weight.detach()\n",
        "                    .cpu().numpy().astype(np.float64)\n",
        "                )\n",
        "                self.item_means = (\n",
        "                    self.decoder.bias.detach()\n",
        "                    .cpu().numpy().astype(np.float64)\n",
        "                )\n",
        "\n",
        "    def _decode_at(self, z: torch.Tensor, items: torch.Tensor) -> torch.Tensor:\n",
        "        \"\"\"Decode only the entries (z[..., n, :], items[n]), without computing every movie.\"\"\"\n",
        "        if self.hidden_dim is None:\n",
        "            return (z * self.decoder.weight[items]).sum(-1) + self.decoder.bias[items]\n",
        "        first, _, last = self.decoder\n",
        "        return (torch.tanh(first(z)) * last.weight[items]).sum(-1) + last.bias[items]\n",
        "\n",
        "    def predict(\n",
        "        self,\n",
        "        user_id: NDArray[np.intp],\n",
        "        item_id: NDArray[np.intp],\n",
        "    ) -> NDArray[np.float64]:\n",
        "        users = torch.as_tensor(user_id, dtype=torch.long)\n",
        "        items = torch.as_tensor(item_id, dtype=torch.long)\n",
        "\n",
        "        with torch.no_grad():\n",
        "            z = torch.as_tensor(self.user_embeddings, dtype=torch.float32)[users]\n",
        "            result = self._decode_at(z, items)\n",
        "\n",
        "        return result.cpu().numpy().astype(np.float64)\n",
        "\n",
        "    def predict_log_proba(\n",
        "        self,\n",
        "        user_id: NDArray[np.intp],\n",
        "        item_id: NDArray[np.intp],\n",
        "        ratings: NDArray[np.float64],\n",
        "        n_samples: int = 100,\n",
        "        chunk: int = 2048,\n",
        "    ) -> NDArray[np.float64]:\n",
        "        \"\"\"Log of the Monte Carlo predictive density, computed in chunks of rating pairs.\"\"\"\n",
        "        users = torch.as_tensor(user_id, dtype=torch.long)\n",
        "        items = torch.as_tensor(item_id, dtype=torch.long)\n",
        "        ratings_t = torch.as_tensor(ratings, dtype=torch.float32)\n",
        "        user_means = torch.as_tensor(self.user_embeddings, dtype=torch.float32)\n",
        "        user_stds = torch.as_tensor(self.user_variances, dtype=torch.float32).sqrt()\n",
        "\n",
        "        result = []\n",
        "        with torch.no_grad():\n",
        "            for start in range(0, len(users), chunk):\n",
        "                u = users[start : start + chunk]\n",
        "                means, stds = user_means[u], user_stds[u]\n",
        "                z = means + stds * torch.randn(n_samples, *means.shape)\n",
        "                predictions = self._decode_at(z, items[start : start + chunk])\n",
        "\n",
        "                log_probs = -0.5 * (\n",
        "                    np.log(2 * np.pi * self.sigma2)\n",
        "                    + (ratings_t[start : start + chunk] - predictions) ** 2 / self.sigma2\n",
        "                )\n",
        "                result.append(torch.logsumexp(log_probs, dim=0) - np.log(n_samples))\n",
        "\n",
        "        return torch.cat(result).cpu().numpy().astype(np.float64)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "55e6f9e0",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Faça o treinamento e a avaliação do VAE aqui"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ffdd7e57",
      "metadata": {},
      "source": [
        "**Discussão:**\n",
        "*Escreva suas conclusões aqui*"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9bfc8f2b",
      "metadata": {},
      "source": [
        "---\n",
        "\n",
        "### Visualização dos espaços latentes\n",
        "\n",
        "Dos modelos treinados, quatro deles aprendem vetores latentes por usuário e por item: PCA com imputação, fatoração de matrizes, PPCA e VAE.\n",
        "Esses vetores podem ser analisados com as mesmas ferramentas, permitindo que possamos interpretar e comparar os espaços latentes aprendidos por cada modelo.\n",
        "Garante-se que pelo menos um modelo de cada tipo tenha sido treinado, e então podemos comparar os espaços latentes aprendidos.\n",
        "\n",
        "> **Tarefa 2.3**\n",
        "Visualize os espaços latentes de usuários e de itens dos quatro modelos, usando um $k$ de sua escolha, e relacione o que aparece nas figuras com informações que não foram usadas no treino: por exemplo, gênero, idade e ocupação dos usuários, e gêneros, ano e popularidade dos filmes.\n",
        ">\n",
        "> O espaço de cada modelo organiza essas informações, e isso pode ser medido além da inspeção visual?\n",
        "Os quatro espaços compartilham a mesma estrutura, ou cada suposição produz uma geometria diferente?\n",
        "\n",
        "\n",
        "**Ideias para começar.**\n",
        "\n",
        "Para visualizar, projete os vetores em duas dimensões com PCA, que é determinístico e preserva a geometria global, e compare com o t-SNE.\n",
        "Colorir os pontos por um metadado mostra se ele aparece, e um modelo supervisionado linear (`Ridge` para variáveis contínuas, `LogisticRegression` com `roc_auc_score` para categorias) pode medir quanto de informação sobre o metadado está presente no espaço latente original, sem precisar de uma projeção.\n",
        "\n",
        "Isso permite verificar se, a partir dos vetores latentes, é possível prever a idade do usuário, o gênero do filme ou a popularidade do filme, quantificando o quanto de informação sobre essas variáveis está presente no espaço latente.\n",
        "\n",
        "Ao comparar os espaços latentes, é importante lembrar que cada modelo aprende uma representação diferente, e que a correspondência entre os vetores latentes de um modelo e de outro não é trivial, pois os espaços latentes podem ser invariantes à rotação, reflexão e escala.\n",
        "Para comparar dois espaços é necessário alinhar os vetores latentes de um modelo com os do outro, ou através de análises de correspondência de vizinhos, como a correlação de distâncias entre pontos, ou através da sobreposição dos vizinhos mais próximos de cada filme em dois modelos."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "c0b9d7eb",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Faça as visualizações e análises dos espaços latentes aqui"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "88d287c4",
      "metadata": {},
      "source": [
        "**Discussão:**\n",
        "*Escreva suas conclusões aqui*"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c45ba2d5",
      "metadata": {},
      "source": [
        "---\n",
        "\n",
        "## Etapa 3: Projeto do grupo\n",
        "\n",
        "Nesta última etapa, cada um de vocês deve propor um projeto de investigação, que pode ser uma melhoria do modelo ou uma análise mais profunda de algum aspecto do problema.\n",
        "Após finalizar as duas primeiras etapas, verifiquem se há alguma pergunta que surgiu durante as análises e ainda não foi explorado, ou alguma hipótese que vocês gostariam de testar, por exemplo, sobre a geometria do espaço latente, uma formulação alternativa, ou um hiperparâmetro que não foi isolado.\n",
        "\n",
        "O importante é que o projeto esteja bem definido, com uma motivação clara, uma hipótese testável e um protocolo experimental que permita avaliar o impacto da mudança proposta.\n",
        "\n",
        "> **Tarefa 3.1**\n",
        "Escolha um projeto e escreva a pergunta, a hipótese, a previsão e o experimento antes de rodar.\n",
        "Defina o modelo referência, qual é a mudança proposta e qual resultado refutaria a hipótese de vocês.\n",
        "Ao final, compare o que foi previsto com o que foi observado e o que o resultado diz sobre a previsão sobre o espaço latente.\n",
        "\n",
        "**Projetos sugeridos.**\n",
        "Aqui está uma lista de ideias que podem ser exploradas, mas vocês podem propor outras, desde que estejam bem justificadas.\n",
        "\n",
        "*1. Ausência não aleatória.*\n",
        "Usuários ativos e filmes populares aparecem mais nos dados, o que cria um viés de observação.\n",
        "Pergunta: corrigir esse viés com pesos pelo inverso da probabilidade de observação [Schnabel et al.](https://arxiv.org/abs/1602.05352) melhora as previsões para filmes pouco avaliados, e o que acontece com os populares?\n",
        "\n",
        "*2. Metadados como features.*\n",
        "Ao invés de tratar cada usuário e item como um identificador único, podemos usar suas características como features para o \"codificador\".\n",
        "Pergunta: usar os metadados para definir a média do prior de $\\mathbf{z}$ melhora as previsões de usuários com pouco histórico?\n",
        "\n",
        "*3. Geometria da compatibilidade.*\n",
        "Todos os modelos utilizam o produto interno para medir a compatibilidade entre usuários e itens.\n",
        "Pergunta: trocar o produto interno por uma distância, ou por uma pequena rede que recebe $\\mathbf{z}_u$ e $\\mathbf{w}_i$, muda o que o espaço codifica e o erro de previsão? Ver [Rendle et al.](https://arxiv.org/abs/2005.09683) para um resultado que motiva a pergunta.\n",
        "\n",
        "*4. Estrutura do espaço latente.*\n",
        "A selecão de hiperparâmetros não observou a geometria do espaço latente, mas apenas o erro de previsão.\n",
        "Pergunta: o que muda no espaço latente aprendido do VAE ao variar $\\beta$, o tamanho do espaço ou a regularização, em termos de estabilidade de clusters, correlação com metadados e erro de previsão?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "05bf7992",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Desenvolva o experimento do seu grupo aqui"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f143a3fb",
      "metadata": {},
      "source": [
        "**Discussão:**\n",
        "*Escreva suas conclusões aqui*"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "base",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.13.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}