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')

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

    return P

P = calculate_transition_probabilities(mdp_data, num_states, num_actions)

def create_feature_matrix(num_states):
    # Define meaningful features for your states
    # For simplicity, using identity matrix here
    return np.eye(num_states)

feature_matrix = create_feature_matrix(num_states)

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

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

def compute_empirical_feature_counts(trajectories, feature_matrix):
    num_features = feature_matrix.shape[1]
    mu_E = np.zeros(num_features)
    for traj in trajectories:
        for s, _ in traj:
            mu_E += feature_matrix[int(s)]
    mu_E /= len(trajectories)
    return mu_E

mu_E = compute_empirical_feature_counts(trajectories, feature_matrix)

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

    # Initialize theta
    theta = np.random.uniform(size=(num_features,))

    for iteration in range(num_iterations):
        # 1. Compute the reward function
        R = feature_matrix.dot(theta)

        # 2. Compute the state-value function V(s) using log-sum-exp
        V = value_iteration(P, R, num_states, num_actions, gamma)

        # 3. Compute the stochastic policy π(s, a)
        pi = compute_stochastic_policy(P, R, V, num_states, num_actions, gamma)

        # 4. Compute expected state visitation frequencies under the current policy
        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()
        expected_svf = compute_expected_state_visitation_freq(P, pi, start_prob, num_states, num_actions, trajectory_length)

        # 5. Compute expected feature counts under the model
        mu_theta = feature_matrix.T.dot(expected_svf)

        # 6. Compute the gradient
        grad = mu_E - mu_theta

        # 7. Update theta
        theta += learning_rate * grad

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

    return theta

def value_iteration(P, R, num_states, num_actions, gamma=0.99):
    V = np.zeros(num_states)
    delta = float('inf')
    threshold = 1e-4

    while delta > threshold:
        V_prev = V.copy()
        for s in range(num_states):
            Q = np.zeros(num_actions)
            for a in range(num_actions):
                Q[a] = R[s] + gamma * np.sum(P[s, a, :] * V_prev)
            V[s] = logsumexp(Q)
        delta = np.max(np.abs(V - V_prev))
    return V

def compute_stochastic_policy(P, R, V, num_states, num_actions, gamma=0.99):
    pi = np.zeros((num_states, num_actions))
    for s in range(num_states):
        Q = np.zeros(num_actions)
        for a in range(num_actions):
            Q[a] = R[s] + gamma * np.sum(P[s, a, :] * V)
        log_pi = Q - V[s]  # Since V[s] = logsumexp(Q)
        pi[s, :] = np.exp(log_pi)
    # pi is already normalized due to logsumexp in V[s]
    return pi

def compute_expected_state_visitation_freq(P, pi, start_prob, num_states, num_actions, trajectory_length):
    D = np.zeros((trajectory_length, num_states))
    D[0, :] = start_prob
    for t in range(1, trajectory_length):
        for s in range(num_states):
            for a in range(num_actions):
                D[t, :] += D[t - 1, s] * pi[s, a] * P[s, a, :]
    expected_svf = D.sum(axis=0)
    return expected_svf

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

# Compute the reward function
reward_function = feature_matrix.dot(theta)

# Plot the learned reward function
plt.bar(range(num_states), reward_function)
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)
    for trajectory in trajectories:
        for state, _ in trajectory:
            svf[int(state)] += 1
    svf /= svf.sum()
    return svf

svf = compute_state_visitation_freq(trajectories, num_states)

# Plot state visitation frequencies
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.scatter(reward_function, 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(reward_function, 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, theta=1e-6):
    V = np.zeros(num_states)
    pi = np.zeros(num_states, dtype=int)
    delta = float('inf')
    while delta > theta:
        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] = np.sum(P[s, a, :] * (rewards + gamma * 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, reward_function, num_states, num_actions)

# Plot the optimal policy
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 kod Ziebart'ın MaxEnt IRL makalesine uygun olacak şekilde hazırlandı. Makaleye olan uygunluğunu kontrol edip hatalı yerleri düzeltir misin?