Laboratorium 12 + 13: grupowanie dokumentów tekstowych

Na dwóch kolejnych laboratoriach przeprowadzimy grupowanie dokumentów tekstowych z użyciem 3 metod:

• k-means
• EM
• DBSCAN (i jego wariantu HDBSCAN)

Narzędzia:

• Docker i docker-compose - do budowy środowiska z zainstalowanymi pakietami i udostępniania jupyter lab
• duckdb - do zapisu danych i wykonania kwerend SQL
• spaCy - do wektoryzacji i tokenizacji tekstów

Etapy:

1. Pobrane zostaną teksty 10000 artykułów z Wikipedii
2. Zbudowane zostaną ich reprezentacje wektorowe (wektory osadzeń i wektory TF-IDF)
3. Zastosujemy UMAP do redukcji wymiarów wektorów na potrzeby klasteryzacji i wizualizacji
4. Do dokumentów zostaną dodane etykiety klasteryzacji dla różnych algorytmów
5. Porównamy ich wyniki

Pobierz i rozpakuj archiwum clustering-jupyter.zip

Znajduje się tam plik pyproject.toml ze specyfikacją pakietów do zainstalowania

[project]
name = "clustering-nlp"
version = "0.1.0"
description = "Projekt na zajęcia: Eksploracja danych, 2026"
readme = "README.md"
requires-python = ">=3.12"

dependencies = [
    "pandas",
    "tabulate",
    "datasets",            # Pobieranie Wikipedii z Hugging Face
    "spacy",               # NLP i embeddingi statyczne
    "scikit-learn",        # K-means, GMM, DBSCAN
    "hdbscan",             # Warian DBSCAN
    "umap-learn",          # Redukcja wymiarów
    "plotly",              # Interaktywne wykresy
    "nbformat",            
    "tqdm",                # Paski postępu
    "matplotlib",          
    "seaborn",
    "wordcloud",             
    "jupyterlab",          
    "ipywidgets", 
    "jupysql",             # SQL w notatniku
    "duckdb",              # Baza OLAP 
    "duckdb-engine",
    "ollama",              # Opcjonalnie: biblioteka do komunikacji z lokalnym LLM
    "pyarrow",             
    "fastparquet",         
]

Plik Dockerfile

FROM python:3.12-slim

RUN apt-get update && apt-get install -y \
    curl \
    build-essential \
#    pandoc \
    && rm -rf /var/lib/apt/lists/*

# Instalacja uv
COPY --from=ghcr.io/astral-sh/uv:latest /uv /uvx /bin/

WORKDIR /app

COPY pyproject.toml .

# Instalacja zależności przez uv 
RUN uv sync

# Pobranie modelu spaCy
RUN uv run python -m spacy download pl_core_news_lg


WORKDIR /app/notebooks

oraz plik docker-compose.yaml

version: '3.8'

services:
  nlp-jupyter:
    build: .
    image: nlp-jupyter
    container_name: eksploracja_danych_nlp
    ports:
      - "9888:8888"
    volumes:
      - ./notebooks:/app/notebooks
      - ./hf_cache:/root/.cache/huggingface
    shm_size: '4gb'
    deploy:
      resources:
        limits:
          memory: 12gb
    command: >
      uv run jupyter lab 
      --ip=0.0.0.0 
      --allow-root 
      --no-browser 
      --IdentityProvider.token='nlp' 
      --notebook-dir=/app/notebooks

Uruchomienie

• aby zbudować obraz wydaj komendę docker compose build
• aby uruchomić kontener wydaj komendę docker compose up
• Po uruchomieniu otwórz link http://localhost:9888/lab?token=nlp
• Defaultowym katalogiem dla jupyetera lab jest notebooks. Znajduje się tam notatnik wikipedia.ipynb, w którym można kontynuować implementację kodu

Uwagi: Jako menadżer pakietów stosowany jest uv. Aby dodać pakiet należy w notatniku uruchomić komendę !uv add package

Kod, który pobiera 10 000 artykułów z Hugging Face znajduje się w notatniku wikipedia.ipynb. Kluczowym elementem jest losowanie wartości OFFSET - numeru dokumentu od którego zostanie ropoczęte pobieranie. Celem jest zróżnicowanie wyników i wizualizacji.

Zapisz korpus po załadowaniu dokumentów.

SpaCy Industrial-Strength Natural Language Processing to zaawansowana biblioteka, stworzona specjalnie z myślą o przemysłowych i produkcyjnych zastosowaniach NLP. W przeciwieństwie do akademickich narzędzi stawiających na mnogość algorytmów, spaCy skupia się na maksymalnej wydajności i dostarczaniu jednego, zoptymalizowanego rozwiązania dla każdego zadania. Za pomocą gotowych, wielojęzycznych modeli statystycznych i neuronowych biblioteka pozwala w ułamku sekundy przeprowadzić pełną analizę tekstu:

• podziału na tokeny (słowa) i zdania
• rozpoznawanie części mowy oraz form podstawowych (lematyzację)
• wykrywanie encji nazwanych (NER), takich jak imiona, daty czy lokalizacje
• związków składniowych

Ze względu na swoją szybkość, stabilność oraz łatwość integracji z frameworkami głębokiego uczenia, jest to obecnie jedna z najpopularniejszych bibliotek do budowania komercyjnych systemów wyszukiwania semantycznego, chatbotów oraz narzędzi do analizy dużych zbiorów tekstowych.

SpaCy oferuje ciągi przetwarzania dostosowane do konkretnego języka. Użyjemy modelu o dużej dokładności pl_core_news_lg

import spacy
 
nlp = spacy.load("pl_core_news_lg")
 
doc = nlp('Ala ma kota')
print(doc.vector.shape)
for t in doc:
    print(t,t.vector.shape)

W wyniku przetwarzania dokumentu dzielony jest on na tokeny, każdy token ma przypisane rózne atrybuty, w tym wektor osadzeń (o długości 300 dla wybranego wcześniej modelu). Wektory osadzeń tokenów po zsumowaniu dają wektor dokumentu. Wynik:

(300,)
Ala (300,)
ma (300,)
kota (300,)

Wyznaczanie wektorów osadzeń za pomocą spaCy może być czasochłonną operacją. SpaCy jest zainstalowane w obrazie dockera bez wsparcia dla GPU.

Zarejestrowane czasy to:

• 5 min z wydajnością 33 dokumentów/sekundę [wariant optymistyczny]
• 20 min z wydajnością 8.2 dokumentów/sekundę
• 47 minut z wydajnością 3.5 dokumentów/sekundę [wariant pesymistyczny, bez opcji n_process=-1]

Użyj poniższego kodu.

spacy_pipeline.py

import numpy as np
from tqdm import tqdm
 
processed_data = []
 
texts = df.text.to_list()
titles = df.title.to_list()
 
data_tuples = list(zip(texts, titles))
 
with nlp.select_pipes(enable=["tok2vec", "attribute_ruler", "lemmatizer"]): 
    for i, doc in enumerate(tqdm(nlp.pipe(texts, batch_size=50, n_process=-1), total=len(texts))):
 
        valid_tokens = [
            t for t in doc if not t.is_stop and t.is_alpha
        ]
 
        if valid_tokens:
            valid_vectors = [t.vector for t in valid_tokens]
            doc_vector = np.mean(valid_vectors, axis=0)
        else:
            doc_vector = np.zeros(nlp.vocab.vectors_length)
 
        processed_data.append(
            {
                "title": titles[i],
                "text": texts[i],
                "tokens": [t.lemma_.lower() for t in valid_tokens],  
                "vector": doc_vector,
            }
        )
 

• W celu przyspieszenia wykonywane są tylko części pipeline, np. nie jest budowany graf syntaktyczny zdania
• Obliczane są wektory osadzeń dla dokumentów z pominięciem stop words i tokenów zawierających liczby
• Ewentualnie można zmniejszyć batch_size, np. ustawiając na 20, 10, 5.
• Opcja n_process=-1 powoduje wykorzystanie wszystkich rdzeni procesora

Zapisz wyniki Po zakończeniu przetwarzania koniecznie zapisz wyniki (format dowolny). SpaCy zwraca wektory numpy. Należy je przekonwertować do postaci list, ponieważ ułatwi to późniejsze przetwarzanie.

df_vec = pd.DataFrame(processed_data)
df_vec['vector'] = df_vec['vector'].apply(lambda x: x.tolist() if isinstance(x, np.ndarray) else x)
 
# Na przykład zapis do Parquet
df_vec.to_parquet('wiki_vect.parquet', index=False)
print("Zapisano do Parquet (wektory jako listy).")

Wyznaczymy współczynniki TF-IDF dla słów w dokumentach, przyjmując wielkość słownika 300. W wyniku transformacji każdemu dokumentowi zostanie przypisany 300-elementowy rzadki wektor. Ta reprezentacja nazywana jest Bag of Words (BoW)

$$\text{tfidf}(t, d, D) = \text{tf}(t, d) \cdot \text{idf}(t, D)$$

• Term Frequency (TF) - to miara częstości termu $t$ w dokumencie $d$

$$\text{tf}(t, d) = \frac{n_{t,d}}{\sum_{k} n_{k,d}}$$

• Inverse Document Frequency (IDF) - to miara zawartości informacji termu na tle całego korpusu $D$

$$\text{idf}(t, D) = \log \frac{N}{1 + |\{d \in D : t \in d\}|}$$

Gdzie:

• $N$ – całkowita liczba dokumentów w korpusie,
• $|\{d \in D : t \in d\}|$ – liczba dokumentów zawierających słowo $t$.

• połącz tokeny spacjami
• zastosuj tfidf_vec = TfidfVectorizer(max_features=300) (z biblioteki scikit-learn)
• wynik transformacji (załóżmy, że jest on w zmiennej tfidf_sparse_matrix) zapisz w kolumnie “vector_sparse”. Przekonwertuj do postaci listy
• zapisz dane

Słowa są uporządkowane alfabetycznie features = tfidf_vec.get_feature_names_out()

Sprawdź wizualnie, że ta reprezentacja jest rzeczywiście rzadka:

import matplotlib.pyplot as plt
import pandas as pd
 
sample_matrix = tfidf_sparse_matrix[:20, :].toarray()
 
plt.figure(figsize=(8, 6))
plt.spy(tfidf_sparse_matrix[:20, :], precision=0.01, aspect='auto', markersize=1)
plt.title("Rozkład niezerowych wartości (Sparsity Pattern)")
plt.show()
 
plt.figure(figsize=(8, 8))
plt.imshow(sample_matrix, aspect='auto', cmap='viridis')
plt.colorbar(label='TF-IDF Score')
plt.xlabel("Indeks Słowa (Feature Index)")
plt.ylabel("Indeks Dokumentu")
plt.title("Fragment macierzy TF-IDF (pierwsze 20 dokumentów)")
plt.show()

Wypisz pierwszych 20 słów o najwyższych wartościach TF-IDF

weights = tfidf_sparse_matrix.sum(axis=0).A1
vocab = tfidf_vec.get_feature_names_out()

Oczekiwany wynik - podobny do poniższej tabeli:

TODO Utwórz obiekt pandas dataframe z kolumnami word_index, word oraz idf Słowa odczytujemy za pomocą funkcji tfidf_vec.get_feature_names_out(). Indeks słowa to numer na liście, a wartość idf odczytamy z atrybutu wynikowego tfidf_vec.idf_.

Przykładowy wynik:

Zakładamy że df (typu pandas.DataFrame) zawiera wynik przetwarzania, a df_vocab słownik z wartościami idf.

Wewnątrz kwerend duckdb mozna odwoływać się do obiektów pandas.DataFrame tak, jakby były widokami. Za pomocą poniższego kodu można zapisać dane w duckdb.

import duckdb
 
con = duckdb.connect('wikipedia.duckdb')
 
# Tworzymy tabelę bezpośrednio z DataFrame
con.execute("CREATE OR REPLACE TABLE wikipedia_corpus AS SELECT * FROM df")
con.execute("CREATE OR REPLACE TABLE vocabulary AS SELECT * FROM df_vocab")
 
print(con.execute("SHOW TABLES").fetchall())
con.close()

Napiszemy funkcję, do której zostanie przekazany tekst jako argument, a następnie zwróci ona wskazaną liczbę n_samples najbardziej podobnych dokumentów. Wpierw dokonywana jest zamiana tekstu na reprezentację spaCy, a następnie odczytywane wektory osadzeń i tokeny. Funkcja ma trzy tryby działania (opis w docstringu):

• vector - wykorzystuje wektory osadzeń
• token - na podstawie zbiorów słów
• hybrydowy - łączący oba

Implementacja trybu hybrydowego jest podana.

TODO Należy zaimplementować tryby: wektorowy i oparty na tokenach.

import pandas as pd
import numpy as np
 
def find_nearest_documents(con, nlp, query_text: str, n_samples: int = 10, mode: str = "vector") -> pd.DataFrame:
    """
    Funkcja wyszukuje najbardziej podobne dokumenty w bazie DuckDB na podstawie tekstu zapytania.
 
    @param con - połączenie do bazy DuckDB zawierającej tabelę 'wikipedia_corpus'
    @param nlp - zainicjalizowany pipeline SpaCy (np. pl_core_news_md lub pl_core_news_trf)
    @param query_text - tekst zapytania, dla którego szukamy podobnych dokumentów
    @param n_samples - liczba zwracanych najbliższych dokumentów
    @param mode - tryb wyszukiwania: 
                  "vector"  - oblicza odległość cosinusową między osadzeniami (embeddings)
                  "tokens"  - wyszukuje na podstawie pokrycia (części wspólnej) tokenów tekstowych
                  "hybrid"  - łączy podejście wektorowe i tokenowe przy użyciu normalizacji rangowej
    @returns pandas.DataFrame zawierający n_samples najbardziej podobnych dokumentów wraz z miarą dopasowania
    """
    # Przetwarzanie tekstu zapytania za pomocą SpaCy
    doc = nlp(query_text)
 
    # Wyciągamy wektor 
    query_vector = doc.vector.tolist()
 
    # Wyciągamy czyste tokeny (odrzucając stop-words i interpunkcję)
    query_tokens = [token.lemma_.lower() for token in doc if not token.is_stop and not token.is_punct]
 
    if mode == "vector":
        # Używamy LIST_COSINE_SIMILARITY, która przyjmuje zwykłe listy (FLOAT[])
        pass
 
    elif mode == "tokens":
        # Szukamy dokumentów, które mają najwięcej wspólnych słów z zapytaniem
        pass
 
    elif mode == "hybrid":
        # Pobieramy szerszy zestaw kandydatów z obu metod (np. 10x więcej niż n_samples)
        oversample = n_samples * 10
 
        # Pobranie danych metodą wektorową (SELECT * już zawiera wszystkie potrzebne kolumny)
        df_vec = con.execute(f"""
            SELECT *,
                   1 - LIST_COSINE_SIMILARITY(vector::FLOAT[], $1::FLOAT[]) as distance
            FROM wikipedia_corpus
            ORDER BY distance ASC
            LIMIT $2
        """, [query_vector, oversample]).df()
 
        # Pobranie danych metodą tokenową
        df_tok = con.execute(f"""
            SELECT *, 
                   LEN(ARRAY_INTERSECT(tokens, $1::VARCHAR[])) as score 
            FROM wikipedia_corpus 
            ORDER BY score DESC 
            LIMIT $2
        """, [query_tokens, oversample]).df()
 
        # Przypisanie rang (pozycji) w obu rankingach
        df_vec['rank_vec'] = df_vec.index + 1
        df_tok['rank_tok'] = df_tok.index + 1
 
        # Łączenie zestawów danych na podstawie kolumny 'title'
        hybrid_df = pd.merge(df_vec[['title', 'rank_vec']], df_tok[['title', 'rank_tok']], on='title', how='outer')
        hybrid_df.fillna(oversample + 1, inplace=True)
 
        # Obliczanie wyniku hybrydowego (Reciprocal Rank Fusion)
        hybrid_df['hybrid_score'] = (1 / (60 + hybrid_df['rank_vec'])) + (1 / (60 + hybrid_df['rank_tok']))
        hybrid_df = hybrid_df.sort_values(by='hybrid_score', ascending=False).head(n_samples)
 
        # Wyciągnięcie listy zwycięskich tytułów
        target_titles = hybrid_df['title'].tolist()
 
        # 7. Szybkie odzyskanie pełnych danych bez ponownego odpytywania bazy SQL:
        # Łączymy wyniki z df_vec (lub df_tok), odfiltrowując je tylko do najlepszych tytułów.
        # Używamy kombinacji z df_vec i df_tok na wypadek, gdyby jakiś dokument był tylko w jednym z nich.
        combined_all_data = pd.concat([df_vec, df_tok]).drop_duplicates(subset=['title'])
        final_df = combined_all_data[combined_all_data['title'].isin(target_titles)].copy()
 
        # Opcjonalnie: sortujemy finalny DataFrame dokładnie w takiej kolejności, jaką wyznaczył algorytm hybrydowy
        final_df['title'] = pd.Categorical(final_df['title'], categories=target_titles, ordered=True)
        return final_df.sort_values('title').reset_index(drop=True)        
    else:
        raise ValueError(f"Nieznany tryb wyszukiwania: {mode}. Wybierz 'vector', 'tokens' lub 'hybrid'.")

W trybie “hybrid” stosowany jest algorytm Reciprocal Rank Fusion (RRF). Jest on często używany w systemach wyszukiwania informacji (Information Retrieval) do łączenia wyników z kilku różnych wyszukiwarek.

Dla każdego dokumentu ze zbioru wszystkich znalezionych dokumentów, oblicza sumaryczną ocenę według poniższego wzoru: $$ RRF\_Score(d \in D) = \sum_{m \in M} \frac{1}{k + r_m(d)} $$ Gdzie:

• $M$ – zestaw systemów wyszukiwania (np. $m_1$ = wyszukiwanie wektorowe, $m_2$ = wyszukiwanie tokenowe),
• $r_m(d)$ – ranga (pozycja) dokumentu $d$ na liście zwróconej przez system $m$ (pierwsze miejsce to $1$, drugie to $2$, itd.)
• $k$ – stała wygładzająca (stała regularyzacyjna). Standardowo w literaturze i systemach takich jak Elasticsearch czy Qdrant przyjmuje się $k = 60$

if 'nlp' not in globals():
    import spacy
    nlp = spacy.load("pl_core_news_lg")
 
import duckdb
 
def check_if_active(connection_var_name):
    if connection_var_name not in globals():
        return False
    try:
        globals()[connection_var_name].execute("SELECT 1;")
        return True
    except:
        return False
 
if not check_if_active("con"):
    con = duckdb.connect("wikipedia.duckdb")
 
query_text = "Kraków to miasto połozone nad rzeką Wisłą w południowej Polsce. Zostało założone przez księcia Kraka w IX wieku. "
df_result  = find_nearest_documents(con, nlp, query_text, mode="hybrid")
df_selected = df_result[["title","text"]]
df_selected['text'] = df_selected.text.apply(lambda r: r[:80])
print(df_selected.to_markdown(index=False))

Przykładowy wynik:

| title                | text                                                                             |
|:---------------------|:---------------------------------------------------------------------------------|
| Grody Czerwieńskie   | Grody Czerwieńskie, Ziemia czerwieńska (, ) – przyjęta w historiografii nazwa zi |
| Goci                 | Goci, Gotowie (, Gutþiuda; , u Pliniusza Starszego Gutones, u Tacyta Gotones; )  |
| Antwerpia            | Antwerpia (, wym. []; , wym. []) – miasto w północnej Belgii, w Regionie Flamand |
| Przemyśl             | Przemyśl (, , ) – miasto na prawach powiatu w południowo-wschodniej Polsce, w wo |
| Baranów Sandomierski | Baranów Sandomierski () – miasto w Polsce, w województwie podkarpackim, w powiec |
| Wolin (wyspa)        | Wolin () – przybrzeżna wyspa należąca do Polski (powiat kamieński oraz Świnoujśc |
| Włocławek            | Włocławek (, ) – miasto w centralnej Polsce, na prawach powiatu. Trzecie co do w |
| Oświęcim             | Oświęcim (, Oszpicin, ) – miasto w województwie małopolskim, siedziba władz powi |
| Konstantynopol       | Konstantynopol () – nazwa Bizancjum nadana miastu przez Konstantyna Wielkiego, k |
| 3 maja               | Święta Imieniny obchodzą: Aleksander, Alodia, Antonina, Diodor, Diodora, Juwenal |

TODO Przygotuj 5 przykładów tekstów o różnej tematyce. Następnie znajdź 10 najbliższych dokumentów w 3 trybach. Oceń poziom dopasowania licząc $precision@10$

$$precision@10=\frac{\text{Liczba dopasowanych dokumentów na liście top 10}}{10}$$

Wyniki zbierz w tabelce.

Zastosujemy redukcję wymiarów, aby:

• wyświetlać punkty reprezentujących dokumenty na wykresach 2D
• dokonać konwersji do postaci wektorów o rozmiarach 10 i 20 w celu dalszego przeprowadzenia na nich grupowania

UMAP (Uniform Manifold Approximation and Projection) to metoda redukcji wymiarowości oparta na założeniu, że dane wysokowymiarowe leżą na rozmaitości (manifold) o niższym wymiarze, którą można lokalnie aproksymować. Algorytm buduje ważony graf sąsiedztwa w przestrzeni wejściowej, wykorzystując najbliższych sąsiadów.

Dla punktów $x_i$ i $x_j$ prawdopodobieństwo sąsiedztwa jest definiowane jako: $$ p_{ij} = \exp\left(-\frac{\max(0, d(x_i,x_j)-\rho_i)}{\sigma_i}\right) $$

gdzie $\rho_i$ odpowiada odległości do najbliższego sąsiada zapewniającej lokalną spójność, a $\sigma_i$ jest parametrem normalizującym dobieranym tak, aby liczba efektywnych sąsiadów była zgodna z parametrem n_neighbors.

W przestrzeni niskowymiarowej UMAP konstruuje analogiczny graf z prawdopodobieństwami: $$ q_{ij} = \frac{1}{1 + a \|y_i - y_j\|^{2b}} $$ gdzie $y_i$ i $y_j$ są reprezentacjami punktów po redukcji wymiaru, a parametry a,b kontrolują kształt funkcji odległości.

Algorytm minimalizuje funkcję celu będącą krzyżową entropią pomiędzy strukturą wysokowymiarową i niskowymiarową:

$$ \mathcal{L} = \sum_{i \neq j} \left[ p_{ij}\log\frac{p_{ij}}{q_{ij}} + (1-p_{ij})\log\frac{1-p_{ij}}{1-q_{ij}} \right] $$

Minimalizacja tej funkcji powoduje, że punkty będące sąsiadami w przestrzeni wejściowej pozostają blisko siebie po projekcji, przy jednoczesnym odpychaniu punktów niesąsiednich.

Wykonaj poniższy kod (po załadowaniu danych do dataframe df.

Wydaje się, że algorytm pozostawia dużo drobnoziarnistych obiektów w pamięci, dlatego po jego wywołaniu został dodany kod do wymuszenia grabage collection.

import umap
reducer = umap.UMAP(n_neighbors=30, n_components=2, metric='cosine', random_state=42)
umap_results = reducer.fit_transform(list(df['vector']))
 
df['umap_x'] = umap_results[:, 0]
df['umap_y'] = umap_results[:, 1]
 
# garbage collection
import gc
del reducer
del umap_results
gc.collect()

W podobny sposób przekształć wektory w kolumnie “vector_sparse”.

Wynik można zwizualizować za pomocą plotly:

import plotly.express as px
 
# Szybki podgląd dla Sparse
fig = px.scatter(df, x='umap_x', y='umap_y', hover_name='title', title='UMAP na SpaCy (Dense)')
fig.show()

Napisz funkcję według specyfikacji. Możesz w niej uwaględnić także zwolnienie obiektów i uruchomienie garbage colleection.

def add_low_dim_features(df, n_comp, source_column, features_columns):
    """
    Przeprowadza redukcję wymiarów i dodaje wynikową kolumnę do dataframe. 
    @param df - pandas.DataFrame zawierająca kolumnę z wektorami poddanymi redukjci wymiarów
    @param n_comp - wymiary wektorów po redukcji
    @param source_column - kolumna zawierająca wektory wejściowe (tablice float)
    @param features_columns - kolumna do której zostaną umieszczone wynikowe wektory (w postaci list)
    """
    print(f"Generuję UMAP (wymiar = {n_comp}) dla {source_column}...")

oraz wywołaj ją dodając kilka wariantów cech:

add_low_dim_features(df,10, 'vector','umap_10_dense')
add_low_dim_features(df,20, 'vector','umap_20_dense')
add_low_dim_features(df,10, 'vector_sparse','umap_10_sparse')
add_low_dim_features(df,20, 'vector_sparse','umap_20_sparse')
add_low_dim_features(df,2, 'vector','umap_2_dense')
add_low_dim_features(df,2, 'vector_sparse','umap_2_sparse')

Zapisz wynik do duckdb :CREATE OR REPLACE TABLE wikipedia_corpus AS SELECT * FROM df

Tabela (dataframe lub duckdb) zawiera informacje o dokumentach:

• tytuły
• treść
• tokeny
• oryginalne cechy - wektory (spaCy i TF-IDF)
• cechy UMAP (wektory o niższych wymiarach)
• wektory używane do wyświetlania

Do tego dodamy kolumny z etykietami grupowania przeprowadzonego za pomoca różnie skonfigurowanych algorytmów z uzyciem różnych cech. Przyjmijmy, że nazwa kolumny z grupowaniem będzie kodowała interesujące nas metadane. Na przykład będzie miała postać:

cluster.{algorithm}.{algorithm_params}.{features_column}

Czyli, np. nazwa kolumny “cluster.kmeans.10.vector” koduje informacje, że zawiera ona etykiety grupowania za pomocą algorytmu k-means z parametrem k=10 na cechach vector (czyli pełnych wektorach ze SpaCy. Z kolei “cluster.kmeans.4.umap_20_dense” to wynik k-means dla k=4 na cechach umap_20_dense.

Zwróć uwagę na obecność kropek w nazwach kolumn. W takim przypadku nazwa kolumny w kwerendach SQL musi być zapisana w podwójnych cudzysłowach.

Dodanie kolumny z grupowaniem

import duckdb
import pandas as pd
import numpy as np
from sklearn.cluster import KMeans
import plotly.express as px
 
con = duckdb.connect('wikipedia.duckdb')
df_tmp = con.execute("SELECT * FROM wikipedia_corpus").df()
con.close()
 
n_clusters = 10
features_column = 'umap_20_dense'
X = np.stack(df_tmp[features_column].values)
kmeans = KMeans(n_clusters=n_clusters, random_state=42, n_init=10)
cluster_column = f'cluster.kmeans.{n_clusters}.{features_column}'
df_tmp[cluster_column] = kmeans.fit_predict(X).astype(str)

Wizualizacja:

fig = px.scatter(
    df_tmp, 
    x='umap_x', 
    y='umap_y', 
    color=str(cluster_column),  # Kolorowanie według klastrów
    hover_name='title',        # Tytuł artykułu po najechaniu
    title=f'K-Means (k={n_clusters}) na {features_column}',
    labels={cluster_column: 'Klaster'},
    template='plotly_white',
    width=1000,
    height=700
)
 
fig.update_traces(marker=dict(size=4, opacity=0.7))
 
fig.show()

Przy wizualizacji za pomocą plotly trzeba ustawić color=str(cluster_column). Wtedy numery klastrów zostaną potraktowane jak wartości kategoryczne i kolory zostaną wybrane z mapy o dyskretnych kolorach. Dla int byłyby to mało wyraziste kolory z mapy ciagłej.

Zaimplementuj funkcję według specyfikacji w docstring

def add_kmeans_clustering(df, n_clusters, features_column, cluster_column):
    """
    Funkcja wywołuje algorytm k-means na wektorach features i dodaje do df kolumnę cluster_column
    z numerami grup
    :param df: wejściowy zbiór danych (pandas DataFrame)
    :param n_clusters: liczba grup
    :param features_column: nazwa kolumny z wektorami obserwacji
    :param cluster_column: nazwa kolumny z numerami grup (zakkodowanymi jako stringi
    :return: df 
    """

Wykorzystaj poniższą funkcję do wyświetlania danych. Oblicza ona centroidy oraz rysuje diagram Woronoja. Zwróć uwagę na to, że grupowanie odbywa się w przestrzeni o większych wymiarach niż 2D, więc grupy mogą nie mieścić się wewnątrz komórek Woronoja 2D.

plot_kmeans_results.py

import numpy as np
import plotly.express as px
import plotly.graph_objects as go
from scipy.spatial import Voronoi
 
 
def plot_kmeans_results(df, cluster_column, title, x_column='umap_x', y_column='umap_y', show_voronoi=False):
 
    x_min, x_max = df[x_column].min(), df[x_column].max()
    y_min, y_max = df[y_column].min(), df[y_column].max()
 
    margin_x = (x_max - x_min) * 0.05
    margin_y = (y_max - y_min) * 0.05
 
    fig = px.scatter(
        df, x=x_column, y=y_column, 
        color=df[cluster_column].astype(str),
        labels={cluster_column: "Klaster"},
        hover_name='title', title=title, template='plotly_white',
        width=1000, height=700
    )
    fig.update_traces(marker=dict(size=4, opacity=0.6))
 
 
    if show_voronoi:
        centroids_2d = df.groupby(cluster_column)[[x_column, y_column]].mean().values
        vor = Voronoi(centroids_2d)
        center = centroids_2d.mean(axis=0)
 
        far_dist = (x_max - x_min) * 5 
 
        for pointidx, simplex in zip(vor.ridge_points, vor.ridge_vertices):
            simplex = np.asarray(simplex)
            if np.all(simplex >= 0):
                v1, v2 = vor.vertices[simplex]
                fig.add_shape(type="line", x0=v1[0], y0=v1[1], x1=v2[0], y1=v2[1],
                              line=dict(color="rgba(0,0,0,0.8)", width=1, dash="dot"))
            else:
                i = simplex[simplex >= 0][0]
                t = centroids_2d[pointidx[1]] - centroids_2d[pointidx[0]]
                t /= np.linalg.norm(t)
                n = np.array([-t[1], t[0]])
                midpoint = centroids_2d[pointidx].mean(axis=0)
                direction = np.sign(np.dot(midpoint - center, n)) * n
                far_point = vor.vertices[i] + direction * far_dist
 
                fig.add_shape(type="line", x0=vor.vertices[i][0], y0=vor.vertices[i][1],
                              x1=far_point[0], y1=far_point[1],
                              line=dict(color="rgba(0,0,0,0.8)", width=1, dash="dot"))
 
 
 
                fig.add_trace(go.Scatter(
                    x=centroids_2d[:, 0], 
                    y=centroids_2d[:, 1],
                    mode='markers', 
                    marker=dict(
                        symbol='x', 
                        size=15, 
                        color='white',      # Białe tło (obrys)
                        line=dict(width=4, color='white') 
                    ),
                    showlegend=False,
                    hoverinfo='skip'
                ))
 
                fig.add_trace(go.Scatter(
                    x=centroids_2d[:, 0], 
                    y=centroids_2d[:, 1],
                    mode='markers', 
                    marker=dict(
                        symbol='x', 
                        size=14, 
                        color='black',      
                        line=dict(width=2) 
                    ),
                    name='Centroidy',
                    showlegend=False,
                    hoverinfo='skip'
                ))
 
                fig.update_layout(
                    xaxis=dict(range=[x_min - margin_x, x_max + margin_x]),
                    yaxis=dict(range=[y_min - margin_y, y_max + margin_y]),
                    # To wymusza, by scatter był pod spodem warstw rysowanych później
                    scattermode='group' 
                )
 
    fig.update_layout(
        xaxis=dict(range=[x_min - margin_x, x_max + margin_x]),
        yaxis=dict(range=[y_min - margin_y, y_max + margin_y])
    )
 
    fig.show()

Typową metodą doboru jest metoda łokcia. Zaimplementuj i wywołaj poniższa funkcję

def display_kmeans_elbow_plot(df, features_column, k_range=range(1,20) ):
    """
    Funkcja wywołuje algorytm k-means na wektorach features i tworzy listę wartości funkcji celu
    (atrybut kmeans.inertia_ obiektu kmeans = KMeans(...)). Następnie sporządza wykres zalezności inertia od k.
 
    :param df: wejściowy zbiór danych (pandas DataFrame)
    :param k_range: lista wartości
    :param features_column: nazwa kolumny z wektorami obserwacji
    :param cluster_column: nazwa kolumny z numerami grup (zakkodowanymi jako stringi
 
    """

Wartość punktu załamania można też próbowac wyznaczyć automatycznie na podstawie odległości od linii łączącej krańcowe punkty:

find_elbow.py

def find_elbow(k_range, inertia):
    # Konwersja wejściowych list na tablice NumPy dla bezpieczeństwa
    k_range = np.asarray(k_range)
    inertia = np.asarray(inertia)
 
    # Współrzędne punktów skrajnych
    x1, y1 = k_range[0], inertia[0]
    x2, y2 = k_range[-1], inertia[-1]
 
    # Równanie prostej przechodzącej przez p1 i p2 w postaci ogólnej: Ax + By + C = 0
    # A = y2 - y1, B = -(x2 - x1), C = x2*y1 - y2*x1
    A = y2 - y1
    B = -(x2 - x1)
    C = x2 * y1 - y2 * x1
 
    # Obliczamy odległości dla wszystkich punktów jednocześnie (bez pętli)
    # Wzór: |A*x + B*y + C| / sqrt(A^2 + B^2)
    numerator = np.abs(A * k_range + B * inertia + C)
    denominator = np.sqrt(A**2 + B**2)
    distances = numerator / denominator
 
    # Zwracamy k, dla którego odległość jest największa
    return k_range[np.argmax(distances)]

W przypadku dokumentów tekstowych granica nie jest bardzo wyrazista. Optymalna wartość (tu 4 grupy) też nie zawsze jest użyteczna.

Liczbę grup można dobrać na podstawie oceny kształtu dla kolejnych wartości k. Wykorzystamy miarę Silhouette. Napisz funkcję według specyfikacji i wyświetl wykres.

from sklearn.metrics import silhouette_score
import matplotlib.pyplot as plt
 
def display_kmeans_silhouette_plot(df, features_column, k_range=range(2,40,2) ):
    """
    Funkcja wywołuje algorytm k-means na wektorach features i tworzy listę wartości współczynnika Silhouette 
    dla róznych wartości k. Następnie wyświetla jej wykres
 
    Parametry:
    @df - dataframe z kolumną features_column
    @features_column - kolumna z wektorami cech
    k_range - sekwencja wartości
 
    """

Dodaj grupowanie k-means dla takiej siatki (iloczynu kartezzjańskiego) parametrów

k_range = [k_z_elbow, 10, k_z_silhouette]
features = ['vector','vector_sparse','umap_10_dense','umap_10_sparse','umap_20_dense', 'umap_20_sparse']

zapisz wyniki w duckdb.

TODO Wyświetl i przeanalizuj kilka (2-3) przykłady grupowania. Interesujące są dokumenty na granicy komórek Woronoja. Oceń ich zgodność tematyczną z dokumentami w grupie.

TODO Dla wybranego grupowania wyświetl chmury słów dla poszczególnych grup. Użyj gotowej funkcji (przeanalizuj zastosowaną kwerendę). UNNEST odpowiada explode() Sparka lub flatten_map() w Ray.

def display_word_cloud(con,cluster_column,cluster_id ):
 
    # 2. Kwerenda: 
    # - UNNEST(tokens) zamienia listę ['a', 'b'] w osobne wiersze
    # - JOIN z vocabulary pobiera idf_score dla każdego słowa
    # - SUM(idf_score) daje nam wagę znaczenia słowa w danym klastrze
    query = f"""
        SELECT 
            t.word, 
            SUM(v.idf_score) as total_weight
        FROM (
            SELECT UNNEST(tokens) as word 
            FROM wikipedia_corpus 
            WHERE "{cluster_column}" = '{cluster_id}'
        ) t
        JOIN vocabulary v ON t.word = v.word
        GROUP BY t.word
        ORDER BY total_weight DESC
        LIMIT 100
    """
    query = f"""
        SELECT 
            t.word, 
            COUNT(*) as count
        FROM (
            SELECT UNNEST(tokens) as word 
            FROM wikipedia_corpus 
            WHERE "{cluster_column}" = '{cluster_id}'
        ) t
        --- JOIN vocabulary v ON t.word = v.word
        GROUP BY t.word 
        ORDER BY count DESC
        LIMIT 100
    """
 
 
    # Pobieramy wynik jako słownik {słowo: waga}
    result = con.execute(query).fetchall()
    word_weights = {row[0]: row[1] for row in result}
 
    # 3. Generowanie chmury na podstawie wag (Frequencies)
    if word_weights:
        wc = WordCloud(
            width=1000, 
            height=500, 
            background_color='black', 
            colormap='spring',
            max_font_size=150
        ).generate_from_frequencies(word_weights)
 
        # 4. Wyświetlanie
        plt.figure(figsize=(15, 8))
        plt.imshow(wc, interpolation='bilinear')
        plt.axis('off')
        plt.title(f"Tematyka klastra {cluster_id} dla grupowania {cluster_column}", fontsize=20)
        plt.show()
    else:
        print(f"Brak danych dla klastra {cluster_id}")

TODO Napisz komentarz jaka jest przypuszczalna tematyka dokumentów w grupach

Użyjemy (wygenerowanych przez gemini :) ) funkcji, których zadaniem jest wydruk charakterystycznych dokumentów w grupie. Wpierw obliczymy centroidy, a następnie pweną liczbę dokumentów położonych najblizej srodka. Obliczenia zostaną przeprowadzone za pomocą kwerend wykonywanych przez duckdb.

get_cluster_centers.py

import pandas as pd
 
def get_cluster_centers(con, features_column: str, cluster_column: str) -> pd.DataFrame:
    """ 
    Funkcja zwraca pandas DataFrame z kolumnami cluster_column i "centroid" 
    zawierającą pary (cluster_id, list of float).
    """
 
    # Kwerenda wykonuje następujące kroki:
    # 1. UNNEST z WITH ORDINALITY rozbija listę na poszczególne wartości (val) i ich indeksy (idx)
    # 2. Grupuje po klastrze i indeksie wymiaru, licząc średnią (AVG)
    # 3. Składa średnie z powrotem w listę (LIST) zachowując kolejność wymiarów (ORDER BY idx)
 
    query = f"""
        WITH exploded AS (
            SELECT 
                "{cluster_column}",
                UNNEST({features_column}) as val,
                GENERATE_SUBSCRIPTS({features_column}, 1) as idx
            FROM wikipedia_corpus
        ),
        averages AS (
            SELECT 
                "{cluster_column}",
                idx,
                AVG(val) as avg_val
            FROM exploded
            GROUP BY "{cluster_column}", idx
        )
        SELECT 
            "{cluster_column}",
            LIST(avg_val ORDER BY idx) as centroid
        FROM averages
        GROUP BY "{cluster_column}"
        ORDER BY "{cluster_column}"
    """
 
    return con.execute(query).df()

Dalej:

documents_by_cluster.py

def documents_by_cluster(features_column,cluster_column, con=None, n_documents = 20,lead_text_len=100):
    if con is None:
        con_ = duckdb.connect('wikipedia.duckdb')
    else:
        con_ =con
    df_centroids = get_cluster_centers(con_,features_column=features_column, cluster_column=cluster_column)
 
    query=f"""WITH distances AS (
        SELECT 
            s."{cluster_column}" AS cluster,
            s.title,
            -- Skracamy tekst do 100 znaków
            LEFT(s.text, {lead_text_len}) AS lead_text,
            -- Obliczamy kwadrat odległości euklidesowej między wektorem a centroidem
            sqrt(list_reduce(list_transform(
                list_zip(s."{features_column}", c.centroid), 
                x -> (x[1] - x[2]) * (x[1] - x[2])
            ), (a, b) -> a + b)) AS dist_to_centroid
        FROM wikipedia_corpus s
        JOIN df_centroids c ON s."{cluster_column}" = c."{cluster_column}" 
    ),
    ranked AS (
        SELECT *,
               ROW_NUMBER() OVER(PARTITION BY cluster ORDER BY dist_to_centroid ASC) as rank
        FROM distances
    )
    SELECT cluster, rank, title, lead_text, dist_to_centroid
    FROM ranked
    WHERE rank <= {n_documents}
    ORDER BY cluster, rank;"""
    df_docs = con_.execute(query).df()
    if con is None:
        con_.close()
    return df_docs
 
def print_documents_by_cluster(df_docs):
    for cluster_id in sorted(df_docs['cluster'].unique(), key=int):
        print(f"\n{'='*80}")
        print(f" KLASTER {cluster_id} - NAJBLIŻEJ CENTROIDU ".center(80, ' '))
        print(f"{'='*80}")
 
        # Filtrujemy dokumenty dla danego klastra
        subset = df_docs[df_docs['cluster'] == cluster_id].sort_values('rank')
 
        for _, row in subset.iterrows():
            print(f"[{row['rank']}] {row['title']}")
            print(f"    {row['lead_text']}...")
            print("-" * 40)

TODO Wybierz jedno z grupowań oraz kolumnę features. Wydrukuj dokumenty należące do klastra, a następnie spróbuj określić tematykę klastra. Przykładowy wynik: documents_by_cluster.txt.zip

Oceń tematykę grup. Wyniki zbierz w tabelce. Do analizy możesz skorzystać z wybranego narzędzia LLM (przesłać plik z odpowiednim promptem).

Napisz funkcję według specyfikacji. Użyj klasy GaussianMixture

def add_em_clustering(df, n_clusters, features_column, cluster_column, covariance_type='full'):
    """
    Funkcja wywołuje algorytm em na wektorach features i dodaje do df kolumnę cluster_column
    z numerami grup
    :param df: wejściowy zbiór danych (pandas DataFrame)
    :param n_clusters: liczba grup
    :param features_column: nazwa kolumny z wektorami obserwacji
    :param cluster_column: nazwa kolumny z numerami grup
    :param covariance_type  {'full', 'tied', 'diag', 'spherical'}
    :return: df 
    """

Dodaj grupowanie z kowariancją full dla takich samych liczb klastrów jak dla k-means, ponieważ ułatwi to porównanie. W nazwie uwzglednij kowariancję:

k_range = [???,10,??]
features = ['umap_10_dense','umap_10_sparse','umap_20_dense', 'umap_20_sparse']
 
for  ...:
   cluster_column = f'cluster.em.{k}.full.{f}' 

Dodaj też grupowanie dla różnych postaci macierzy kowariancji i np. k=10

cov_types = ['tied', 'diag', 'spherical']
k=10
features = 'umap_10_dense'
# zakoduj kowariancję w nazwie
cluster_column = f'cluster.em.{k}.{cov}.{features}' 

Użyj poniższej funkcji do wyświetlania grup oraz aproksymacji elips definiujących kształt klastrów.

plot_em_results.py

import numpy as np
import plotly.graph_objects as go
import plotly.express as px
from sklearn.mixture import GaussianMixture
 
def plot_em_results(df, cluster_column, title, x_column='umap_x', y_column='umap_y', 
                    show_ellipses=True, covariance_type='full'):
    """
    Jeśli parametr show_ellipses=True funkcja aproksymuje elipsę odpowiadającą macierzy kowariancji.
    - covariance_type: 'full' - na podstawie rzeczywistych etykiet
    - covariance_type: ''tied', 'diag', 'spherical' - buduje model GMM na danych dwuwymiarowych
    """
 
    # 2. Podstawowy wykres punktowy (z nowymi etykietami klastrów z przestrzeni 2D)
    fig = px.scatter(
        df, x=x_column, y=y_column, 
        color=df[cluster_column].astype(str),
        labels={cluster_column: "Klaster"},
        hover_name='title', title=f"{title} ({covariance_type})", template='plotly_white',
        width=1000, height=700
    )
    fig.update_traces(marker=dict(size=4, opacity=0.5))
 
    if show_ellipses:
        if covariance_type == 'full':
            clusters = sorted(df[cluster_column].unique(), key=int)
            t = np.linspace(0, 2*np.pi, 100) # Parametr do rysowania obwodu
 
            for cluster_id in clusters:
                points = df[df[cluster_column] == cluster_id][[x_column, y_column]].values
                if len(points) < 5: continue
 
                mean = points.mean(axis=0)
                cov = np.cov(points, rowvar=False)
                vals, vecs = np.linalg.eigh(cov)
 
                # Skalowanie dla 2 odchyleń standardowych (ok. 95% punktów)
                # Używamy pierwiastka z wartości własnych jako promieni
                width, height = 2 * np.sqrt(vals) 
 
                # Obliczanie punktów elipsy z uwzględnieniem obrotu
                # x = mean_x + a*cos(t)*cos(theta) - b*sin(t)*sin(theta)
                # y = mean_y + a*cos(t)*sin(theta) + b*sin(t)*cos(theta)
                theta = np.arctan2(vecs[1, 0], vecs[0, 0])
 
                x_ellipse = mean[0] + width * np.cos(t) * np.cos(theta) - height * np.sin(t) * np.sin(theta)
                y_ellipse = mean[1] + width * np.cos(t) * np.sin(theta) + height * np.sin(t) * np.cos(theta)
 
                # Dodanie elipsy jako ścieżki (go.Scatter z wypełnieniem lub linią)
                fig.add_trace(go.Scatter(
                    x=x_ellipse, y=y_ellipse,
                    mode='lines',
                    line=dict(color='rgba(0,0,0,0.4)', width=1.5, dash='dot'),
                    showlegend=False,
                    hoverinfo='skip'
                ))
 
                # Centroid
                fig.add_trace(go.Scatter(
                    x=[mean[0]], y=[mean[1]],
                    mode='markers',
                    marker=dict(symbol='x', size=10, color='black'),
                    showlegend=False, hoverinfo='skip'
                ))
        else:
            # Pobieramy dane 2D do trenowania EM
            X_2d = df[[x_column, y_column]].values
 
            # Dynamicznie sprawdzamy liczbę klastrów na podstawie unikalnych wartości w kolumnie
            # (Zakładamy, że kolumna zawiera poprawne identyfikatory, np. liczby)
            n_clusters = len(df[cluster_column].unique())
 
            # 1. Uruchamiamy Gaussian Mixture bezpośrednio na danych 2D
            gmm = GaussianMixture(n_components=n_clusters, covariance_type=covariance_type, random_state=42)
            gmm.fit_predict(X_2d)
            t = np.linspace(0, 2*np.pi, 100) # Parametr obwodu
 
            for cluster_id in range(n_clusters):
                # Pobieramy średnią dla danego klastra z modelu GMM
                mean = gmm.means_[cluster_id]
 
                # Ekstrakcja macierzy kowariancji w zależności od covariance_type
                if covariance_type == 'full':
                    cov = gmm.covariances_[cluster_id]
                elif covariance_type == 'tied':
                    cov = gmm.covariances_  # Wspólna macierz dla wszystkich
                elif covariance_type == 'diag':
                    cov = np.diag(gmm.covariances_[cluster_id]) # Tylko wariancje na przekątnej
                elif covariance_type == 'spherical':
                    cov = np.eye(2) * gmm.covariances_[cluster_id] # Taka sama wariancja w każdym kierunku
 
                # Geometria elipsy na podstawie macierzy kowariancji
                vals, vecs = np.linalg.eigh(cov)
 
                # Skalowanie dla 2 odchyleń standardowych (95% prawdopodobieństwa rozkładu)
                width, height = 2 * np.sqrt(vals) 
                theta = np.arctan2(vecs[1, 0], vecs[0, 0])
 
                x_ellipse = mean[0] + width * np.cos(t) * np.cos(theta) - height * np.sin(t) * np.sin(theta)
                y_ellipse = mean[1] + width * np.cos(t) * np.sin(theta) + height * np.sin(t) * np.cos(theta)
 
                # Rysowanie elipsy
                fig.add_trace(go.Scatter(
                    x=x_ellipse, y=y_ellipse,
                    mode='lines',
                    line=dict(color='rgba(0,0,0,0.5)', width=1.5, dash='dot'),
                    showlegend=False, hoverinfo='skip'
                ))
 
                # Rysowanie centroidu wyznaczonego przez EM
                fig.add_trace(go.Scatter(
                    x=[mean[0]], y=[mean[1]],
                    mode='markers',
                    marker=dict(symbol='x', size=10, color='black'),
                    showlegend=False, hoverinfo='skip'
                ))
 
    # Skalowanie osi (dodaje margines, żeby elipsy nie wychodziły poza wykres)
    x_min, x_max = df[x_column].min(), df[x_column].max()
    y_min, y_max = df[y_column].min(), df[y_column].max()
    margin = 0.15
    fig.update_layout(
        xaxis=dict(range=[x_min - (x_max-x_min)*margin, x_max + (x_max-x_min)*margin]),
        yaxis=dict(range=[y_min - (y_max-y_min)*margin, y_max + (y_max-y_min)*margin])
    )
 
    fig.show()

TODO

• Wyświetl i porównaj kilka przykładów grupowania dla róznych postaci macierzy kowariancji.
• Porównaj wizualnie aproksymacje klastrów na sparse i dense

Log-Likelihood Score

W modelu mieszanin gaussowskich (GMM) gęstość prawdopodobieństwa dla pojedynczego punktu $x_i$ definiuje się jako sumę ważoną składowych gaussowskich. Metoda .score(X) w scikit-learn oblicza średni logarytm wiarygodności (log-likelihood) przypadający na jedną próbkę: $$ \text{score}(X) = \frac{1}{N} \sum_{i=1}^{N} \ln \left( \sum_{m=1}^{M} \pi_m \mathcal{N}(x_i \mid \mu_m, \Sigma_m) \right) $$

gdzie:

• $N$ to liczba próbek,
• $M$ to liczba klastrów,
• $\pi_m$ to waga danego klastra,
• $\mathcal{N}$ to wielowymiarowy rozkład normalny ze średnią $\mu_m$ i macierzą kowariancji $\Sigma_m$.

Im wyższa wartość (bliższa zeru lub dodatnia), tym lepiej model opisuje rozkład danych.

AIC (Akaike Information Criterion)

Kryterium informacyjne Akaikego ocenia jakość modelu, łacząc w ocenie jego dopasowanie i poziom skomplikowania (liczbą parametrów). Preferowane są modele dobrze dopasowane do dancyh, a równocześnie prostsze.

$$ \text{AIC} = 2k - 2\ln(L) $$

gdzie:

• $L$ to całkowity log-likelihood (czyli $\text{score}(X) \times N$),
• $k$ to łączna liczba estymowanych parametrów modelu (która rośnie wraz z liczbą klastrów oraz zależy od wybranego covariance_type).

AIC nakłada relatywnie łagodną karę za nadmiar parametrów.

BIC (Bayesian Information Criterion)

Bayesowskie kryterium informacyjne działa podobnie do AIC, lecz wprowadza znacznie surowszą karę za złożoność modelu, która rośnie logarytmicznie wraz z wielkością zbioru danych ($N$):

$$ \text{BIC} = k\ln(N) - 2\ln(L) $$

gdzie

• $k$ to liczba parametrów,
• $N$ to liczba próbek,
• $L$ to całkowity log-likelihood.

Z uwagi na kary, BIC ma tendencję do wskazywania prostszych modeli (mniejszej liczby klastrów) niż AIC.

Metody doboru

• Minimalizacja AIC/BIC
• Metoda łokcia zastosowana na AIC lub BIC lub LogLikelihood

Napisz funkcję według specyfikacji

import numpy as np
from sklearn.mixture import GaussianMixture
from tqdm import tqdm
 
def compute_gmm_scores(X, k_range, X_val=None):
    """
    Oblicza metryki dopasowania modelu GMM dla podanego zakresu liczby klastrów.
 
    Funkcja dopasowuje model GaussianMixture dla każdej wartości z zakresu `k_range`
    i zbiera kryteria informacyjne (AIC, BIC) oraz średnie wartości log-likelihood
    zarówno dla zbioru treningowego, jak i opcjonalnego zbioru walidacyjnego.
 
    Args:
        X (array-like of shape (n_samples, n_features)):
            Dane treningowe używane do dopasowania modeli GMM.
        k_range (iterable):
            Zakres lub lista liczb całkowitych określających liczbę komponentów (klastrów)
            do przetestowania (np. `range(2, 20)`).
        X_val (array-like of shape (n_samples_val, n_features), optional):
            Opcjonalne dane walidacyjne używane do obliczenia log-likelihood poza próbą treningową
            w celu wykrycia overfittingu. Domyślnie None.
 
    Returns:
        dict: Słownik zawierający listy obliczonych metryk dla każdej liczby klastrów:
            - 'n': Lista testowanych liczb komponentów (skopiowana z `k_range`).
            - 'aic': Wartości kryterium Akaikego (AIC) dla zbioru treningowego (im mniej, tym lepiej).
            - 'bic': Wartości kryterium Bayesowskiego (BIC) dla zbioru treningowego (im mniej, tym lepiej).
            - 'll': Średnie wartości log-likelihood na próbkę dla zbioru treningowego (im więcej, tym lepiej).
            - 'll_val': Średnie wartości log-likelihood na próbkę dla zbioru walidacyjnego 
              (obecne tylko wtedy, gdy przekazano `X_val`).
    """

TODO

• Wyświetl wykresy AIC i BIC (dla wybranej kolumny z cechami)
• Wyświetl wykresy Log Likelihood na zbiorze treningowym i walidacyjnym (zastsosuj podział typu 0.7:0.3)
• Wybierz liczbę grup k na podstawie miejsca minimum AIC/BIC lub maksimum/nasycenie LL. Jezeli takiego puntu nie będzie zastosuj metodę łokcia. Możesz użyc funkcji ustalającej połozenie punktu łokcia dla k-means
• Dodaj grupowanie typu full dla wyznaczonego k

Wyniki DBSCAN zależą od doboru dwóch parametrów:

• min_points - minimalnej liczby punktów w sąsiedztwie
• $\epsilon$ - promienia sąsiedztwa

W celu ich wyznaczenia tworzony jest wykres k-distance.

Zaimplementuj funkcję według specyfikacji:

def plot_multi_k_dist(X, min_pts_values=[10, 15, 20, 25, 30, 35, 40, 45, 50]):
    """
    Funkcja:
    1. Buduje reprezentację ball tree: tree = BallTree(X, metric="euclidean")
    2. Wykonuje: tree.query(X, k=?) aby pobrać informacje o odległosciach do k najbliższych sąsiadów
    3. Sortuje odległosci (uwaga, uwzględnia min_pts-1 odległości)
    4. Sporządza wykres posortowanych odległosci
 
    @param: X - macierz punktów (obserwacje w wierszach)
    @param: min_pts_values - wartosci min_pts, dla których należy sporządzić wykresy 
    """

Przykładowy wynik (z użyciem plotly)

Następnie dla wybranej liczby min_points:

• znajdź odpowiedni wykres
• odczytaj $\epsilon$ z punktu przegięcia (elbow)

Napisz funkcję według specyfikacji:

def add_DBSCAN_clustering(df, cluster_column,eps,min_samples):
    """
    Funkcja wywołuje algorytm DBSCAN na wektorach features i dodaje do df kolumnę cluster_column
    z numerami grup
    :param df: wejściowy zbiór danych (pandas DataFrame)
    :param cluster_column: nazwa kolumny z numerami grup
    :param eps - wielkość sąsiedztwa
    :param min_samples - minimalna liczba obserwacji w sąsiedztwie
    :return: df 
    """

Za jej pomocą dodaj kilka grupowań dla wybranego wektora cech features. Pamiętaj aby odpowiednio zakodować metadane w nazwie kolumny z etykietami grupowania:

cluster_column = f"cluster.dbscan.k_{p['min_samples']}.eps_{p['eps']}.{features_column}" 

Wyniki zapisz w bazie duckdb.

Wyświetlimy wykresy za pomocą poniższej funkcji. Jest ona przystosowana do wyświetlania etykiet -1 kolorem spoza mapy kolorów (jasnoszarym). DBSCAN oznacza etykietami -1 wartości odstające, bez przypisanej grupy.

def plot_results(df, cluster_column, title="", x_column='umap_x', y_column='umap_y'):
    n_clusters = len(df[cluster_column].unique()) - (1 if '-1' in df[cluster_column].values else 0)
    noise_pct = (df[cluster_column] == -1).mean() * 100
 
    fig = px.scatter(
        df, 
        x=x_column, 
        y=y_column, 
        color=df[cluster_column].astype(str),
        color_discrete_map={'-1': 'rgb(200,200,200)'}, 
        hover_name='title',        # Tytuł artykułu po najechaniu
        title=title 
            + f'({cluster_column})' 
            + f"<br><sup>Liczba klastrów: {n_clusters} | Szum: {noise_pct:.1f}%</sup>",
        labels={cluster_column: 'Klaster'},
        template='plotly_white',
        width=1000,
        height=700
    )
 
    fig.update_traces(marker=dict(size=4, opacity=0.7))
 
    fig.show()

TODO Wyświetl wykresy dla kilku wyników grupowania. Porównaj z wynikami dla k-means i EM.

HDBSCAN to wariant algorytmu DBSCAN wykorzystujący hierachiczne grupowanie aglomeracyjne.

DBSCAN zakłada, że wszystkie klastry w zbiorze mają zbliżoną gęstość zdefiniowaną sztywno za pomocą paramentów min_points i $\epsilon$. Jeśli w danych wystęopują klastry bardzo ciasne oraz klastry mocno rozproszone, DBSCAN nie wykryje ich jednocześnie – albo te rzadkie uzna za szum, albo te gęste zleje w jeden wielki konglomerat.

HDBSCAN potrafi wykrywać klastry o różnej gęstości. Nie potrzebuje parametru $\epsilon$. Zamiast tego buduje drzewo hierarchiczne metodą aglomeracyjną (od pojedynczych punktów do jednego wielkiego klastra) i automatycznie sprawdza, które zagęszczenia na różnych poziomach drzewa są stabilne i zasługują na przekształcenie w grupę.

Podstawowe parametry:

• min_cluster_size - minimalny rozmiar grupy
• min_samples - odpowiada za obliczenie (zmiennych) promieni sąsiedztwa dla indywidualnych punktów. Zamiast sztywnego $\epsilon$ dla wszystkich punktów, jego wartość obliczana jest indywidualnie. Jeżeli nie podamy min_samples, algorytm przyjmuje min_samples=min_cluster_size

Zaimplementuj funkcję według specyfikacji:

import hdbscan
 
def add_HDBSCAN_clustering(df, cluster_column,min_cluster_size,min_samples):
    """
    Funkcja wywołuje algorytm em na wektorach features i dodaje do df kolumnę cluster_column
    z numerami grup
    :param df: wejściowy zbiór danych (pandas DataFrame)
    :param cluster_column: nazwa kolumny z numerami grup
    :param min_cluster_size - minimalny rozmiar klastra
    :param min_samples - minimalna liczba obserwacji w sąsiedztwie
    :return: df 
    """

TODO Następnie użyj funkcji add_HDBSCAN_clustering do dodania kilku grupowań HDBSCAN z różnymi parametrami min_cluster_size,min_samples (ten drugi może być None) dla wybranych cech. Nie zapomnij o kodowaniu:

cluster_column = f“cluster.hdbscan.k_{p['min_samples']}.s_{p['min_cluster_size']}.{features_column}”

Wyświetl wykresy za pomocą funkcji plot_results()

TODO Porównaj rezultaty zastosowania DBSCAN i HDBSCAN

Zastosujemy metody oceny wewnętrznej wykorzystujące wartości cech ($X\in \mathbb{R}^n$) oraz etykiety grupowania.

Dla każdej kolumny z etykietami grupowania (zaczynającymi się od cluster) wyekstrahujemuy z jej nazwy:

• nazwę metody grupowania
• nazwę kolumny z cechami
• nazwę kolumny z etykietami

Oczekiwany wynik:

[{'method': 'kmeans',
  'features': 'vector',
  'clustering': 'cluster.kmeans.4.vector'},
 {'method': 'kmeans',
  'features': 'vector',
  'clustering': 'cluster.kmeans.10.vector'},
  ...

Dalej wykorzystamy podana poniżej funkcję evaluate_clusterings() do wyznaczenia trzech metryk:

• Silhouette
• Daviesa-Bouldina
• Calinskiego Harabasza

W zasadzie każda z nich wyraża oczekiwanie, że grupy powinny skupiać punkty połozone blisko ich środków (centroidów), natomiast same grupy powinny być oddalone od siebie.

evaluate_clusterings.py

from sklearn.metrics import silhouette_score, davies_bouldin_score, calinski_harabasz_score
import pandas as pd
import numpy as np
from tqdm import tqdm
 
def evaluate_clusterings(df, clustering_info):
    results = []
 
    for config in tqdm(clustering_info, desc="Obliczanie metryk"):
        method = config['method']
        feat_col = config['features']
        clust_col = config['clustering']
 
        # 1. Przygotowanie cech (X) - upewniamy się, że to macierz numpy
        if isinstance(df[feat_col].iloc[0], (list, np.ndarray)):
            X = np.stack(df[feat_col].values)
        else:
            # Jeśli to pojedyncze kolumny (np. po UMAP 2D), musisz je najpierw połączyć
            # Zakładając, że feat_col to nazwa kolumny z listą/wektorem:
            X = np.stack(df[feat_col].values)
 
        # 2. Pobranie etykiet klastrów
        labels = df[clust_col].astype(int).values
 
        # 3. Obsługa szumu (-1) dla DBSCAN/HDBSCAN
        # Metryki wymagają co najmniej 2 klastrów i nie radzą sobie dobrze z "szumem" jako klastrem
        # Częstą praktyką jest obliczanie metryk tylko dla punktów, które NIE są szumem
        mask = labels != -1
        if mask.sum() < 2 or len(np.unique(labels[mask])) < 2:
            print(f"Pominięto {clust_col}: zbyt mało klastrów po usunięciu szumu.")
            continue
 
        X_valid = X[mask]
        labels_valid = labels[mask]
 
        # 4. Obliczanie metryk
        # Silhouette: wyższa = lepiej (-1 do 1)
        s_score = silhouette_score(X_valid, labels_valid)
 
        # Davies-Bouldin: niższa = lepiej (lepiej odseparowane i zwarte klastry)
        db_score = davies_bouldin_score(X_valid, labels_valid)
 
        # Calinski-Harabasz: wyższa = lepiej (stosunek dyspersji między i wewnątrz klastrów)
        ch_score = calinski_harabasz_score(X_valid, labels_valid)
 
        results.append({
            'Grupowanie': clust_col,
            'Metoda': method,
            'Cechy': feat_col,
            'Silhouette': round(s_score, 4),
            'Davies-Bouldin': round(db_score, 4),
            'Calinski-Harabasz': round(ch_score, 2),
            'Liczba_punktów': len(labels_valid),
            'Liczba_klastrów': len(np.unique(labels_valid)),
            'Procent_szumu': round((1 - mask.mean()) * 100, 2)
        })
 
    return pd.DataFrame(results)

Rezultatem jest tablica zawierająca różne oceny:

• print(df_results.to_markdown())

TODO Utwórz listy rankingowe i wyznacz najlepsze grupowania metodą Reciprocal Rank Fusion (RRF). Uwzględnij metryki grupowania oraz szum (im wiecej tym gorzej).

Zaimplementuj funkcję add_rrf_ranking() według poniżeszej specyfikacji i wyświetl jej wyniki (np. top 10 w rankingu)

def add_rrf_ranking(df_eval: pd.DataFrame, k: int = 60) -> pd.DataFrame:
    """
    Funkcja oblicza Reciprocal Rank Fusion (RRF) na podstawie indeksów:
    Silhouette, Davies-Bouldin oraz Calinski-Harabasz oraz szumu.
 
    Użyj metody pandas.DataFrame.rank()
 
    * Silhouette: im WIĘCEJ, tym lepiej (ascending=False)
    * Davies-Bouldin: im MNIEJ, tym lepiej (ascending=True)
    * Calinski-Harabasz: im WIĘCEJ, tym lepiej (ascending=False)
    * Szum: im MNIEJ, tym lepiej (ascending=True)
    Dodaje kolumnę 'RRF_Score' oraz 'RRF_Rank' i sortuje DataFrame.
    """

Przykładowy wynik (cechy umap_2_dense niekoniecznie muszą być brane pod uwagę)

V-measure zakłada, że znamy etykiety klas i na tej podstawie oceniamy jakość grupowania. Zauważmy, że grupowanie nie jest klasyfikacją. Problemem metod oceny zewnętrznej jest permutacja etykiet grupowania względem etykiet klas.

W naszym przypadku nie mamy etykiet klas, ale różne grupowania i możemy okreslić ich podobienstwo za pomocą wartości V-measure.

plot_clustering_similarity_matrix.py

from sklearn.metrics import v_measure_score
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from tqdm import tqdm
 
def plot_clustering_similarity_matrix(df, clustering_cols):
    # 1. Inicjalizacja pustej macierzy
    n = len(clustering_cols)
    v_matrix = np.zeros((n, n))
 
    # 2. Obliczanie V-measure dla każdej pary (używa pełnych nazw kolumn)
    for i in tqdm(range(n), desc="Obliczanie V-measure"):
        for j in range(i, n):
            score = v_measure_score(
                df[clustering_cols[i]].astype(str), 
                df[clustering_cols[j]].astype(str)
            )
            v_matrix[i, j] = score
            v_matrix[j, i] = score
 
    # 3. Przygotowanie skróconych nazw poprzez usunięcie "cluster." z początku
    # Jeśli nazwa zaczyna się od 'cluster.', odcinamy to, w innym wypadku zostawiamy oryginalną.
    short_labels = [col.replace('cluster.', '', 1) if col.startswith('cluster.') else col for col in clustering_cols]
 
    # 4. Konwersja do DataFrame dla Seaborn z nowymi, skróconymi etykietami
    v_df = pd.DataFrame(v_matrix, index=short_labels, columns=short_labels)
 
    # 5. Rysowanie wykresu
    plt.figure(figsize=(22, 22)) # Delikatnie powiększone, aby pomieścić 32 kolumny
    sns.set_theme(style="white")
 
    # Tworzymy maskę, aby ukryć górny trójkąt dla czytelności
    mask = np.triu(np.ones_like(v_df, dtype=bool), k=1)
 
    heatmap = sns.heatmap(
        v_df, 
        mask=mask,
        annot=True,           # Wyświetla wartości liczbowe
        fmt=".2f",            # Formatowanie do 2 miejsc po przecinku
        cmap="YlGnBu",        # Przyjemna dla oka paleta kolorów
        linewidths=.5, 
        cbar_kws={"shrink": .7, "label": "V-measure Score"},
        annot_kws={"size": 9} # Zmniejszony font podpisów w kwadracikach, żeby się nie zbiły
    )
 
    plt.title('Podobieństwo klastrowań (V-measure Matrix)', fontsize=18, pad=25, fontweight='bold')
 
    # Obrót etykiet dla lepszej czytelności przy 32 kolumnach
    # plt.xticks(rotation=45, ha='right', fontsize=10)
    plt.xticks(rotation=90,  ha='right', fontsize=10)
    plt.yticks(rotation=0, fontsize=10)
 
    plt.tight_layout()
    plt.show()

Przykładowy wynik:

TODO Przeanalizuj macierz dla Twojego grupowania i określ, które algorytmy lub rodziny algorytmów są do siebie podobne.

TODO Napisz funkcję według specyfikacji i wyświetl top 30 podobnych grupowań

from sklearn.metrics import v_measure_score
import pandas as pd
import numpy as np
from tqdm import tqdm
 
def get_ranked_clustering_pairs(df, clustering_cols):
    """
    Funkcja oblicza podobieństwo V-measure dla każdej unikalnej pary kolumn,
    wykorzystując tqdm do wizualizacji postępu.
    Zwraca DataFrame posortowany od najbardziej do najmniej podobnych par.
    """

Przykładowy wynik:

TODO Oceń co raczej przesądza o podobieństwie:

• algorytm
• jego parametry
• cechy (rozmiary wektorów, dense vs. sparse)

Które algorytmy dają podobne rezultaty?

W notatniku można łączyć sie z bazą danych za pomocą sql magic (obecnie jupysql). Konieczny jest pakiet SQLAlchemy i odpowiedni driver bazy danych. W obrazie Dockera jest on zainstalowany.

%load_ext sql
%sql duckdb:///wikipedia.duckdb
%config SqlMagic.named_parameters = "enabled"

Następnie można wykonywać kwerendy inline

%sql select * from vocabulary limit 10

lub wieloliniowe

%%sql
-- Sprawdzamy w ilu dokumentach pojawia się "piłkarz" vs "pedagog"
SELECT 
    v.word, 
    count(t.title) as doc_count
FROM (
    SELECT title, unnest(vector_sparse) as val, generate_subscripts(vector_sparse, 1) - 1 as idx 
    FROM wikipedia_corpus
) t
JOIN vocabulary v ON t.idx = v.word_index
WHERE v.word IN ('piłkarz', 'pedagog', 'lekkoatleta', 'wykładowca') AND t.val > 0
GROUP BY v.word;

Przykład wyszukiwania wektorowego

if 'nlp' not in globals():
    import spacy
    nlp = spacy.load("pl_core_news_lg")
 
query_text = """
Tensor – obiekt matematyczny, będący – w pewien szczególny sposób określonym – uogólnieniem pojęcia wektora[a][1]. 
Zbiór wszystkich tensorów wraz z odpowiednimi działaniami dodawania i mnożenia przez skalar, 
nazywa się przestrzenią tensorową. Tensory, podobnie jak wektory, mogą być swobodne i zaczepione. 
Rozważa się pola tensorowe (nazywane również w skrócie tensorami), czyli pola, które każdemu punktowi 
przestrzeni przypisują pewien tensor. Tensory, które zmieniają się przy zmianie skali, ściśle nazywa się 
gęstościami tensorowymi.
"""
 
doc = nlp(query_text)
query_vector = doc.vector.tolist()

Wywołanie SQL

%%sql top_documents <<
SELECT 
    title, 
    CASE 
        WHEN LENGTH(text) > 200 THEN LEFT(text, 200) || '...'
        ELSE text
    END AS text_truncated,
    -- Operator <=> oblicza odległość cosinusową (im mniejsza, tym bardziej podobne dokumenty)
    (vector <=> :query_vector) AS cosine_distance
FROM 
    wikipedia_corpus
ORDER BY 
    cosine_distance ASC
LIMIT 10;

Wynik (top_documents typu ResultSet)

Operator ⇔ zwraca $$\text{cosine_distance} = 1 - \text{cosine_similarity}$$

Przykład wyszukiwania pełnotekstowego

W wyszukiwaniu pełnotekstowym (ang. Full Text Search (FTS)) korzysta się z indeksu token → dokument. Mimo jego braku w naszej bazie danych dla 10_000 dokumentów można je całkiem efektywnie wykonać…

Rozszerzenie FTS do duckdb istnieje. Można je załadować i zbudować indeks https://duckdb.org/docs/current/core_extensions/full_text_search

query_words=[t.lemma_.lower() for t in doc if not t.is_stop and t.is_alpha]

w osobnej komórce

%%sql top_documents_fts <<
SELECT 
    title, 
    CASE 
        WHEN LENGTH(text) > 200 THEN LEFT(text, 200) || '...'
        ELSE text
    END AS text_truncated,
    len(
        list_intersect(tokens, :query_words)
    ) AS token_matches
FROM 
    wikipedia_corpus
WHERE 
    token_matches > 0
ORDER BY 
    token_matches DESC
OFFSET 1
LIMIT 10;

Czasem FTS jest bardziej skuteczny i dlatego stosuje się metody hybrydowe. Jak widać, w rzypadku dokumentów wybranych na podstawie tokenów odległość cosinusowa jest znacznie większa.

%%sql top_documents_fts <<
SELECT 
    title, 
    CASE 
        WHEN LENGTH(text) > 200 THEN LEFT(text, 200) || '...'
        ELSE text
    END AS text_truncated,
    len(
        list_intersect(tokens, :query_words)
    ) AS token_matches,
    (vector <=> :query_vector) AS cosine_distance
FROM 
    wikipedia_corpus
WHERE 
    token_matches > 0
ORDER BY 
    token_matches DESC,
    cosine_distance ASC
OFFSET 1
LIMIT 10;

Aby zamknąć połaczenie

%sql --close duckdb:///wikipedia.duckdb

TODO Napisz (lub wygeneruj za pomocą LLMa) następującą kwerendę.

• W zmiennych title_a i title_b zdefiniowano tytuły dokumentów. Tytuły są unikalne w tabeli wikipedia_corpus
• W zmiennej feature_column podano nazwę kolumny z cechami. Ma postać wyrażenia regularnego, np. '^vector$'

Kwerenda znajduje w bazie wektory z kolumny feature_column . Następnie tworzy 8 punktów pośrednich na linii prostej łączącej cechy pierwszego i ostatniego puntu. Dla każdego z tych punktów znajduje dokładnie 1 najbliżej położony dokument. Ostatecznie zwraca listę tytułów znalezionych kolejno dokumentów i około 100 pierwszych znaków tych dokumentów. Lista ma zawierać początkowy i końcowy dokument.

Zaimplementu kwerendę jako %%sql (sql magic/jupysql)

Aby przekazać nazwę kolumny możesz skorzystać z poniższego przykładu:

%sql select columns(:features_column) from wikipedia_corpus where title=:title_a

Jeżeli generujesz kod za pomocą LLM - napisz jakim i w ilu iteracjach otrzymałeś ostateczny wynik

TODO Przetestuj dwa przypadki:

• punkty końcowe wewnątrz wybranego klastra
• punkty końcowe skrajnie oddalone

Dla nich porównaj trajektorie na cechach o wymiarach 300, 10 i 2. Czym się różnią i dlaczego?

Przykładowy wynik (trajektoria na umap_2_dense):