import numpy as np
import pandas as pd
from scipy.special import logsumexp
import matplotlib.pyplot as plt

def load_mdp_data(file_path):
    return pd.read_csv(file_path)

# Load your MDP data
mdp_data = load_mdp_data('/content/drive/MyDrive/BIST_10/MDP_IRL/MLB/AKBNK.csv')

# Define the number of states and actions based on your data
num_states = 15  # Update based on your data
num_actions = 3  # Update based on your data

def calculate_transition_probabilities(mdp_data, num_states, num_actions):
    P = np.zeros((num_states, num_actions, num_states))
    for _, row in mdp_data.iterrows():
        s, a, s_next = int(row['State']), int(row['Action']), int(row['Next_State'])
        P[s, a, s_next] += 1

    # Normalize
    for s in range(num_states):
        for a in range(num_actions):
            total = np.sum(P[s, a, :])
            if total > 0:
                P[s, a, :] /= total
            else:
                # If no transitions for this state-action pair, make it uniform
                P[s, a, :] = 1.0 / num_states

    return P

P = calculate_transition_probabilities(mdp_data, num_states, num_actions)

def create_feature_matrix(num_states, num_actions):
    # Create state-action feature matrix: one-hot encoding of state-action pairs
    num_features = num_states * num_actions
    feature_matrix = np.identity(num_features)
    return feature_matrix

feature_matrix = create_feature_matrix(num_states, num_actions)

def generate_trajectories(mdp_data, num_trajectories, trajectory_length):
    trajectories = []
    data_length = len(mdp_data)
    for _ in range(num_trajectories):
        start_idx = np.random.randint(0, data_length - trajectory_length)
        trajectory = mdp_data.iloc[start_idx:start_idx+trajectory_length][['State', 'Action']].values.tolist()
        trajectories.append(trajectory)
    return trajectories

num_trajectories = 20
trajectory_length = 60
trajectories = generate_trajectories(mdp_data, num_trajectories, trajectory_length)

def compute_empirical_feature_counts(trajectories, num_states, num_actions):
    num_features = num_states * num_actions
    mu_E = np.zeros(num_features)
    total_timesteps = 0
    for traj in trajectories:
        for s, a in traj:
            idx = s * num_actions + a
            mu_E[idx] += 1  # Feature is one-hot on state-action pair
            total_timesteps += 1
    mu_E /= total_timesteps
    return mu_E

mu_E = compute_empirical_feature_counts(trajectories, num_states, num_actions)

def maxent_irl(feature_matrix, P, mu_E, num_states, num_actions, trajectories, gamma=0.99,
               learning_rate=0.1, num_iterations=100, lambda_reg=0.1):
    num_features = feature_matrix.shape[1]

    # Initialize theta
    theta = np.zeros(num_features)

    # Precompute start state probabilities from trajectories
    start_state_counts = np.zeros(num_states)
    for traj in trajectories:
        start_state_counts[int(traj[0][0])] += 1
    start_prob = start_state_counts / start_state_counts.sum()

    for iteration in range(num_iterations):
        # 1. Compute the reward function
        rewards = feature_matrix.dot(theta)  # Size: num_states * num_actions

        # 2. Compute the policy using soft value iteration
        policy = soft_value_iteration(P, rewards, num_states, num_actions, gamma)

        # 3. Compute expected feature counts under the model
        mu_theta = compute_expected_feature_counts(policy, P, num_states, num_actions, start_prob, gamma)

        # 4. Compute the gradient
        grad = mu_E - mu_theta - lambda_reg * theta

        # 5. Update theta
        theta += learning_rate * grad

        # Optionally, print the gradient norm to monitor convergence
        print(f"Iteration {iteration+1}, Gradient Norm: {np.linalg.norm(grad)}")

    return theta

def soft_value_iteration(P, rewards, num_states, num_actions, gamma=0.99, threshold=1e-6):
    # Rewards should be of shape (num_states, num_actions)
    rewards = rewards.reshape((num_states, num_actions))

    # Initialize value function
    V = np.zeros(num_states)

    delta = float('inf')
    while delta > threshold:
        V_prev = V.copy()

        # Compute Q-values
        Q = np.zeros((num_states, num_actions))
        for a in range(num_actions):
            Q[:, a] = rewards[:, a] + gamma * P[:, a, :].dot(V_prev)

        # Update value function V(s) using log-sum-exp for numerical stability
        V = logsumexp(Q, axis=1)

        delta = np.max(np.abs(V - V_prev))

    # Compute policy: π(a|s) ∝ exp(Q(s, a) - V(s))
    policy = np.exp(Q - V[:, None])

    return policy

def compute_expected_feature_counts(policy, P, num_states, num_actions, start_prob, gamma=0.99, threshold=1e-6):
    # Compute state transition matrix under policy
    P_pi = np.zeros((num_states, num_states))
    for s in range(num_states):
        for a in range(num_actions):
            P_pi[s, :] += policy[s, a] * P[s, a, :]

    # Compute expected state visitation frequencies using fixed-point iteration
    D_s = np.copy(start_prob)
    D_s = D_s * (1 - gamma)
    delta = float('inf')
    while delta > threshold:
        D_s_prev = D_s.copy()
        D_s = (1 - gamma) * start_prob + gamma * P_pi.T.dot(D_s_prev)
        delta = np.max(np.abs(D_s - D_s_prev))

    # Compute expected feature counts: μ = ∑_s D_s[s] ∑_a π(a|s) * f(s,a)
    mu_theta = np.zeros(num_states * num_actions)
    for s in range(num_states):
        for a in range(num_actions):
            idx = s * num_actions + a
            mu_theta[idx] = D_s[s] * policy[s, a]

    return mu_theta

# Run MaxEnt IRL to get theta
theta = maxent_irl(feature_matrix, P, mu_E, num_states, num_actions, trajectories)

# Compute the reward function
rewards = feature_matrix.dot(theta)
rewards = rewards.reshape((num_states, num_actions))
# For visualization, compute state rewards by averaging over actions
state_rewards = rewards.mean(axis=1)

# Plot the learned reward function
plt.figure(figsize=(10, 6))
plt.bar(range(num_states), state_rewards)
plt.xlabel('State')
plt.ylabel('Reward')
plt.title('Learned Reward Function')
plt.show()

# Compute state visitation frequencies from trajectories
def compute_state_visitation_freq(trajectories, num_states):
    svf = np.zeros(num_states)
    total_timesteps = 0
    for trajectory in trajectories:
        for state, _ in trajectory:
            svf[int(state)] += 1
            total_timesteps += 1
    svf /= total_timesteps
    return svf

svf = compute_state_visitation_freq(trajectories, num_states)

# Plot state visitation frequencies
plt.figure(figsize=(10, 6))
plt.bar(range(num_states), svf)
plt.xlabel('State')
plt.ylabel('Visitation Frequency')
plt.title('State Visitation Frequencies')
plt.show()

# Plot reward vs visitation frequency
plt.figure(figsize=(10, 6))
plt.scatter(state_rewards, svf)
plt.xlabel('Learned Reward')
plt.ylabel('Visitation Frequency')
plt.title('Reward vs Visitation Frequency')
plt.show()

# Compute and print the correlation
correlation = np.corrcoef(state_rewards, svf)[0, 1]
print(f"Correlation between rewards and visitation frequencies: {correlation:.2f}")

# Compute the optimal policy and value function using the learned rewards
def value_iteration_optimal(P, rewards, num_states, num_actions, gamma=0.99, threshold=1e-6):
    V = np.zeros(num_states)
    pi = np.zeros(num_states, dtype=int)
    rewards = rewards.reshape((num_states, num_actions))
    delta = float('inf')
    while delta > threshold:
        V_prev = V.copy()
        for s in range(num_states):
            Q_sa = np.zeros(num_actions)
            for a in range(num_actions):
                Q_sa[a] = rewards[s, a] + gamma * np.sum(P[s, a, :] * V_prev)
            V[s] = np.max(Q_sa)
            pi[s] = np.argmax(Q_sa)
        delta = np.max(np.abs(V - V_prev))
    return pi, V

optimal_policy, V = value_iteration_optimal(P, rewards, num_states, num_actions)

# Plot the optimal policy
plt.figure(figsize=(10, 6))
plt.bar(range(num_states), optimal_policy)
plt.xlabel('State')
plt.ylabel('Optimal Action')
plt.title('Optimal Policy')
plt.yticks(range(num_actions))
plt.show()

def policy_agreement(trajectories, optimal_policy):
    agreement = 0
    total = 0
    for trajectory in trajectories:
        for state, action in trajectory:
            if int(action) == optimal_policy[int(state)]:
                agreement += 1
            total += 1
    return agreement / total

agreement_ratio = policy_agreement(trajectories, optimal_policy)
print(f"Agreement between optimal policy and demonstrated trajectories: {agreement_ratio:.2f}")

Yukarıdaki kodun Ziebart'ın MaxEnt IRL makalesine uygun halde olup olmadığını kontrol edebilir misin?