import numpy as np
import pandas as pd
from scipy.optimize import minimize

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

mdp_data = load_mdp_data('/content/drive/MyDrive/BIST_10/MDP_IRL/MLB/AKBNK.csv')

num_states = 15
num_actions = 3

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 = row['State'], row['Action'], 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):
    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 maxent_irl(feature_matrix, P, trajectories, learning_rate=0.1, num_iterations=100):
    num_states, num_features = feature_matrix.shape
    num_actions = P.shape[1]

    # Initialize theta randomly
    theta = np.random.rand(num_features)

    for _ in range(num_iterations):
        # Forward pass: calculate state visitation frequencies
        D = np.zeros((num_states, num_actions, num_states))
        for trajectory in trajectories:
            for (s, a), (next_s, _) in zip(trajectory[:-1], trajectory[1:]):
                D[s, a, next_s] += 1

        # Backward pass: compute the gradient
        expected_svf = np.zeros(num_states)
        for s in range(num_states):
            for a in range(num_actions):
                expected_svf += D[s, a].sum() * P[s, a]

        grad = feature_matrix.T.dot(expected_svf)

        # Extract initial states from trajectories
        initial_states = [traj[0][0] for traj in trajectories]
        grad -= feature_matrix[initial_states].sum(axis=0)

        # Update theta
        theta += learning_rate * grad

    return theta

# Example usage:
# Assuming 'feature_matrix', 'P', and 'trajectories' are defined appropriately
theta = maxent_irl(feature_matrix, P, trajectories)

reward_function = feature_matrix.dot(theta)

import matplotlib.pyplot as plt

plt.bar(range(num_states), reward_function)
plt.xlabel('State')
plt.ylabel('Reward')
plt.title('Learned Reward Function')
plt.show()

def compute_state_visitation_freq(trajectories, num_states):
    svf = np.zeros(num_states)
    for trajectory in trajectories:
        for state, _ in trajectory:
            svf[state] += 1
    svf /= svf.sum()
    return svf

svf = compute_state_visitation_freq(trajectories, num_states)

plt.bar(range(num_states), svf)
plt.xlabel('State')
plt.ylabel('Visitation Frequency')
plt.title('State Visitation Frequencies')
plt.show()

plt.scatter(reward_function, svf)
plt.xlabel('Learned Reward')
plt.ylabel('Visitation Frequency')
plt.title('Reward vs Visitation Frequency')
plt.show()

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

def value_iteration(P, rewards, gamma=0.99, theta=1e-6):
    V = np.zeros(num_states)
    while True:
        delta = 0
        for s in range(num_states):
            v = V[s]
            V[s] = max([sum([P[s, a, s1] * (rewards[s1] + gamma * V[s1]) for s1 in range(num_states)]) for a in range(num_actions)])
            delta = max(delta, abs(v - V[s]))
        if delta < theta:
            break

    # Compute the policy
    policy = np.zeros(num_states, dtype=int)
    for s in range(num_states):
        policy[s] = np.argmax([sum([P[s, a, s1] * (rewards[s1] + gamma * V[s1]) for s1 in range(num_states)]) for a in range(num_actions)])

    return policy, V

optimal_policy, V = value_iteration(P, reward_function)

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 action == optimal_policy[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}")

Kodu detaylı bir şekilde kontrol edebilir misin? Düzeltmesi ve geliştirilmesi gereken yerleri belirt. Kod ileZiebart'ın MaxEnt makalesi uygulanmaya çalışıldı.