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
            else:
                # If no transitions for this state-action pair, make it a uniform distribution
                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):
    # 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 = []
    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, feature_matrix):
    num_features = feature_matrix.shape[1]
    mu_E = np.zeros(num_features)
    total_timesteps = 0
    for traj in trajectories:
        for s, _ in traj:
            mu_E += feature_matrix[int(s)]
            total_timesteps += 1
    mu_E /= total_timesteps
    return mu_E

mu_E = compute_empirical_feature_counts(trajectories, feature_matrix)

def maxent_irl(feature_matrix, P, mu_E, num_states, num_actions, trajectories, trajectory_length,
               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,))

    # 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
        R = feature_matrix.dot(theta)

        # 2. Compute the state partition function (logZ) via backward pass
        V = compute_state_partition_function(P, R, num_states, num_actions, gamma)

        # 3. Compute the policy
        pi = compute_policy(P, R, V, num_states, num_actions, gamma)

        # 4. Compute expected state visitation frequencies under the current policy
        expected_svf = compute_expected_state_visitation_freq(P, pi, start_prob, num_states, num_actions, trajectory_length, gamma)

        # 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+1}, Gradient Norm: {np.linalg.norm(grad)}")

    return theta

def compute_state_partition_function(P, R, num_states, num_actions, gamma):
    V = np.zeros(num_states)
    V_old = np.zeros(num_states)
    max_iters = 1000
    tol = 1e-6
    for iteration in range(max_iters):
        V_old[:] = V[:]
        for s in range(num_states):
            Vs = []
            for a in range(num_actions):
                Qsa = R[s]
                Qsa += gamma * np.sum(P[s, a, :] * V_old)
                Vs.append(Qsa)
            V[s] = logsumexp(Vs)
        if np.max(np.abs(V - V_old)) < tol:
            break
    return V

def compute_policy(P, R, V, num_states, num_actions, gamma):
    pi = np.zeros((num_states, num_actions))
    for s in range(num_states):
        Qa = np.zeros(num_actions)
        for a in range(num_actions):
            Qa[a] = R[s] + gamma * np.sum(P[s, a, :] * V)
        log_pi_s = Qa - V[s]
        pi[s, :] = np.exp(log_pi_s)
        # Normalize to ensure sum over actions equals 1
        pi[s, :] /= np.sum(pi[s, :])
    return pi

def compute_expected_state_visitation_freq(P, pi, start_prob, num_states, num_actions, trajectory_length, gamma):
    # Initialize state visitation frequencies
    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 = np.sum(D, axis=0)
    return expected_svf

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

# 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)
    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.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, threshold=1e-6):
    V = np.zeros(num_states)
    pi = np.zeros(num_states, dtype=int)
    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] + 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, 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 kodun Ziebart'ın MaxEnt IRL makalesine olan uygunluğunu en ince ayrıntısına kadar kontrol eder misin? Eksik ya da hatalı yerleri tek bir kod parçası olarak yazmanı istiyorum..