Capstone - Complete Analytics System

# Python: Extracting Possession Sequences
import pandas as pd
import numpy as np

def extract_sequences(events):
    """Extract and characterize possession sequences."""
    events = events.sort_values("index").copy()

    # Detect possession changes
    events["team_change"] = events["team"] != events["team"].shift(1)
    events["new_sequence"] = (
        events["team_change"] |
        events["type"].isin(["Half Start", "Starting XI", "Kick Off"])
    )
    events["sequence_id"] = events["new_sequence"].cumsum()

    # Summarize each sequence
    def summarize_sequence(group):
        return pd.Series({
            "team": group["team"].iloc[0],
            "start_minute": group["minute"].min(),
            "start_second": group["second"].min(),
            "end_minute": group["minute"].max(),
            "end_second": group["second"].max(),
            "n_events": len(group),
            "n_passes": (group["type"] == "Pass").sum(),
            "n_carries": (group["type"] == "Carry").sum(),
            "n_dribbles": (group["type"] == "Dribble").sum(),
            "ends_with_shot": (group["type"] == "Shot").any(),
            "ends_with_goal": ((group["type"] == "Shot") &
                              (group["shot_outcome"] == "Goal")).any(),
            "total_xg": group.loc[group["type"] == "Shot", "shot_statsbomb_xg"].sum(),
            "start_x": group["location_x"].dropna().iloc[0] if len(group["location_x"].dropna()) > 0 else np.nan,
            "end_x": group["location_x"].dropna().iloc[-1] if len(group["location_x"].dropna()) > 0 else np.nan,
            "max_x": group["location_x"].max(),
        })

    sequences = events.groupby("sequence_id").apply(summarize_sequence).reset_index()

    # Derived metrics
    sequences["duration"] = (
        (sequences["end_minute"] * 60 + sequences["end_second"]) -
        (sequences["start_minute"] * 60 + sequences["start_second"])
    )
    sequences["territory_gained"] = sequences["end_x"] - sequences["start_x"]
    sequences["reached_final_third"] = sequences["max_x"] >= 80
    sequences["reached_box"] = sequences["max_x"] >= 102

    # Sequence type classification
    def classify_sequence(n_passes):
        if n_passes <= 2:
            return "Counter/Direct"
        elif n_passes <= 5:
            return "Quick Attack"
        elif n_passes <= 10:
            return "Build-up"
        else:
            return "Extended Build-up"

    sequences["sequence_type"] = sequences["n_passes"].apply(classify_sequence)

    return sequences

sequences = extract_sequences(events)

print("Sequence Type Distribution:")
print(sequences["sequence_type"].value_counts())

print("\nAverage sequence by type:")
print(sequences.groupby("sequence_type").agg({
    "sequence_id": "count",
    "n_passes": "mean",
    "ends_with_shot": "mean",
    "total_xg": "mean"
}).rename(columns={"sequence_id": "count"}))

# Python: Multi-dimensional Sequence Classification
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans

def classify_sequences_advanced(events):
    """Classify sequences using multiple dimensions."""
    events = events.sort_values(["sequence_id", "index"]).copy()

    def compute_features(group):
        locs_x = group["location_x"].dropna()
        locs_y = group["location_y"].dropna()
        times = group["minute"] * 60 + group["second"]

        duration = times.max() - times.min() if len(times) > 1 else 0
        n_events = len(group)

        # Spatial features
        start_x = locs_x.iloc[0] if len(locs_x) > 0 else 60
        end_x = locs_x.iloc[-1] if len(locs_x) > 0 else 60
        max_x = locs_x.max() if len(locs_x) > 0 else 60
        min_x = locs_x.min() if len(locs_x) > 0 else 60

        # Width
        min_y = locs_y.min() if len(locs_y) > 0 else 40
        max_y = locs_y.max() if len(locs_y) > 0 else 40
        width_used = max_y - min_y

        # Directness
        net_progress = end_x - start_x
        x_changes = locs_x.diff().abs().sum() if len(locs_x) > 1 else 0
        directness = net_progress / x_changes if x_changes > 0 else 0

        # Pass analysis
        passes = group[group["type"] == "Pass"]
        n_passes = len(passes)

        return pd.Series({
            "team": group["team"].iloc[0],
            "duration": duration,
            "n_events": n_events,
            "events_per_second": n_events / duration if duration > 0 else n_events,
            "start_x": start_x,
            "end_x": end_x,
            "max_x": max_x,
            "net_progress": net_progress,
            "directness": directness,
            "width_used": width_used,
            "n_passes": n_passes,
            "ends_with_shot": (group["type"] == "Shot").any(),
            "total_xg": group.loc[group["type"] == "Shot", "shot_statsbomb_xg"].sum(),
        })

    seq_features = events.groupby("sequence_id").apply(compute_features).reset_index()

    # Multi-dimensional classification
    def tempo_class(d):
        if d <= 4: return "Lightning"
        if d <= 10: return "Fast"
        if d <= 20: return "Medium"
        if d <= 40: return "Patient"
        return "Extended"

    def width_class(w):
        if w < 30: return "Narrow"
        if w < 45: return "Medium"
        if w < 55: return "Wide"
        return "Full"

    def directness_class(d):
        if d > 0.7: return "Very Direct"
        if d > 0.4: return "Direct"
        if d > 0.1: return "Balanced"
        if d > -0.1: return "Recycling"
        return "Retreating"

    seq_features["tempo_class"] = seq_features["duration"].apply(tempo_class)
    seq_features["width_class"] = seq_features["width_used"].apply(width_class)
    seq_features["directness_class"] = seq_features["directness"].apply(directness_class)

    return seq_features

# Apply classification
sequences_classified = classify_sequences_advanced(events)

# Analyze effectiveness
effectiveness = sequences_classified.groupby(["tempo_class", "directness_class"]).agg({
    "sequence_id": "count",
    "ends_with_shot": "mean",
    "total_xg": "mean",
    "max_x": "mean"
}).rename(columns={"sequence_id": "count"}).reset_index()

effectiveness = effectiveness[effectiveness["count"] >= 10].sort_values("total_xg", ascending=False)
print("Effectiveness by Tempo and Directness:")
print(effectiveness)

# Python: K-means Clustering of Sequences
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

def cluster_sequences(sequence_features, n_clusters=6):
    """Cluster sequences by style characteristics."""
    # Select features
    cluster_vars = ["duration", "n_passes", "directness",
                    "width_used", "net_progress"]

    # Prepare data
    cluster_data = sequence_features[["sequence_id"] + cluster_vars].dropna()

    # Standardize
    scaler = StandardScaler()
    scaled_data = scaler.fit_transform(cluster_data[cluster_vars])

    # Find optimal k
    silhouette_scores = []
    for k in range(2, 11):
        km = KMeans(n_clusters=k, random_state=42, n_init=25)
        labels = km.fit_predict(scaled_data)
        score = silhouette_score(scaled_data, labels)
        silhouette_scores.append(score)

    optimal_k = np.argmax(silhouette_scores) + 2
    print(f"Optimal clusters: {optimal_k}")

    # Final clustering
    final_km = KMeans(n_clusters=n_clusters, random_state=42, n_init=25)
    cluster_data["cluster"] = final_km.fit_predict(scaled_data)

    # Analyze clusters
    profiles = cluster_data.groupby("cluster").agg({
        "sequence_id": "count",
        "duration": "mean",
        "n_passes": "mean",
        "directness": "mean",
        "width_used": "mean",
        "net_progress": "mean"
    }).rename(columns={"sequence_id": "count"}).reset_index()

    # Name clusters
    def name_cluster(row):
        if row["n_passes"] < 3 and row["directness"] > 0.5:
            return "Counter Attack"
        if row["n_passes"] > 10 and row["width_used"] > 45:
            return "Tiki-Taka"
        if row["n_passes"] > 8 and row["directness"] < 0.2:
            return "Patient Buildup"
        if row["width_used"] > 50 and row["directness"] > 0.3:
            return "Wing Play"
        if row["directness"] > 0.6:
            return "Direct Attack"
        return "Mixed"

    profiles["style_name"] = profiles.apply(name_cluster, axis=1)

    return cluster_data, profiles

cluster_data, profiles = cluster_sequences(sequences_classified)
print(profiles)

# Python: Counter-Attack Detection
import pandas as pd
import numpy as np

def detect_counter_attacks(events):
    """Detect and classify counter-attacking sequences."""
    events = events.sort_values("index").copy()

    # Identify possession changes
    events["team_prev"] = events["team"].shift(1)
    events["possession_change"] = events["team"] != events["team_prev"]
    events["sequence_id"] = events["possession_change"].cumsum()

    # Get trigger info (previous event before possession won)
    events["prev_type"] = events["type"].shift(1)
    events["prev_outcome"] = events.get("pass_outcome", pd.Series()).shift(1)
    events["prev_x"] = events["location_x"].shift(1)

    # Calculate sequence characteristics
    def analyze_sequence(group):
        locs_x = group["location_x"].dropna()
        times = group["minute"] * 60 + group["second"]
        passes = group[group["type"] == "Pass"]

        start_x = locs_x.iloc[0] if len(locs_x) > 0 else 60
        end_x = locs_x.iloc[-1] if len(locs_x) > 0 else 60
        max_x = locs_x.max() if len(locs_x) > 0 else 60
        duration = times.max() - times.min() if len(times) > 1 else 0

        # Count forward passes
        n_forward = 0
        for i in range(len(group) - 1):
            if group.iloc[i]["type"] == "Pass":
                curr_x = group.iloc[i]["location_x"]
                next_x = group.iloc[i + 1]["location_x"]
                if pd.notna(curr_x) and pd.notna(next_x) and next_x - curr_x > 5:
                    n_forward += 1

        return pd.Series({
            "team": group["team"].iloc[0],
            "start_x": start_x,
            "end_x": end_x,
            "max_x": max_x,
            "duration": duration,
            "x_progress": end_x - start_x,
            "n_passes": len(passes),
            "n_forward_passes": n_forward,
            "ends_shot": (group["type"] == "Shot").any(),
            "total_xg": group.loc[group["type"] == "Shot", "shot_statsbomb_xg"].sum(),
            "prev_type": group["prev_type"].iloc[0],
            "prev_outcome": group["prev_outcome"].iloc[0] if "prev_outcome" in group else None,
        })

    seq_chars = events.groupby("sequence_id").apply(analyze_sequence).reset_index()
    seq_chars["pace"] = np.where(
        seq_chars["duration"] > 0,
        seq_chars["x_progress"] / seq_chars["duration"],
        seq_chars["x_progress"]
    )

    # Counter-attack criteria
    counters = seq_chars[
        (seq_chars["start_x"] < 60) &
        (seq_chars["duration"] <= 12) &
        (seq_chars["x_progress"] >= 35) &
        (seq_chars["n_passes"] <= 5) &
        ((seq_chars["n_forward_passes"] >= seq_chars["n_passes"] * 0.6) |
         (seq_chars["n_passes"] <= 2))
    ].copy()

    # Classify counter type
    def classify_counter(row):
        if row["prev_type"] == "Interception":
            return "Interception Counter"
        if row["prev_type"] == "Tackle":
            return "Tackle Counter"
        if row["prev_outcome"] == "Incomplete":
            return "Failed Pass Counter"
        if row["prev_type"] == "Clearance":
            return "Transition Counter"
        if row["prev_type"] in ["Goal Keeper", "Save"]:
            return "GK Distribution"
        return "Other"

    counters["counter_type"] = counters.apply(classify_counter, axis=1)

    return counters

counters = detect_counter_attacks(events)

print("Counter-attack Statistics:")
print(f"Total counters detected: {len(counters)}")
print(f"Shot rate: {counters[\"ends_shot\"].mean():.3f}")
print(f"Average xG: {counters[\"total_xg\"].mean():.4f}")
print()

print("By Trigger Type:")
print(counters.groupby("counter_type").agg({
    "sequence_id": "count",
    "ends_shot": "mean",
    "total_xg": "mean",
    "pace": "mean"
}).rename(columns={"sequence_id": "count"}))

# Python: Counter-Pressing and PPDA Analysis
import pandas as pd
import numpy as np

def analyze_pressing(events):
    """Analyze counter-pressing and calculate PPDA."""
    events = events.sort_values("index").copy()

    # Identify possession changes
    events["next_team"] = events["team"].shift(-1)
    events["prev_team"] = events["team"].shift(1)

    events["lost_possession"] = (
        (events["team"] != events["next_team"]) &
        events["type"].isin(["Pass", "Dribble", "Carry"])
    )
    events["won_possession"] = events["team"] != events["prev_team"]

    # Track losses and recoveries
    losses = events[events["lost_possession"]][
        ["index", "team", "location_x", "location_y", "minute", "second"]
    ].copy()
    losses.columns = ["loss_index", "loss_team", "loss_x", "loss_y", "loss_min", "loss_sec"]
    losses["loss_time"] = losses["loss_min"] * 60 + losses["loss_sec"]

    recoveries = events[
        events["won_possession"] |
        events["type"].isin(["Interception", "Tackle", "Ball Recovery"])
    ][["index", "team", "location_x", "location_y", "minute", "second", "type"]].copy()
    recoveries.columns = ["rec_index", "rec_team", "rec_x", "rec_y", "rec_min", "rec_sec", "rec_type"]
    recoveries["rec_time"] = recoveries["rec_min"] * 60 + recoveries["rec_sec"]

    # Match losses to recoveries
    counter_press_data = []
    for _, loss in losses.iterrows():
        team_recoveries = recoveries[
            (recoveries["rec_team"] == loss["loss_team"]) &
            (recoveries["rec_index"] > loss["loss_index"])
        ]

        if len(team_recoveries) > 0:
            recovery = team_recoveries.iloc[0]
            time_diff = recovery["rec_time"] - loss["loss_time"]

            loss_y = loss["loss_y"] if pd.notna(loss["loss_y"]) else 40
            rec_y = recovery["rec_y"] if pd.notna(recovery["rec_y"]) else 40
            distance = np.sqrt(
                (recovery["rec_x"] - loss["loss_x"])**2 +
                (rec_y - loss_y)**2
            )

            is_counter_press = (time_diff <= 5) and (distance <= 25)

            counter_press_data.append({
                "loss_team": loss["loss_team"],
                "loss_x": loss["loss_x"],
                "recovery_time": time_diff,
                "recovery_distance": distance,
                "is_counter_press": is_counter_press
            })

    cp_df = pd.DataFrame(counter_press_data)

    # PPDA calculation
    defensive_actions = ["Tackle", "Interception", "Foul Committed", "Duel"]

    teams = events["team"].unique()
    ppda_data = []

    for team in teams:
        # Opponent passes in their attacking half (our defensive half)
        opponent_passes = len(events[
            (events["team"] != team) &
            (events["type"] == "Pass") &
            (events["location_x"] > 60)  # Opponent attacking half
        ])

        # Our defensive actions in opponent half
        our_def_actions = len(events[
            (events["team"] == team) &
            (events["type"].isin(defensive_actions)) &
            (events["location_x"] > 60)
        ])

        ppda = opponent_passes / max(our_def_actions, 1)
        ppda_data.append({"team": team, "ppda": ppda})

    return cp_df, pd.DataFrame(ppda_data)

cp_df, ppda_df = analyze_pressing(events)

print("Counter-Press Analysis:")
print(f"Total losses: {len(cp_df)}")
print(f"Counter-presses: {cp_df[\"is_counter_press\"].sum()}")
print(f"Success rate: {cp_df[\"is_counter_press\"].mean():.3f}")
successful_cp = cp_df[cp_df["is_counter_press"]]
if len(successful_cp) > 0:
    print(f"Avg recovery time: {successful_cp[\"recovery_time\"].mean():.2f}s")

print("\nPPDA by Team:")
print(ppda_df)

# Python: Buildup Pattern Analysis
import pandas as pd
import numpy as np

def analyze_buildup_patterns(events):
    """Analyze team buildup patterns in attacking sequences."""
    # Focus on sequences reaching final third
    sequence_max_x = events.groupby("sequence_id")["location_x"].transform("max")
    attacking_events = events[sequence_max_x >= 80].copy()

    # Categorize passes
    passes = attacking_events[attacking_events["type"] == "Pass"].copy()

    def safe_subtract(end, start):
        if pd.isna(end) or pd.isna(start):
            return 0
        return end - start

    passes["x_progress"] = passes.apply(
        lambda r: safe_subtract(
            r["pass_end_location"][0] if isinstance(r["pass_end_location"], list) else None,
            r["location_x"]
        ), axis=1
    )

    # Pass direction
    passes["pass_direction"] = pd.cut(
        passes["x_progress"],
        bins=[-np.inf, -10, 10, np.inf],
        labels=["Backward", "Lateral", "Forward"]
    )

    # Pass zone
    passes["pass_zone"] = pd.cut(
        passes["location_x"],
        bins=[0, 40, 80, 120],
        labels=["Defensive Third", "Middle Third", "Final Third"]
    )

    # Pass type
    def classify_pass_type(row):
        if row.get("pass_cross"):
            return "Cross"
        if row.get("pass_through_ball"):
            return "Through Ball"
        y = row.get("location_y", 40)
        if y < 20 or y > 60:
            return "Wide"
        return "Central"

    passes["pass_type"] = passes.apply(classify_pass_type, axis=1)

    # Sequence-level buildup profile
    def sequence_profile(group):
        return pd.Series({
            "team": group["team"].iloc[0],
            "total_passes": len(group),
            "pct_forward": (group["pass_direction"] == "Forward").mean(),
            "pct_backward": (group["pass_direction"] == "Backward").mean(),
            "pct_lateral": (group["pass_direction"] == "Lateral").mean(),
            "pct_wide": (group["pass_type"] == "Wide").mean(),
            "pct_central": (group["pass_type"] == "Central").mean(),
            "through_balls": (group["pass_type"] == "Through Ball").sum(),
            "crosses": (group["pass_type"] == "Cross").sum(),
        })

    buildup_profile = passes.groupby("sequence_id").apply(sequence_profile).reset_index()

    # Team-level summary
    team_buildup = buildup_profile.groupby("team").agg({
        "sequence_id": "count",
        "total_passes": "mean",
        "pct_forward": "mean",
        "pct_wide": "mean",
        "through_balls": "mean",
        "crosses": "mean"
    }).rename(columns={"sequence_id": "sequences"}).reset_index()

    return {"sequences": buildup_profile, "teams": team_buildup}

buildup = analyze_buildup_patterns(events)
print(buildup["teams"])

# Python: Visualize Common Buildup Paths
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.sankey import Sankey

def visualize_buildup_flow(events, team_name):
    """Visualize passing flow between zones."""
    # Get team passes
    team_passes = events[
        (events["team"] == team_name) &
        (events["type"] == "Pass")
    ].copy()

    # Define zone assignment function
    def assign_zone(x, y):
        if x < 40:
            prefix = "Def"
        elif x < 80:
            prefix = "Mid"
        else:
            prefix = "Att"

        if y < 26.7:
            suffix = "Right"
        elif y < 53.3:
            suffix = "Center"
        else:
            suffix = "Left"

        return f"{prefix} {suffix}"

    team_passes["zone_from"] = team_passes.apply(
        lambda r: assign_zone(r["location_x"], r["location_y"]), axis=1
    )

    # Get end location
    def get_end_coords(row):
        loc = row.get("pass_end_location")
        if isinstance(loc, list) and len(loc) >= 2:
            return loc[0], loc[1]
        return None, None

    team_passes["end_x"], team_passes["end_y"] = zip(
        *team_passes.apply(get_end_coords, axis=1)
    )
    team_passes = team_passes.dropna(subset=["end_x", "end_y"])

    team_passes["zone_to"] = team_passes.apply(
        lambda r: assign_zone(r["end_x"], r["end_y"]), axis=1
    )

    # Count zone transitions
    zone_flow = team_passes.groupby(["zone_from", "zone_to"]).size().reset_index(name="passes")
    zone_flow = zone_flow[zone_flow["passes"] >= 5]

    # Create chord-style visualization
    fig, ax = plt.subplots(figsize=(14, 10))

    # Define zone positions for visualization
    zone_positions = {
        "Def Left": (20, 67), "Def Center": (20, 40), "Def Right": (20, 13),
        "Mid Left": (60, 67), "Mid Center": (60, 40), "Mid Right": (60, 13),
        "Att Left": (100, 67), "Att Center": (100, 40), "Att Right": (100, 13),
    }

    # Draw pitch background
    ax.set_facecolor("#1a472a")
    ax.set_xlim(0, 120)
    ax.set_ylim(0, 80)

    # Draw zone points
    for zone, (x, y) in zone_positions.items():
        ax.scatter(x, y, s=500, c="white", zorder=10)
        ax.annotate(zone.split()[1], (x, y), ha="center", va="center",
                    fontsize=8, fontweight="bold", zorder=11)

    # Draw flows
    max_passes = zone_flow["passes"].max()
    for _, row in zone_flow.iterrows():
        if row["zone_from"] in zone_positions and row["zone_to"] in zone_positions:
            x1, y1 = zone_positions[row["zone_from"]]
            x2, y2 = zone_positions[row["zone_to"]]
            width = row["passes"] / max_passes * 5

            ax.annotate("", xy=(x2, y2), xytext=(x1, y1),
                        arrowprops=dict(
                            arrowstyle="-|>",
                            connectionstyle="arc3,rad=0.1",
                            lw=width,
                            color="yellow",
                            alpha=0.6
                        ))

    ax.set_aspect("equal")
    ax.set_title(f"{team_name} - Passing Flow Between Zones", fontsize=14)
    ax.axis("off")
    plt.show()

visualize_buildup_flow(events, "Barcelona")

# Python: Markov Chain Possession Model
import pandas as pd
import numpy as np
from collections import defaultdict

def build_possession_markov(events, n_x_zones=4, n_y_zones=3):
    """Build Markov chain model for possession flow."""
    zone_width = 120 / n_x_zones
    zone_height = 80 / n_y_zones

    # Assign zones
    events = events.dropna(subset=["location_x", "location_y"]).copy()
    events["zone"] = events.apply(
        lambda r: f"Z{int(np.ceil(r['location_x']/zone_width))}_{int(np.ceil(r['location_y']/zone_height))}",
        axis=1
    )

    # Define states
    def get_state(row, next_row):
        if row["type"] == "Shot":
            if row.get("shot_outcome") == "Goal":
                return "GOAL"
            return "SHOT_NO_GOAL"
        if next_row is not None and next_row.get("team") != row["team"]:
            return "TURNOVER"
        return row["zone"]

    # Build transition counts
    events = events.sort_values("index").reset_index(drop=True)
    transitions = defaultdict(lambda: defaultdict(int))

    for i in range(len(events) - 1):
        current = events.iloc[i]
        next_event = events.iloc[i + 1]

        current_state = get_state(current, next_event)
        next_state = get_state(next_event, events.iloc[i + 2] if i + 2 < len(events) else None)

        transitions[current_state][next_state] += 1

    # Build transition matrix
    states = list(set(transitions.keys()) | set(
        s for d in transitions.values() for s in d.keys()
    ))
    n_states = len(states)
    state_idx = {s: i for i, s in enumerate(states)}

    trans_matrix = np.zeros((n_states, n_states))
    for from_state, to_dict in transitions.items():
        for to_state, count in to_dict.items():
            trans_matrix[state_idx[from_state], state_idx[to_state]] = count

    # Normalize rows
    row_sums = trans_matrix.sum(axis=1, keepdims=True)
    row_sums[row_sums == 0] = 1  # Avoid division by zero
    trans_matrix = trans_matrix / row_sums

    return trans_matrix, states, state_idx

def calculate_goal_probability(trans_matrix, states, state_idx):
    """Calculate probability of reaching GOAL from each zone."""
    n = len(states)
    goal_idx = state_idx.get("GOAL")

    if goal_idx is None:
        return {}

    # Identify absorbing states
    absorbing = ["GOAL", "SHOT_NO_GOAL", "TURNOVER"]
    absorbing_idx = [state_idx[s] for s in absorbing if s in state_idx]
    transient_idx = [i for i in range(n) if i not in absorbing_idx]

    if not transient_idx:
        return {}

    # Extract Q matrix (transient to transient transitions)
    Q = trans_matrix[np.ix_(transient_idx, transient_idx)]

    # Extract R matrix (transient to absorbing transitions)
    R = trans_matrix[np.ix_(transient_idx, absorbing_idx)]

    # Fundamental matrix: N = (I - Q)^-1
    try:
        N = np.linalg.inv(np.eye(len(transient_idx)) - Q)
    except np.linalg.LinAlgError:
        return {}

    # Absorption probabilities: B = N * R
    B = N @ R

    # Get goal probability (column index in absorbing)
    goal_col = absorbing.index("GOAL") if "GOAL" in absorbing else None
    if goal_col is None:
        return {}

    goal_probs = {}
    transient_states = [states[i] for i in transient_idx]
    for i, state in enumerate(transient_states):
        goal_probs[state] = B[i, goal_col]

    return goal_probs

trans_matrix, states, state_idx = build_possession_markov(events)
goal_probs = calculate_goal_probability(trans_matrix, states, state_idx)

print("Goal probability from each zone:")
for zone, prob in sorted(goal_probs.items(), key=lambda x: -x[1])[:10]:
    print(f"  {zone}: {prob:.3f}")

# Python: Sequential Pattern Mining
import pandas as pd
import numpy as np
from collections import defaultdict
from itertools import combinations

def prepare_sequence_data(events):
    """Prepare events as action sequences."""
    events = events.sort_values(["sequence_id", "index"]).copy()

    # Create action labels
    def create_action(row):
        zone = "DefThird" if row["location_x"] < 40 else \
               "MidThird" if row["location_x"] < 80 else "AttThird"
        return f"{row['type']}_{zone}"

    events["action"] = events.apply(create_action, axis=1)

    # Build sequence strings
    sequence_data = events.groupby("sequence_id").agg(
        sequence=("action", lambda x: " -> ".join(x)),
        team=("team", "first"),
        ended_shot=("type", lambda x: (x == "Shot").any()),
        ended_goal=("shot_outcome", lambda x: (x == "Goal").any() if "Goal" in x.values else False)
    ).reset_index()

    return sequence_data

def find_frequent_patterns(sequence_data, min_support=0.05, max_length=3):
    """Find frequently occurring action patterns."""
    pattern_sequences = defaultdict(list)

    for idx, row in sequence_data.iterrows():
        actions = row["sequence"].split(" -> ")

        # Extract n-grams of different lengths
        for length in range(2, min(len(actions), max_length) + 1):
            for i in range(len(actions) - length + 1):
                pattern = " -> ".join(actions[i:i + length])
                pattern_sequences[pattern].append(idx)

    # Calculate support and other metrics
    n_sequences = len(sequence_data)
    pattern_stats = []

    for pattern, indices in pattern_sequences.items():
        unique_indices = list(set(indices))
        support = len(unique_indices) / n_sequences

        if support >= min_support:
            pattern_stats.append({
                "pattern": pattern,
                "support": support,
                "count": len(unique_indices),
                "shot_rate": sequence_data.loc[unique_indices, "ended_shot"].mean(),
                "goal_rate": sequence_data.loc[unique_indices, "ended_goal"].mean()
            })

    patterns = pd.DataFrame(pattern_stats).sort_values("support", ascending=False)
    return patterns

seq_data = prepare_sequence_data(events)
patterns = find_frequent_patterns(seq_data, min_support=0.03)

print("Most Common Patterns:")
print(patterns.head(10))

# Python: Building xGChain-style Metrics
import pandas as pd
import numpy as np

def calculate_sequence_value(events):
    """Calculate xGChain and xGBuildup for players."""
    # Get sequence xG totals
    sequence_xg = events[events["type"] == "Shot"].groupby("sequence_id").agg(
        sequence_xg=("shot_statsbomb_xg", "sum"),
        sequence_goal=("shot_outcome", lambda x: (x == "Goal").any())
    ).reset_index()

    # Join back to all events
    events_with_xg = events.merge(sequence_xg, on="sequence_id", how="left")
    events_with_xg["sequence_xg"] = events_with_xg["sequence_xg"].fillna(0)
    events_with_xg["sequence_goal"] = events_with_xg["sequence_goal"].fillna(False)

    # Calculate xGChain per player
    def player_xgchain(group):
        unique_sequences = group.drop_duplicates(subset=["sequence_id"])
        return pd.Series({
            "sequences_involved": group["sequence_id"].nunique(),
            "xgchain": unique_sequences["sequence_xg"].sum()
        })

    xgchain = events_with_xg.groupby(["player", "team"]).apply(
        player_xgchain
    ).reset_index()

    # Identify assists and shots to exclude for xGBuildup
    events_sorted = events_with_xg.sort_values("index")
    events_sorted["next_type"] = events_sorted["type"].shift(-1)

    assists = events_sorted[
        (events_sorted["type"] == "Pass") &
        (events_sorted["next_type"] == "Shot")
    ]["id"].tolist()

    shots = events_sorted[events_sorted["type"] == "Shot"]["id"].tolist()

    # Calculate xGBuildup (excluding shots and assists)
    buildup_events = events_with_xg[~events_with_xg["id"].isin(assists + shots)]

    def player_xgbuildup(group):
        unique_sequences = group.drop_duplicates(subset=["sequence_id"])
        return unique_sequences["sequence_xg"].sum()

    xgbuildup = buildup_events.groupby(["player", "team"]).apply(
        player_xgbuildup
    ).reset_index(name="xgbuildup")

    # Combine
    player_value = xgchain.merge(xgbuildup, on=["player", "team"], how="left")
    player_value["xgbuildup"] = player_value["xgbuildup"].fillna(0)

    return player_value

player_sequence_value = calculate_sequence_value(events)
print(player_sequence_value.sort_values("xgchain", ascending=False).head(10))

# Python: Sequence Similarity Using Edit Distance
import pandas as pd
import numpy as np
from Levenshtein import distance as levenshtein_distance
from scipy.cluster.hierarchy import linkage, fcluster
from scipy.spatial.distance import squareform

def calculate_sequence_similarity(events, max_sequences=500):
    """Calculate pairwise similarity between sequences."""
    events = events.sort_values(["sequence_id", "index"]).copy()

    # Convert to action codes
    def get_action_code(row):
        if row["type"] == "Pass":
            if row["location_x"] < 40:
                return "PD"
            elif row["location_x"] < 80:
                return "PM"
            return "PA"
        elif row["type"] == "Carry":
            return "C"
        elif row["type"] == "Dribble":
            return "D"
        elif row["type"] == "Shot":
            return "S"
        return "O"

    events["action_code"] = events.apply(get_action_code, axis=1)

    # Build action strings
    seq_strings = events.groupby("sequence_id").agg(
        team=("team", "first"),
        action_string=("action_code", lambda x: "-".join(x)),
        n_actions=("action_code", "count"),
        ends_shot=("type", lambda x: (x == "Shot").any()),
        total_xg=("shot_statsbomb_xg", "sum")
    ).reset_index()

    # Filter minimum length
    seq_strings = seq_strings[seq_strings["n_actions"] >= 3]

    # Sample if too large
    if len(seq_strings) > max_sequences:
        seq_strings = seq_strings.sample(n=max_sequences, random_state=42)

    # Calculate pairwise distances
    n_seq = len(seq_strings)
    distance_matrix = np.zeros((n_seq, n_seq))

    strings = seq_strings["action_string"].tolist()
    for i in range(n_seq):
        for j in range(i + 1, n_seq):
            dist = levenshtein_distance(strings[i], strings[j])
            max_len = max(len(strings[i]), len(strings[j]))
            normalized_dist = dist / max_len if max_len > 0 else 0
            distance_matrix[i, j] = normalized_dist
            distance_matrix[j, i] = normalized_dist

    # Convert to similarity
    similarity_matrix = 1 - distance_matrix

    return seq_strings, similarity_matrix

def find_similar_sequences(seq_strings, similarity_matrix, target_idx, top_n=5):
    """Find sequences most similar to target."""
    sims = similarity_matrix[target_idx].copy()
    sims[target_idx] = -1  # Exclude self

    top_indices = np.argsort(sims)[::-1][:top_n]

    similar = seq_strings.iloc[top_indices].copy()
    similar["similarity_to_target"] = sims[top_indices]

    return similar

# Run analysis
seq_strings, sim_matrix = calculate_sequence_similarity(events)

# Find sequences similar to high-xG sequence
goal_seqs = seq_strings[
    seq_strings["ends_shot"] &
    (seq_strings["total_xg"] > 0.5)
].index.tolist()

if goal_seqs:
    target_idx = seq_strings.index.get_loc(goal_seqs[0]) if goal_seqs[0] in seq_strings.index else 0
    target_row = seq_strings[seq_strings.index == goal_seqs[0]]
    print(f"Target sequence: {target_row[\"action_string\"].values[0]}")

    similar = find_similar_sequences(seq_strings.reset_index(drop=True), sim_matrix, target_idx)
    print("Similar sequences found:")
    print(similar[["action_string", "ends_shot", "similarity_to_target"]])

# Python: Hierarchical Clustering of Sequences
import pandas as pd
import numpy as np
from scipy.cluster.hierarchy import linkage, fcluster, dendrogram
from scipy.spatial.distance import squareform
import matplotlib.pyplot as plt

def cluster_similar_sequences(seq_strings, similarity_matrix, k_clusters=8):
    """Cluster sequences by similarity."""
    # Convert similarity to condensed distance matrix
    distance_matrix = 1 - similarity_matrix
    np.fill_diagonal(distance_matrix, 0)

    # Ensure symmetry and non-negative
    distance_matrix = (distance_matrix + distance_matrix.T) / 2
    distance_matrix = np.maximum(distance_matrix, 0)

    condensed_dist = squareform(distance_matrix)

    # Hierarchical clustering
    linkage_matrix = linkage(condensed_dist, method="ward")

    # Cut into clusters
    clusters = fcluster(linkage_matrix, k_clusters, criterion="maxclust")

    # Add to dataframe
    result = seq_strings.copy().reset_index(drop=True)
    result["cluster"] = clusters

    # Analyze clusters
    cluster_summary = result.groupby("cluster").agg({
        "sequence_id": "count",
        "n_actions": "mean",
        "ends_shot": "mean",
        "total_xg": "mean",
        "action_string": "first"
    }).rename(columns={
        "sequence_id": "count",
        "n_actions": "avg_actions",
        "ends_shot": "shot_rate",
        "total_xg": "avg_xg",
        "action_string": "example_sequence"
    }).reset_index().sort_values("shot_rate", ascending=False)

    # Find cluster archetypes
    archetypes = []
    for cl in range(1, k_clusters + 1):
        cluster_mask = clusters == cl
        cluster_indices = np.where(cluster_mask)[0]

        if len(cluster_indices) <= 1:
            archetype_idx = cluster_indices[0] if len(cluster_indices) > 0 else 0
            centrality = 1.0
        else:
            cluster_sims = similarity_matrix[np.ix_(cluster_indices, cluster_indices)]
            centrality_scores = cluster_sims.mean(axis=1)
            most_central = cluster_indices[np.argmax(centrality_scores)]
            archetype_idx = most_central
            centrality = centrality_scores.max()

        arch = result.iloc[archetype_idx].copy()
        arch["centrality"] = centrality
        archetypes.append(arch)

    archetypes_df = pd.DataFrame(archetypes)

    return result, cluster_summary, archetypes_df, linkage_matrix

# Cluster sequences
result, summary, archetypes, linkage_mat = cluster_similar_sequences(
    seq_strings.reset_index(drop=True), sim_matrix, k_clusters=6
)

print("Sequence Cluster Summary:")
print(summary)

print("\nCluster Archetypes:")
print(archetypes[["cluster", "action_string", "n_actions", "ends_shot", "centrality"]])

# Python: Analyzing Temporal Dynamics
import pandas as pd
import numpy as np

def analyze_temporal_dynamics(events):
    """Analyze speed and rhythm of sequences."""
    events = events.sort_values(["sequence_id", "index"]).copy()
    events["time_seconds"] = events["minute"] * 60 + events["second"]

    def compute_temporal_features(group):
        times = group["time_seconds"].values
        gaps = np.diff(times)

        if len(gaps) == 0:
            return pd.Series({
                "team": group["team"].iloc[0],
                "total_duration": 0,
                "n_events": len(group),
                "avg_gap": 0,
                "median_gap": 0,
                "gap_sd": 0,
                "first_half_avg_gap": 0,
                "second_half_avg_gap": 0,
                "ends_shot": (group["type"] == "Shot").any(),
                "total_xg": group.loc[group["type"] == "Shot", "shot_statsbomb_xg"].sum()
            })

        half = len(gaps) // 2
        first_half_gaps = gaps[:max(half, 1)]
        second_half_gaps = gaps[half:] if half < len(gaps) else gaps

        return pd.Series({
            "team": group["team"].iloc[0],
            "total_duration": times[-1] - times[0],
            "n_events": len(group),
            "avg_gap": gaps.mean(),
            "median_gap": np.median(gaps),
            "gap_sd": gaps.std() if len(gaps) > 1 else 0,
            "max_gap": gaps.max(),
            "min_gap": gaps[gaps > 0].min() if (gaps > 0).any() else 0,
            "first_half_avg_gap": first_half_gaps.mean(),
            "second_half_avg_gap": second_half_gaps.mean(),
            "ends_shot": (group["type"] == "Shot").any(),
            "total_xg": group.loc[group["type"] == "Shot", "shot_statsbomb_xg"].sum()
        })

    temporal = events.groupby("sequence_id").apply(compute_temporal_features).reset_index()

    # Tempo change
    temporal["tempo_change"] = temporal["second_half_avg_gap"] - temporal["first_half_avg_gap"]

    # Classify rhythm
    def classify_rhythm(sd):
        if sd < 1: return "Steady"
        if sd < 2: return "Variable"
        return "Erratic"

    temporal["rhythm"] = temporal["gap_sd"].apply(classify_rhythm)

    # Classify tempo
    def classify_tempo(avg):
        if avg < 2: return "Very Fast"
        if avg < 3: return "Fast"
        if avg < 5: return "Medium"
        if avg < 8: return "Slow"
        return "Very Slow"

    temporal["tempo"] = temporal["avg_gap"].apply(classify_tempo)

    return temporal

temporal = analyze_temporal_dynamics(events)

# Effectiveness by tempo and rhythm
effectiveness = temporal.groupby(["tempo", "rhythm"]).agg({
    "sequence_id": "count",
    "ends_shot": "mean",
    "total_xg": "mean"
}).rename(columns={"sequence_id": "count"}).reset_index()

effectiveness = effectiveness[effectiveness["count"] >= 10].sort_values("total_xg", ascending=False)
print("Effectiveness by Tempo and Rhythm:")
print(effectiveness)

# Tempo acceleration effect
def classify_tempo_pattern(change):
    if change < -1: return "Accelerating"
    if change > 1: return "Decelerating"
    return "Steady"

temporal["tempo_pattern"] = temporal["tempo_change"].apply(classify_tempo_pattern)

acceleration_effect = temporal.groupby("tempo_pattern").agg({
    "sequence_id": "count",
    "ends_shot": "mean",
    "total_xg": "mean"
}).rename(columns={"sequence_id": "count"}).reset_index()

print("\nEffect of Tempo Changes:")
print(acceleration_effect)

# Python: Comprehensive Team Sequence Profile
import pandas as pd
import numpy as np

def create_team_sequence_profile(events, team_name):
    """Create comprehensive sequence profile for a team."""
    team_events = events[events["team"] == team_name].copy()
    team_events = team_events.sort_values("index")

    # Detect sequences
    team_events["team_prev"] = team_events["team"].shift(1)
    team_events["possession_change"] = (
        (team_events["team"] != team_events["team_prev"]) |
        team_events["type"].isin(["Half Start", "Starting XI", "Kick Off"])
    )
    team_events["sequence_id"] = team_events["possession_change"].cumsum()

    # Calculate sequence metrics
    def compute_seq(group):
        times = group["minute"] * 60 + group["second"]
        locs_x = group["location_x"].dropna()
        locs_y = group["location_y"].dropna()

        return pd.Series({
            "duration": times.max() - times.min() if len(times) > 1 else 0,
            "n_passes": (group["type"] == "Pass").sum(),
            "n_events": len(group),
            "start_x": locs_x.iloc[0] if len(locs_x) > 0 else 60,
            "end_x": locs_x.iloc[-1] if len(locs_x) > 0 else 60,
            "max_x": locs_x.max() if len(locs_x) > 0 else 60,
            "width_used": (locs_y.max() - locs_y.min()) if len(locs_y) > 0 else 0,
            "net_progress": (locs_x.iloc[-1] - locs_x.iloc[0]) if len(locs_x) > 1 else 0,
            "ends_shot": (group["type"] == "Shot").any(),
            "total_xg": group.loc[group["type"] == "Shot", "shot_statsbomb_xg"].sum()
        })

    sequences = team_events.groupby("sequence_id").apply(compute_seq).reset_index()

    # Profile summary
    profile = {
        "team": team_name,
        "total_sequences": len(sequences),
        "avg_passes": sequences["n_passes"].mean(),
        "median_passes": sequences["n_passes"].median(),
        "pct_short": (sequences["n_passes"] <= 3).mean(),
        "pct_long": (sequences["n_passes"] >= 10).mean(),
        "avg_duration": sequences["duration"].mean(),
        "pct_fast": (sequences["duration"] < 10).mean(),
        "avg_progress": sequences["net_progress"].mean(),
        "avg_width": sequences["width_used"].mean(),
        "pct_reach_box": (sequences["max_x"] >= 102).mean(),
        "shot_rate": sequences["ends_shot"].mean(),
        "avg_xg_per_seq": sequences["total_xg"].mean(),
        "xg_per_shot": (sequences["total_xg"].sum() /
                       max(sequences["ends_shot"].sum(), 1))
    }

    return profile, sequences

# Compare two teams
profile1, seq1 = create_team_sequence_profile(events, "Barcelona")
profile2, seq2 = create_team_sequence_profile(events, "Real Madrid")

# Create comparison
metrics = ["total_sequences", "avg_passes", "pct_short", "pct_long",
           "avg_duration", "pct_fast", "avg_progress", "avg_width",
           "pct_reach_box", "shot_rate", "avg_xg_per_seq", "xg_per_shot"]

labels = ["Total Sequences", "Avg Passes/Seq", "% Short (≤3)",
          "% Long (≥10)", "Avg Duration (s)", "% Fast (<10s)",
          "Avg X Progress", "Avg Width Used", "% Reach Box",
          "Shot Rate", "xG per Sequence", "xG per Shot"]

comparison = pd.DataFrame({
    "Metric": labels,
    "Barcelona": [profile1[m] for m in metrics],
    "Real Madrid": [profile2[m] for m in metrics]
})

# Format percentages
for col in ["Barcelona", "Real Madrid"]:
    for i, m in enumerate(metrics):
        if "pct" in m or m == "shot_rate":
            comparison.loc[i, col] = f"{comparison.loc[i, col] * 100:.1f}%"
        elif m in ["avg_xg_per_seq"]:
            comparison.loc[i, col] = f"{comparison.loc[i, col]:.4f}"
        elif m == "xg_per_shot":
            comparison.loc[i, col] = f"{comparison.loc[i, col]:.3f}"
        else:
            comparison.loc[i, col] = f"{comparison.loc[i, col]:.2f}"

print(comparison.to_string(index=False))