# Import libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import seaborn as sns
import zipfile, gc, time, warnings, os
from collections import defaultdict, Counter
from itertools import combinations
import lightgbm as lgb
from xgboost import XGBClassifier
from sklearn.model_selection import GroupShuffleSplit, train_test_split
from sklearn.metrics import roc_auc_score, average_precision_score
from scipy.sparse import csr_matrix
from sklearn.preprocessing import normalize
from google.colab import drive
from matplotlib.gridspec import GridSpec
from joblib import Parallel, delayed

warnings.filterwarnings('ignore')
plt.style.use('seaborn-v0_8-whitegrid')
sns.set_palette('husl')
pd.set_option('display.float_format', '{:.4f}'.format)
pd.set_option('display.max_columns', None)
print('✅  Libraries loaded')

# Read the data
drive.mount('/content/drive/')

# Define functions to compute metrics
def precision_at_k(recs, true_set, k):

    return len(set(recs[:k]) & true_set) / k if k else 0.0

def recall_at_k(recs, true_set, k):

    return len(set(recs[:k]) & true_set) / len(true_set) if true_set else 0.0

def f1_at_k(recs, true_set, k):

    p = precision_at_k(recs, true_set, k)
    r = recall_at_k(recs, true_set, k)

    return 2*p*r/(p+r) if (p+r) > 0 else 0.0

def ndcg_at_k(recs, true_set, k):

    dcg   = sum(1/np.log2(i+2) for i, m in enumerate(recs[:k]) if m in true_set)
    ideal = sum(1/np.log2(i+2) for i in range(min(len(true_set), k)))

    return dcg/ideal if ideal else 0.0

def score_recs(recs, true_set, K):

    return {
        'precision': precision_at_k(recs, true_set, K),
        'recall'   : recall_at_k(recs,    true_set, K),
        'f1'       : f1_at_k(recs,        true_set, K),
        'ndcg'     : ndcg_at_k(recs,      true_set, K),
        'hit_rate' : 1 if any(m in true_set for m in recs[:K]) else 0,
    }

def evaluate_strategy(fn, data_df, K=10, n=None, seed=42):

    rows = data_df if n is None else data_df.sample(n=n, random_state=seed)
    m = {k: [] for k in ('precision','recall','f1','ndcg','hit_rate')}

    for _, row in rows.iterrows():

        s = score_recs(fn(row['userId'], K=K), row['true_items'], K)

        for key in m: m[key].append(s[key])

    return {k: float(np.mean(v)) for k, v in m.items()}

# ── Lookup caches ─────────────────────────────────────────────────────────
# _seen_cache  : user_id -> set of seen movie IDs
# _s2_score_cache : user_id -> precomputed S2 genre score vector (n_movies,)
# _s4_score_cache : user_id -> precomputed S4 taste score vector (n_movies,)
#
# All three are initialised empty here and populated before they are needed.
# When a score cache is empty (e.g. during training pair construction, which
# uses its own batch arrays), the functions fall back to computing on the fly.
_seen_cache     = {}
_s2_score_cache = {}
_s4_score_cache = {}

def get_seen(user_id):

    return _seen_cache.get(user_id, set())

# Function to compute global popularity
def s1_popularity(user_id, K=10, **_):

    seen = get_seen(user_id)

    return [m for m in POPULARITY_POOL if m not in seen][:K]

# Function to compute genre-aware personal history
def s2_personal(user_id, K=10, **_):

    seen = get_seen(user_id)

    if user_id not in user_genre_affinity.index:

        return s1_popularity(user_id, K=K)

    # Use precomputed score vector when available (evaluation path); fall back to on-the-fly matmul otherwise (any ad-hoc call).
    if user_id in _s2_score_cache:

        scores = _s2_score_cache[user_id]

    else:

        u_genre = user_genre_affinity.loc[user_id].values
        u_norm  = np.linalg.norm(u_genre)

        if u_norm == 0:

            return s1_popularity(user_id, K=K)

        scores = movie_genre_norm @ (u_genre / u_norm)

    top_k  = np.argpartition(scores, -K)[-K:]
    top_k  = top_k[np.argsort(scores[top_k])[::-1]]
    ranked = [genre_index_order[i] for i in top_k if genre_index_order[i] not in seen]

    if len(ranked) >= K:

        return ranked[:K]

    # Pad with additional candidates if top-K after filtering is short
    all_ranked = [genre_index_order[i] for i in np.argsort(scores)[::-1] if genre_index_order[i] not in seen]

    return all_ranked[:K]

# Function to compute collaborative filtering
def s3_item_cf(user_id, K=10, **_):

    seen  = get_seen(user_id)
    liked = set(user_liked.get(user_id, []))

    if not liked:

        return s1_popularity(user_id, K=K)

    scores = defaultdict(float)

    for mid in liked:

        for neighbour, sim in item_sim_lookup.get(mid, []):

            if neighbour not in seen:

                scores[neighbour] += sim
    ranked = sorted(scores, key=scores.get, reverse=True)

    if len(ranked) >= K:

        return ranked[:K]

    extra = [m for m in s2_personal(user_id, K=K*3) if m not in set(ranked)]
    return (ranked + extra)[:K]

# Function to compute temporal taste drift
def s4_temporal(user_id, K=10, **_):

    seen = get_seen(user_id)
    # Use precomputed score vector when available (evaluation path); fall back to on-the-fly matmul otherwise.
    if user_id in _s4_score_cache:

        sims = _s4_score_cache[user_id]

    else:

        taste_vec = user_taste_vec.get(user_id)

        if taste_vec is None:

            return s1_popularity(user_id, K=K)

        sims = movie_genre_norm @ taste_vec

    top_k  = np.argpartition(sims, -K)[-K:]
    top_k  = top_k[np.argsort(sims[top_k])[::-1]]
    ranked = [genre_index_order[i] for i in top_k if genre_index_order[i] not in seen]

    if len(ranked) >= K:

        return ranked[:K]

    all_ranked = [genre_index_order[i] for i in np.argsort(sims)[::-1] if genre_index_order[i] not in seen]

    return all_ranked[:K]

# Function to build features
def build_lgb_features(user_id, candidate_movies):

    if not candidate_movies:

        return pd.DataFrame(columns = FEATURE_COLS)

    mf = movie_feat.reindex(candidate_movies)
    mg = movie_genre.reindex(candidate_movies).fillna(0).values.astype('float32')

    mg_norms = np.linalg.norm(mg, axis=1, keepdims=True).clip(min=1e-9)
    mg_norm  = mg / mg_norms

    taste_v        = user_taste_vec.get(user_id, np.zeros(len(genre_cols), dtype='float32'))
    genre_affinity = mg_norm @ taste_v

    u_row = user_stats.loc[user_id] if user_id in user_stats.index else None

    return pd.DataFrame({'userId'            : user_id,
                         'movieId'           : candidate_movies,
                         'u_rating_count'    : float(u_row['rating_count'])   if u_row is not None else 0.0,
                         'u_mean_rating'     : float(u_row['mean_rating'])    if u_row is not None else GLOBAL_MEAN,
                         'u_rating_freq'     : float(u_row['rating_freq'])    if u_row is not None else 0.0,
                         'u_rating_std'      : float(u_row['rating_std'])     if u_row is not None else 0.0,
                         'm_global_count'    : mf['global_count'].fillna(0).values,
                         'm_log_count'       : mf['log_count'].fillna(0).values,
                         'm_global_mean'     : mf['global_mean'].fillna(GLOBAL_MEAN).values,
                         'm_bayesian'        : mf['bayesian_score'].fillna(GLOBAL_MEAN).values,
                         'm_rating_variance' : mf['rating_variance'].fillna(0).values,
                         'm_tag_count'       : mf['tag_count'].fillna(0).values,
                         'genre_affinity'    : genre_affinity,})

# Function to perform lgb recommendations
def s5_lgbm(user_id, K = 10, n_candidates = 200, **_):

    seen = get_seen(user_id)
    candidates = list(dict.fromkeys(s2_personal(user_id, K = n_candidates) + s3_item_cf( user_id, K = n_candidates) + s4_temporal(user_id, K=n_candidates)))

    candidates = [m for m in candidates if m not in seen][:n_candidates]

    if not candidates:

        return s1_popularity(user_id, K=K)

    feat_df = build_lgb_features(user_id, candidates)
    feat_df['score'] = lgb_model.predict(feat_df[FEATURE_COLS])

    return feat_df.sort_values('score', ascending = False)['movieId'].tolist()[:K]

# Function to build second-stage features
def build_stage2_features(user_id, candidates, s2_list, s3_list, s4_list):

    s2r = {m: r for r, m in enumerate(s2_list)}
    s3r = {m: r for r, m in enumerate(s3_list)}
    s4r = {m: r for r, m in enumerate(s4_list)}

    base = build_lgb_features(user_id, candidates)

    base['in_s2']       = base['movieId'].isin(s2r).astype('int8')
    base['in_s3']       = base['movieId'].isin(s3r).astype('int8')
    base['in_s4']       = base['movieId'].isin(s4r).astype('int8')
    base['rank_s2']     = base['movieId'].map(s2r).fillna(N_STAGE1).astype('int16')
    base['rank_s3']     = base['movieId'].map(s3r).fillna(N_STAGE1).astype('int16')
    base['rank_s4']     = base['movieId'].map(s4r).fillna(N_STAGE1).astype('int16')
    base['n_retrievers']= base[['in_s2','in_s3','in_s4']].sum(axis=1).astype('int8')
    base['s5_score']    = lgb_model.predict(base[FEATURE_COLS]).astype('float32')

    return base

# Create a check to see if there is user history
def is_cold(user_id):

    return len(user_liked.get(user_id, [])) < COLD_START_THRESHOLD

# Function to get stage 1 candidates
def stage1_candidates(user_id):

    seen = get_seen(user_id)
    pool = list(dict.fromkeys(s2_personal(user_id, K=N_STAGE1) + s3_item_cf( user_id, K=N_STAGE1) + s4_temporal(user_id, K = N_STAGE1)))

    return [m for m in pool if m not in seen]

# Function to get ensemble recommendations
def s7_ensemble(user_id, K = 10, model = 'lgb', **_):

    if is_cold(user_id):

        return s1_popularity(user_id, K = K)

    s2l = s2_personal(user_id, K=N_STAGE1)
    s3l = s3_item_cf( user_id, K=N_STAGE1)
    s4l = s4_temporal(user_id, K=N_STAGE1)
    candidates = list(dict.fromkeys(s2l + s3l + s4l))

    if not candidates:

        return s1_popularity(user_id, K = K)

    feat = build_stage2_features(user_id, candidates, s2l, s3l, s4l)

    if model == 'lgb':

        feat['score'] = meta_lgb.predict(feat[STAGE2_FEATURE_COLS])

    else:

        feat['score'] = meta_xgb.predict_proba(feat[STAGE2_FEATURE_COLS])[:, 1]

    return feat.sort_values('score', ascending=False)['movieId'].tolist()[:K]

def s6_xgboost(user_id, K = 10, **_):

    return s7_ensemble(user_id, K = K, model = 'xgb')

# Load the data
data = {}

files = [i for i in  os.listdir('/content/drive/MyDrive/ml-32m/')]
file_list = [f for f in files if f.endswith('.csv')]
print(f'Files in archive: {file_list}\n')

for fn in file_list:

    key = fn.replace('.csv', '')
    data[key] = pd.read_csv(f"/content/drive/MyDrive/ml-32m/{fn}")
    print(f'  ✓ {fn}: {data[key].shape[0]:,} rows × {data[key].shape[1]} cols')

print('\n✅  All datasets loaded')

# Downcast to save ~40% RAM
for col in ['userId', 'movieId']:

    data['ratings'][col] = data['ratings'][col].astype('int32')
    data['tags'][col]    = data['tags'][col].astype('int32')

data['movies']['movieId'] = data['movies']['movieId'].astype('int32')
data['ratings']['rating'] = data['ratings']['rating'].astype('float32')

# Extract the data into objects
ratings = data['ratings']
movies = data['movies']
tags = data['tags']

movies.head()

tags.head()

ratings.head()

%%time

# Create a temporal training split for modelling
RATING_THRESHOLD = 4.0

# Cutoff users based on time interval, 90%
user_cutoff = (ratings.groupby('userId')['timestamp'].quantile(0.90).rename('cutoff').astype('int32'))
ratings = ratings.join(user_cutoff, on = 'userId')

# Split the data for training
train_ratings = ratings[ratings['timestamp'] <  ratings['cutoff']]
test_ratings  = ratings[ratings['timestamp'] >= ratings['cutoff']]

# Drop cutoff column immediately - saves ~100 MB across both halves
train_ratings.drop(columns = ['cutoff'], inplace = True)
test_ratings.drop(columns = ['cutoff'], inplace = True)

del ratings, user_cutoff; gc.collect()

# Ground truth: movies rated >= threshold in test, not seen in training
train_seen = train_ratings.groupby('userId')['movieId'].apply(set)

# Function to extract the truth
def get_true_items(grp):

    uid   = grp.name
    liked = set(grp[grp['rating'] >= RATING_THRESHOLD]['movieId'])

    return liked - train_seen.get(uid, set())

# Apply func to get data
ground_truth = (test_ratings.groupby('userId').apply(get_true_items).rename('true_movies'))
ground_truth = ground_truth[ground_truth.map(len) > 0]

# Inspect
print(f'Train rows : {len(train_ratings):,}  | {train_ratings.memory_usage(deep = True).sum()/1e6:.0f} MB')
print(f'Test rows  : {len(test_ratings):,}   | {test_ratings.memory_usage(deep = True).sum()/1e6:.0f} MB')
print(f'Eval users : {len(ground_truth):,}')

ground_truth.head()

# Create a global popularity pool to use as a cold-start ratings guide
pop_stats = (train_ratings[train_ratings['rating'] >= RATING_THRESHOLD].groupby('movieId').agg(pos_count = ('rating','count'), mean_rating = ('rating','mean')).sort_values(['pos_count','mean_rating'], ascending = False))
POPULARITY_POOL = pop_stats.index.tolist()
print(f'Popularity pool: {len(POPULARITY_POOL):,} movies')
pop_stats.head(10).join(movies.set_index('movieId')['title'])

# Create a global popularity pool
GLOBAL_MEAN = float(train_ratings['rating'].mean())

pop_stats = (train_ratings[train_ratings['rating'] >= RATING_THRESHOLD].groupby('movieId').agg(pos_count = ('rating','count'), mean_rating = ('rating','mean')).sort_values(['pos_count','mean_rating'], ascending = False))
POPULARITY_POOL = pop_stats.index.tolist()

print(f'Global mean       : {GLOBAL_MEAN:.4f}')
print(f'Popularity pool   : {len(POPULARITY_POOL):,} movies')
pop_stats.head(8).join(movies.set_index('movieId')['title'])

# Combine train + test for full-dataset EDA stats
n_tr, n_te    = len(train_ratings), len(test_ratings)
total_ratings = n_tr + n_te

# Union of user/movie IDs across splits without creating a combined DataFrame
users_tr  = set(train_ratings['userId'].unique())
users_te  = set(test_ratings['userId'].unique())
movies_tr = set(train_ratings['movieId'].unique())
movies_te = set(test_ratings['movieId'].unique())

n_users_total  = len(users_tr | users_te)
n_movies_rated = len(movies_tr | movies_te)
del users_tr, users_te, movies_tr, movies_te; gc.collect()

global_sum = float(train_ratings['rating'].sum()) + float(test_ratings['rating'].sum())
global_mean_all = global_sum / total_ratings
ts_min = min(int(train_ratings['timestamp'].min()), int(test_ratings['timestamp'].min()))
ts_max = max(int(train_ratings['timestamp'].max()), int(test_ratings['timestamp'].max()))
r_min  = min(float(train_ratings['rating'].min()), float(test_ratings['rating'].min()))
r_max  = max(float(train_ratings['rating'].max()), float(test_ratings['rating'].max()))
sparsity = 1 - total_ratings / (n_users_total * n_movies_rated)

print('=' * 55)
print('  MOVIELENS 32M  -  DATASET SNAPSHOT')
print('=' * 55)
print(f'  Total ratings      : {total_ratings:>12,}')
print(f'  Unique users       : {n_users_total:>12,}')
print(f'  Unique movies rated: {n_movies_rated:>12,}')
print(f'  Movies in catalogue: {len(movies):>12,}')
print(f'  Date range         : {pd.to_datetime(ts_min,unit="s").date()} -> {pd.to_datetime(ts_max,unit="s").date()}')
print(f'  Rating scale       : {r_min} - {r_max} stars')
print(f'  Mean rating        : {global_mean_all:.3f}')
print(f'  Matrix sparsity    : {sparsity:.4%}')
print('=' * 55)

# Combine data for vis
all_ratings_eda = pd.concat([train_ratings.assign(split='train'),test_ratings.assign(split='test')], ignore_index = True)
all_ratings_eda['date'] = pd.to_datetime(all_ratings_eda['timestamp'], unit='s')
all_ratings_eda['year'] = all_ratings_eda['date'].dt.year
all_ratings_eda['year_month'] = all_ratings_eda['date'].dt.to_period('M')
print(f'all_ratings_eda: {all_ratings_eda.shape}  |  '
      f'{all_ratings_eda.memory_usage(deep=True).sum()/1e6:.0f} MB')

# Ratings distribution
fig, axes = plt.subplots(1, 2, figsize=(20, 5))

# Left: count per star value
ax = axes[0]
rating_counts = all_ratings_eda['rating'].value_counts().sort_index()
bars = ax.bar(rating_counts.index, rating_counts.values,
              width=0.35, color=sns.color_palette('husl', len(rating_counts)))
ax.axvline(4.0, color='crimson', lw=2, ls='--', label='Liked threshold (4.0)')
ax.set_xlabel('Star Rating')
ax.set_ylabel('Number of Ratings')
ax.set_title('Rating Distribution (raw counts)')
ax.legend()
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'{x/1e6:.1f}M'))
for bar in bars:
    h = bar.get_height()
    ax.text(bar.get_x()+bar.get_width()/2, h*1.01,
            f'{h/1e6:.1f}M', ha='center', va='bottom', fontsize=7)

# Right: positive label rate vs threshold
ax2 = axes[1]
thresholds = np.arange(0.5, 5.5, 0.5)
pct_liked  = [(all_ratings_eda['rating'] >= t).mean() * 100 for t in thresholds]
ax2.plot(thresholds, pct_liked, marker='o', lw=2.5, color='steelblue')
ax2.axvline(4.0, color='crimson', lw=2, ls='--', label='Our threshold (4.0)')
chosen_pct = pct_liked[int((4.0-0.5)/0.5)]
ax2.axhline(chosen_pct, color='crimson', lw=1, ls=':')
ax2.annotate(f'{chosen_pct:.1f}% liked at 4.0',
             xy = (4.0, chosen_pct),
             xytext = (2.5, chosen_pct-16),
             arrowprops = dict(arrowstyle='->', color = 'crimson'),
             color = 'crimson',
             fontsize = 10)

ax2.set_xlabel('Threshold')
ax2.set_ylabel('% of Ratings Classified as Liked')
ax2.set_title('Positive Label Rate vs Threshold')
ax2.legend()
ax2.set_ylim(0, 105)

plt.suptitle('Rating Scale Analysis', fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig('eda_01_rating_distribution.png', dpi=150, bbox_inches='tight')
plt.show()
print(f'Ratings >= 4.0 stars: {(all_ratings_eda["rating"]>=4.0).mean():.1%} of all ratings')

# User actvitity distribution
user_rating_counts = all_ratings_eda.groupby('userId')['rating'].count()

fig, axes = plt.subplots(1, 3, figsize=(20, 4))

# Histogram (log y)
ax = axes[0]
ax.hist(user_rating_counts, bins=80, color='steelblue', edgecolor='white', lw=0.3)
ax.set_yscale('log')
ax.set_xlabel('Ratings per User')
ax.set_ylabel('Number of Users (log)')
ax.set_title('User Activity Distribution')

for p in [10, 50, 90]:

    v = np.percentile(user_rating_counts, p)
    ax.axvline(v, ls='--', lw=1.2, label=f'p{p}={v:.0f}')

ax.legend(fontsize=8)

# CCDF — power-law tail
ax2 = axes[1]
sorted_counts = np.sort(user_rating_counts)[::-1]
rank = np.arange(1, len(sorted_counts)+1)
ax2.loglog(sorted_counts, rank / len(sorted_counts), lw=2, color='darkorange')
ax2.set_xlabel('Ratings per User (log)')
ax2.set_ylabel('P(X >= x)  [log]')
ax2.set_title('CCDF — Heavy Tail in User Activity')
ax2.grid(True, which='both', alpha=0.3)

# Activity segments
ax3 = axes[2]
bins_   = [0, 20, 50, 100, 250, 500, int(user_rating_counts.max())+1]
labels_ = ['<20', '20-50', '51-100', '101-250', '251-500', '500+']
seg    = pd.cut(user_rating_counts, bins=bins_, labels=labels_)
seg_ct = seg.value_counts().sort_index()
bars3  = ax3.bar(range(len(seg_ct)), seg_ct.values,
                 color=sns.color_palette('husl', len(seg_ct)),
                 tick_label=labels_)
ax3.set_xlabel('Ratings per User')
ax3.set_ylabel('Number of Users')
ax3.set_title('User Activity Segments')

for i, v in enumerate(seg_ct.values):

    ax3.text(i, v + 200, f'{v:,}', ha='center', va='bottom', fontsize=8)

plt.suptitle('User Activity Analysis', fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig('eda_02_user_activity.png', dpi=150, bbox_inches='tight')
plt.show()

print(f'Median ratings/user  : {user_rating_counts.median():.0f}')
print(f'Mean   ratings/user  : {user_rating_counts.mean():.0f}')
top1 = user_rating_counts.nlargest(int(len(user_rating_counts)*0.01)).sum()
print(f'Top 1% of users account for {top1/user_rating_counts.sum():.1%} of all ratings')

# Investigate movie popularities
movie_rating_counts = (all_ratings_eda.groupby('movieId')['rating'].count().sort_values(ascending=False))
movie_pop_named = (movie_rating_counts.to_frame('n_ratings').join(movies.set_index('movieId')['title']))

fig, axes = plt.subplots(1, 2, figsize=(20, 5))

# Log-log rank curve
ax = axes[0]
ranks = np.arange(1, len(movie_rating_counts)+1)
ax.loglog(ranks, movie_rating_counts.values, lw=2, color='steelblue')
ax.fill_between(ranks, movie_rating_counts.values, alpha=0.15, color='steelblue')
top1pct = int(len(movie_rating_counts)*0.01)
ax.axvline(top1pct, color='crimson', ls='--', lw=1.5,
           label=f'Top 1% = {top1pct:,} movies')
ax.set_xlabel('Movie Rank (log)')
ax.set_ylabel('Number of Ratings (log)')
ax.set_title('The Long Tail: Movie Popularity')
ax.legend(fontsize=9)

# Rating concentration curve
ax2 = axes[1]
cum_pct = movie_rating_counts.cumsum() / movie_rating_counts.sum() * 100
ax2.plot(range(len(cum_pct)), cum_pct.values, lw=2, color='darkorange')
for target in [50, 80, 95]:
    idx_pos = int((cum_pct >= target).argmax())
    ax2.axhline(target, ls=':', color='grey', lw=1)
    ax2.axvline(idx_pos, ls=':', color='grey', lw=1)
    ax2.text(idx_pos + 500, target - 4, f'{target}%: {idx_pos:,} movies',
             fontsize=8, color='dimgrey')
ax2.set_xlabel('Number of Movies (ranked by popularity)')
ax2.set_ylabel('Cumulative % of All Ratings')
ax2.set_title('Rating Concentration Curve')

plt.suptitle('Movie Popularity Analysis', fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig('eda_03_movie_longtail.png', dpi=150, bbox_inches='tight')
plt.show()

print('Top 10 most-rated movies:')
print(movie_pop_named.head(10).to_string())

%%time

# Analyze the genres
movies_exploded = (movies.assign(genre=movies['genres'].str.split('|')).explode('genre').query("genre != '(no genres listed)'"))
genre_ratings = (all_ratings_eda.merge(movies_exploded[['movieId','genre']], on='movieId', how='left').dropna(subset=['genre']))

genre_stats = genre_ratings.groupby('genre').agg(n_ratings     = ('rating','count'),
                                                 mean_rating   = ('rating','mean'),
                                                 median_rating = ('rating','median'),
                                                 n_movies      = ('movieId','nunique'),).sort_values('n_ratings', ascending = False)

fig, axes = plt.subplots(1, 2, figsize=(20, 6))

# Volume per genre
ax = axes[0]
palette = sns.color_palette('husl', len(genre_stats))
ax.barh(genre_stats.index, genre_stats['n_ratings'], color=palette[::-1])
ax.invert_yaxis()
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'{x/1e6:.0f}M'))
ax.set_xlabel('Total Ratings')
ax.set_title('Rating Volume by Genre')

# Mean rating per genre (dot plot)
ax2 = axes[1]
gs_sorted = genre_stats.sort_values('mean_rating')
g_mean = all_ratings_eda['rating'].mean()
colors = ['#2ecc71' if v >= g_mean + 0.15
          else '#e74c3c' if v < g_mean - 0.15
          else '#3498db' for v in gs_sorted['mean_rating']]
ax2.scatter(gs_sorted['mean_rating'], range(len(gs_sorted)),
            s=gs_sorted['n_movies']/15, c=colors, alpha=0.85, zorder=3)
ax2.set_yticks(range(len(gs_sorted)))
ax2.set_yticklabels(gs_sorted.index)
ax2.axvline(g_mean, color='black', ls='--', lw=1.2,
            label=f'Global mean ({g_mean:.2f})')
ax2.set_xlabel('Mean Star Rating')
ax2.set_title('Mean Rating by Genre\n(bubble size = # movies in genre)')
ax2.legend(fontsize=9)

plt.suptitle('Genre Analysis', fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig('eda_04_genre_analysis.png', dpi=150, bbox_inches='tight')
plt.show()
print(genre_stats[['n_movies','n_ratings','mean_rating','median_rating']].to_string())

# Genre stats
def movie_level_agg(df):

    return df.groupby('movieId').agg(n = ('rating','count'), s=('rating','sum'))

ml_tr = movie_level_agg(train_ratings)
ml_te = movie_level_agg(test_ratings)
ml    = ml_tr.add(ml_te, fill_value=0).fillna(0)
ml['mean_rating'] = ml['s'] / ml['n'].clip(lower=1)
del ml_tr, ml_te; gc.collect()

# Explode the movie catalogue by genre
movies_exp = (movies.assign(genre = movies['genres'].str.split('|')).explode('genre').query("genre != '(no genres listed)'")[['movieId','genre']])
genre_movie = movies_exp.merge(ml, on='movieId', how='left').fillna(0)
del ml; gc.collect()

# Genre stats
genre_stats = genre_movie.groupby('genre').agg(n_ratings   = ('n', 'sum'),
                                               sum_rating  = ('s', 'sum'),
                                               n_movies    = ('movieId','nunique'),)

genre_stats['mean_rating']   = genre_stats['sum_rating'] / genre_stats['n_ratings'].clip(lower=1)
genre_stats['median_rating'] = (genre_movie.groupby('genre')['mean_rating'].median())

genre_stats = genre_stats.drop(columns='sum_rating').sort_values('n_ratings', ascending=False)
del genre_movie; gc.collect()

# Generate co-occurence rates for the entire catalogue
genre_list_per_movie = (
        movies.assign(genre=movies['genres'].str.split('|'))
        .set_index('movieId')['genre']
    )
top_genres = genre_stats.index[:15].tolist()
g2i = {g: i for i, g in enumerate(top_genres)}

cooc = np.zeros((len(top_genres), len(top_genres)), dtype=np.int32)

for glist in genre_list_per_movie:

    gs = [g for g in glist if g in g2i]

    for g in gs:

        cooc[g2i[g], g2i[g]] += 1

    for a, b in combinations(gs, 2):

        cooc[g2i[a], g2i[b]] += 1
        cooc[g2i[b], g2i[a]] += 1

cooc_df = pd.DataFrame(cooc, index=top_genres, columns=top_genres)
diag = np.diag(cooc).astype(float)
jaccard = cooc / (diag[:, None] + diag[None, :] - cooc + 1e-9)
jaccard_df = pd.DataFrame(jaccard, index=top_genres, columns=top_genres)

fig, axes = plt.subplots(1, 2, figsize = (20, 6))

sns.heatmap(cooc_df/1000, ax=axes[0], cmap='YlOrRd', annot = False, linewidths = 0.3, cbar_kws = {'label': 'Co-occurrences (thousands)'})

axes[0].set_title('Genre Co-occurrence Count')
mask = np.eye(len(top_genres), dtype=bool)

sns.heatmap(jaccard_df, ax=axes[1], cmap='Blues', mask=mask,
                annot=True, fmt='.2f', annot_kws={'size':7},
                linewidths=0.3, vmin=0, vmax=0.5)
axes[1].set_title('Genre Co-occurrence (Jaccard) 1.0 = always appear together')

plt.suptitle('Genre Co-occurrence in the Catalogue', fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig('eda_05_genre_cooccurrence.png', dpi=150, bbox_inches='tight')
plt.show()

# Temporal trends
SECS_PER_YEAR = 365.2425 * 86400

def year_agg(df):

    year = ((df['timestamp'].astype('float32') / SECS_PER_YEAR) + 1970).astype('int16')
    tmp  = pd.DataFrame({'year': year, 'r': df['rating'], 'u': df['userId']})
    out  = tmp.groupby('year').agg(n=('r','count'), s=('r','sum'), nu=('u','nunique'))
    del tmp, year
    return out

yr_tr = year_agg(train_ratings)
yr_te = year_agg(test_ratings)
yearly = pd.DataFrame({
    'n_ratings'  : yr_tr['n'].add(yr_te['n'],  fill_value=0).astype(int),
    'sum_rating' : yr_tr['s'].add(yr_te['s'],  fill_value=0),
    'n_users'    : yr_tr['nu'].add(yr_te['nu'], fill_value=0).astype(int),  # approx
}).reset_index()
yearly['mean_rating'] = yearly['sum_rating'] / yearly['n_ratings'].clip(lower=1)
del yr_tr, yr_te; gc.collect()

# Monthly (only recent years, parse subset to reduce datetime overhead)
RECENT_TS = int((2013 - 1970) * SECS_PER_YEAR)
def monthly_recent(df, cutoff_ts):
    sub = df[df['timestamp'] >= cutoff_ts][['timestamp','rating']].copy()
    sub['period'] = pd.to_datetime(sub['timestamp'], unit='s').dt.to_period('M')
    out = sub.groupby('period')['rating'].count()
    del sub; return out

mo_tr = monthly_recent(train_ratings, RECENT_TS)
mo_te = monthly_recent(test_ratings,  RECENT_TS)
monthly = mo_tr.add(mo_te, fill_value=0).reset_index()
monthly.columns = ['period','n_ratings']
monthly['dt'] = monthly['period'].dt.to_timestamp()
del mo_tr, mo_te; gc.collect()

fig = plt.figure(figsize = (24, 9))
gs2 = GridSpec(2, 2, figure=fig, hspace=0.4, wspace=0.32)

ax1 = fig.add_subplot(gs2[0, :])
bars = ax1.bar(yearly['year'], yearly['n_ratings']/1e6,
               color=sns.color_palette('husl', len(yearly)))
ax1.set_xlabel('Year'); ax1.set_ylabel('Ratings (millions)')
ax1.set_title('Annual Rating Volume (1995-2023)')
ax1.bar_label(bars, fmt='{:.1f}M', padding=2, fontsize=7)

ax2 = fig.add_subplot(gs2[1, 0])
ax2.plot(monthly['dt'], monthly['n_ratings']/1e3, lw=1.5, color='steelblue')
ax2.fill_between(monthly['dt'], monthly['n_ratings']/1e3, alpha=0.15)
ax2.set_xlabel('Month'); ax2.set_ylabel('Ratings (thousands)')
ax2.set_title('Monthly Volume (2013-2023)')

ax3 = fig.add_subplot(gs2[1, 1])
ax3.plot(yearly['year'], yearly['mean_rating'], marker='o', lw=2, color='darkorange')
ax3.axhline(global_mean_all, ls='--', color='grey', lw=1.2,
            label=f'Overall mean ({global_mean_all:.2f})')
ax3.set_ylim(3.2, 4.1); ax3.set_xlabel('Year'); ax3.set_ylabel('Mean Rating')
ax3.set_title('Mean Rating Trend Over Time'); ax3.legend(fontsize=9)

plt.suptitle('Temporal Trends in Rating Activity', fontsize=14, y=1.01)
plt.savefig('eda_06_temporal_trends.png', dpi=150, bbox_inches='tight')
plt.show()

# Perform cold start analysis
lk_tr = train_ratings[train_ratings['rating'] >= RATING_THRESHOLD].groupby('userId')['rating'].count()
lk_te = test_ratings[ test_ratings['rating']  >= RATING_THRESHOLD].groupby('userId')['rating'].count()
user_liked_counts_all = lk_tr.add(lk_te, fill_value=0).astype(int).rename('liked_count')
del lk_tr, lk_te; gc.collect()

COLD_START_THRESHOLD = 5

fig, axes = plt.subplots(1, 2, figsize=(24, 5))

ax = axes[0]
ax.hist(user_liked_counts_all.clip(upper=300), bins=60,
        color='mediumseagreen', edgecolor='white', lw=0.3)
for thresh in [5, 10, 20]:
    ax.axvline(thresh, ls='--', lw=1.5,
               label=f'thresh={thresh}  ({(user_liked_counts_all<thresh).mean():.1%} cold)')
ax.set_xlabel('Liked Ratings per User (clipped at 300)')
ax.set_ylabel('Users'); ax.set_title('Liked-Rating Count Distribution')
ax.legend(fontsize=8)

ax2 = axes[1]
thr_range = np.arange(1, 101)
pct_warm  = [(user_liked_counts_all >= t).mean()*100 for t in thr_range]
ax2.plot(thr_range, pct_warm, lw=2.5, color='mediumseagreen')
ax2.axvline(COLD_START_THRESHOLD, color='crimson', ls='--', lw=2,
            label=f'Chosen = {COLD_START_THRESHOLD}')
chosen_warm = pct_warm[COLD_START_THRESHOLD-1]
ax2.axhline(chosen_warm, color='crimson', ls=':', lw=1)
ax2.annotate(f'{chosen_warm:.1f}% warm',
             xy=(COLD_START_THRESHOLD, chosen_warm),
             xytext=(COLD_START_THRESHOLD+6, chosen_warm-8),
             arrowprops=dict(arrowstyle='->', color='crimson'),
             color='crimson', fontsize=10)
ax2.set_xlabel('Minimum Liked Ratings')
ax2.set_ylabel('% Users Classified as Warm')
ax2.set_title('Warm User Rate vs Cold-Start Threshold')
ax2.legend(fontsize=9)

plt.suptitle('Cold-Start Analysis', fontsize = 14, y = 1.02)
plt.tight_layout()
plt.savefig('eda_08_cold_start.png', dpi = 150, bbox_inches = 'tight')
plt.show()
del user_liked_counts_all; gc.collect()

# Compile user stats
liked_mask = train_ratings['rating'] >= RATING_THRESHOLD
user_liked = (train_ratings[liked_mask].groupby('userId')['movieId'].apply(list))
user_stats = train_ratings.groupby('userId').agg(
    rating_count = ('rating','count'),
    mean_rating  = ('rating','mean'),
    rating_std   = ('rating','std'),
    last_ts      = ('timestamp','max'),
    first_ts     = ('timestamp','min'),)

user_stats['active_days'] = ((user_stats['last_ts'] - user_stats['first_ts']) / 86_400).clip(lower=1).astype('float32')
user_stats['rating_freq'] = (user_stats['rating_count'] / user_stats['active_days']).astype('float32')
user_stats['mean_rating'] = user_stats['mean_rating'].astype('float32')

# Users with a single rating have undefined std; fill with 0
user_stats['rating_std']  = user_stats['rating_std'].fillna(0).astype('float32')
print('user_liked:', len(user_liked), 'users')
print('user_stats:', user_stats.shape)

# Extract genres and perform movie one hot encoding
genres_exploded = (movies.assign(genre = movies['genres'].str.split('|')).explode('genre').query("genre != '(no genres listed)'"))

ALL_GENRES = sorted(genres_exploded['genre'].unique())
print(f'Genres ({len(ALL_GENRES)}):', ALL_GENRES)

movie_genre = pd.get_dummies(genres_exploded.set_index('movieId')['genre'], prefix='g')
movie_genre = movie_genre.groupby(level=0).max()
genre_cols  = movie_genre.columns.tolist()

print('movie_genre shape:', movie_genre.shape)

# Generate Bayesian scores and movie statistics
GLOBAL_MEAN = train_ratings['rating'].mean()
C = 50

movie_stats = (
    train_ratings
    .groupby('movieId')
    .agg(global_count    = ('rating','count'),
         global_mean     = ('rating','mean'),
         rating_variance = ('rating','var'))
)
movie_stats['bayesian_score'] = (
    (movie_stats['global_count'] * movie_stats['global_mean'] + C * GLOBAL_MEAN)
    / (movie_stats['global_count'] + C)
)
# log1p compresses the very wide count range so the model sees a
# smooth popularity signal rather than a skewed raw count
movie_stats['log_count'] = np.log1p(movie_stats['global_count']).astype('float32')
# Single-rating movies have undefined variance; treat as 0
movie_stats['rating_variance'] = movie_stats['rating_variance'].fillna(0).astype('float32')

# Tag count from tags.csv (loaded in Section 1) — log1p scaled.
# Tags are user-applied descriptors; a high tag count signals that a
# movie has attracted community engagement beyond a simple star rating.
movie_tag_count = tags.groupby('movieId')['tag'].count().rename('tag_count')
movie_stats = movie_stats.join(movie_tag_count, how='left')
movie_stats['tag_count'] = np.log1p(movie_stats['tag_count'].fillna(0)).astype('float32')
del movie_tag_count

movie_feat = movie_stats.join(movie_genre, how='left').fillna(0)
print(f'movie_feat shape: {movie_feat.shape}   global mean: {GLOBAL_MEAN:.4f}')
movie_feat[['global_count','log_count','global_mean','bayesian_score','rating_variance','tag_count']].describe()

# User genre affinity matrix
liked_fe = (train_ratings[train_ratings['rating'] >= RATING_THRESHOLD][['userId','movieId','rating']])
uid_uniq = liked_fe['userId'].unique()
mid_uniq = liked_fe['movieId'].unique()

u2i = {u: i for i, u in enumerate(uid_uniq)}
m2i = {m: i for i, m in enumerate(mid_uniq)}

rows_fe = liked_fe['userId'].map(u2i).values
cols_fe = liked_fe['movieId'].map(m2i).values

del liked_fe; gc.collect()

# Binary matrix: 1 where user liked a movie, 0 elsewhere.
# Using rating values in the numerator (previous approach) inflates genres
# where the user gave 5-star ratings vs 4-star ratings for the same number
# of films. Count-based gives the clean interpretation: fraction of a
# user's liked movies that belong to each genre.
R_bin = csr_matrix(
    (np.ones(len(rows_fe), dtype='float32'), (rows_fe, cols_fe)),
    shape=(len(uid_uniq), len(mid_uniq)),
)
del rows_fe, cols_fe; gc.collect()

# Genre matrix aligned to the liked-movie order (float32)
G = movie_genre.reindex(mid_uniq).fillna(0).values.astype('float32')

# Sparse-binary x dense multiply -> (n_users x n_genres) count matrix
uga_counts   = R_bin @ G
liked_counts = np.asarray(R_bin.sum(axis=1)).flatten()
del R_bin, G; gc.collect()

user_genre_affinity = pd.DataFrame(
    uga_counts / liked_counts[:, None].clip(min=1),
    index=uid_uniq,
    columns=genre_cols,
    dtype='float32',
)
del uga_counts, liked_counts; gc.collect()

print(f'user_genre_affinity: {user_genre_affinity.shape}  |  '
      f'{user_genre_affinity.memory_usage(deep=True).sum()/1e6:.0f} MB')