comparison_methods_resource_sharing.py

import math
import time
import copy
import random

import multiprocessing as mp

import os
os.environ["OMP_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"
os.environ["OPENBLAS_NUM_THREADS"] = "1"
os.environ["NUMEXPR_NUM_THREADS"] = "1"
os.environ["VECLIB_MAXIMUM_THREADS"] = "1"
os.environ["BLIS_NUM_THREADS"] = "1"
os.environ["PYTORCH_NUM_THREADS"] = "1"

import numpy as np
import matplotlib.pyplot as plt


from _mdp_single_device import SingleDevice
from _mdp_centralized_agent import CentralizedSystem


from _RL_helper import *
from _value_iteration_helper import *
from stable_baselines3 import A2C, PPO

import pickle as pk


def stochastic_policy_evaluation(list_agents, list_policies, n_episodes = 100):

    global_discounted_rewards = np.zeros((n_episodes, ))

    discounted_returns = np.zeros((n_episodes, len(list_agents)))
    discounted_approximated_rewards = np.zeros((n_episodes, len(list_agents)))
    discounted_costs = np.zeros((n_episodes, len(list_agents)))

    average_returns = np.zeros((n_episodes, len(list_agents)))
    average_approximated_rewards = np.zeros((n_episodes, len(list_agents)))
    average_costs = np.zeros((n_episodes, len(list_agents)))

    centralized_system = CentralizedSystem(list_agents)

    for episode in range(n_episodes):
        state = dict()
        for agent in list_agents:
            state['Agent {}'.format(agent.index)] = agent.reset()
        
        global_discounted_reward = 0
        discounted_reward = np.zeros((len(list_agents), ))
        discounted_approximated_reward = np.zeros((len(list_agents), ))
        
        global_discounted_cost = 0
        discounted_cost = np.zeros((len(list_agents), ))
        
        global_average_reward = 0
        average_reward = np.zeros((len(list_agents), ))
        average_approximated_reward = np.zeros((len(list_agents), ))
        
        global_average_cost = 0
        average_cost = np.zeros((len(list_agents), ))

        for t in range(list_agents[0].max_episodes_steps):
            global_action = np.zeros((len(list_agents), ), dtype=int)
            cost_vector = np.zeros((len(list_agents), ))
            next_state = dict()
            for agent in list_agents:
                agent_state = state['Agent {}'.format(agent.index)]
                index = agent.compute_state_index(agent_state['x'], agent_state['e'])
                # choose the values from 1 to 3 with the probabilities given by the policy
                global_action[agent.index - 1] = random.choices(range(agent.action_space_size), weights=list_policies[agent.index-1][index])[0]
                # we could probably get the approximated reward for each agent from here and save a step...
                next_state_single_agent, individual_reward, _, info = agent.step(state['Agent {}'.format(agent.index)], global_action[agent.index -1], training=False)
                discounted_approximated_reward[agent.index - 1] += agent.gamma**t * individual_reward
                average_approximated_reward[agent.index - 1] += individual_reward
                next_state['Agent {}'.format(agent.index)] = next_state_single_agent
                cost_vector[agent.index - 1] = info['cost']
            # global evaluation
            global_reward, global_reward_vector = centralized_system.global_reward(global_action, state = state, vector_rewards=True)

            global_discounted_reward += list_agents[0].gamma**t * global_reward
            discounted_reward += list_agents[0].gamma**t * global_reward_vector

            global_discounted_cost += list_agents[0].gamma**t * np.sum(cost_vector)
            discounted_cost += list_agents[0].gamma**t * cost_vector
            
            global_average_reward += global_reward_vector
            average_reward += global_reward_vector

            global_average_cost += cost_vector
            average_cost += cost_vector


            state = next_state
        global_discounted_rewards[episode] = global_discounted_reward
        discounted_returns[episode] = discounted_reward
        discounted_approximated_rewards[episode] = discounted_approximated_reward
        discounted_costs[episode] = discounted_cost

        average_returns[episode] = average_reward / (t + 1)
        average_approximated_rewards[episode] = average_approximated_reward / (t + 1)
        average_costs[episode] = average_cost / (t + 1)

    return np.mean(discounted_returns, axis = 0), np.mean(discounted_approximated_rewards, axis = 0), np.mean(discounted_costs, axis = 0), np.mean(average_returns, axis = 0), np.mean(average_approximated_rewards, axis = 0), np.mean(average_costs, axis = 0)



# def single_agent_policy_learning(list_agents, algorithm, noisy_average_constraint_list = None, noisy_discounted_constraint_list = None, noise_position = None):
def single_agent_policy_learning(list_agents, algorithm, local_noise = None, other_noise = None, noise_position = None):
    
    if algorithm != 'RL' and algorithm != 'VI':
        raise ValueError("Method not implemented yet")

    possible_noise_positions = ['local', 'other', 'none']
    if noise_position not in possible_noise_positions:
        raise ValueError("Noise position not implemented yet")
    
    local_optimized_discounted_rewards = np.zeros((n_players, ))
    local_optimized_discounted_costs = np.zeros((n_players, ))

    local_optimized_average_rewards = np.zeros((n_players, ))
    local_optimized_average_costs = np.zeros((n_players, ))
    policy_list = []


    for agent in list_agents: 

        # we copy here the local constraints of the agent, so that we can restore them after the optimization
        agent_local_average_constraint_copy = copy.copy(agent.average_constraint)
        agent_local_discounted_constraint_copy = copy.copy(agent.discounted_constraint)
        agent_local_sum_average_constraints_copy = copy.copy(agent.sum_average_constraints)

        agent.sum_average_constraints = np.sum([other_agent.average_constraint for other_agent in list_agents]) - agent.average_constraint # sum_average_constraints

        if noise_position == 'local':
            # agent.average_constraint += local_noise[list_agents.index(agent)] # noisy_average_constraint_list[list_agents.index(agent)]
            # agent.discounted_constraint = agent.average_constraint * (1 - agent.gamma**(agent.max_episodes_steps+1))/(1 - agent.gamma)
            # we update the sum of constraints, which will be used by the approximation of the immediate reward function 
            agent.sum_average_constraints = agent.sum_average_constraints # - agent_local_average_constraint_copy + agent.average_constraint
        elif noise_position == 'other':
            # we only update the sum of the other values, with the noisy values (except for the local value, which is not noisy)
            agent.sum_average_constraints = agent.sum_average_constraints + other_noise[list_agents.index(agent)] # - noisy_average_constraint_list[list_agents.index(agent)] + agent.average_constraint
        # else:
        #     agent.average_constraint = average_constraint_list[list_agents.index(agent)]
        #     agent.discounted_constraint = discounted_constraint_list[list_agents.index(agent)]

        # we now update the sum of the average constraints
        # if only_local_value:
        #     agent.sum_average_constraints = np.sum(average_constraint_list) - average_constraint_list[list_agents.index(agent)] #+ agent.average_constraint
        # else:
        #     agent.sum_average_constraints = np.sum(noisy_average_constraint_list) - noisy_average_constraint_list[list_agents.index(agent)] # + agent.average_constraint

        # print("Agent = {}. Average constraint = {}. Discounted constraint = {}. Sum other constraints = {}".format(list_agents.index(agent)+1, agent.average_constraint, agent.discounted_constraint, sum_average_constraints))
        
        # optimization of Lagrange multiplier for the agent
        agent.lagrange_multiplier = 0
        LM_agent_optimized = False
        LM_optimization_steps = 1
        LM_optimization_total_steps = 0
        constraint_satisfied = False
        consecutive_constraint_not_satisfied = 0

        LM_optimization_learning_rate = 10
        old_LM = math.inf
        q_table = None

        min_distance_LM = math.inf

        best_policy = None
        
        best_discounted_policy_reward = None
        best_discounted_penalized_reward = None
        best_average_policy_reward = None

        best_discounted_policy_cost = None
        best_average_policy_cost = None

        while not LM_agent_optimized and LM_optimization_steps < 50 and LM_optimization_total_steps < 50:
            # in this case we can simply apply the Q-learning policy learning (this is possible because it is a simple environment)
            if algorithm == 'RL':
                agent_optimal_solution = q_learning(env = agent, q_table = q_table, learning_rate = .5, exploration_rate = .05, exploration_decay = 0.95, num_episodes = 20, episodes_with_no_exploration_rate = 0)
                q_table = agent_optimal_solution[0]
            elif algorithm == 'LP':
                policy_list = compute_optimal_solution(env = agent, cost_distribution = None, harvesting_distribution = None, return_v_estimate = False, n_iterations=10)
            # print('Agent {}. Optimal policy found'.format(list_agents.index(agent)+1))
            

            # we save a measure of the difference from he saturation of the constraint 

            # evaluation of the policy obtained (stochastic evaluation)
            if algorithm == 'RL':
                # single_agent_evaluation = evaluate_single_agent_policy(agent, agent_optimal_solution[1])
                discounted_policy_reward = agent_optimal_solution[2]
                discounted_policy_cost = agent_optimal_solution[3]
                discounted_penalized_reward = agent_optimal_solution[6]

                average_policy_reward = agent_optimal_solution[4]
                average_policy_cost = agent_optimal_solution[5]

                if noise_position == 'local':
                    error_constraint = np.sum(np.abs(average_policy_cost - (agent.average_constraint + local_noise[list_agents.index(agent)]))) / n_players
                else:
                    error_constraint = np.sum(np.abs(average_policy_cost - agent.average_constraint)) / n_players


            elif algorithm == 'LP':
                discounted_policy_reward = evaluate_policy(agent, agent_optimal_solution[1])

            # if math.fabs(average_policy_cost - agent.average_constraint) <= min_distance_LM:
            min_distance_LM = math.fabs(average_policy_cost - agent.average_constraint)
            best_policy = agent_optimal_solution[1]
            best_discounted_policy_reward = discounted_policy_reward
            best_discounted_penalized_reward = discounted_penalized_reward
            best_average_policy_reward = average_policy_reward
            best_discounted_policy_cost = discounted_policy_cost
            best_average_policy_cost = average_policy_cost
                

            # print(discounted_policy_reward)
            # print("Agent {}. Discounted reward = {}. Discounted cost = {}".format(list_agents.index(agent)+1, discounted_policy_reward[0], discounted_policy_reward[1]))
            # print(agent.lagrange_multiplier, discounted_policy_reward, average_policy_cost)
            
            # change the value of the LM to reflect the results
            if noise_position == 'local':
                new_constraint_satisfied = (average_policy_cost - (agent.average_constraint + local_noise[list_agents.index(agent)]) <= 0)
            else:
                new_constraint_satisfied = (average_policy_cost - agent.average_constraint <= 0)
            LM_optimization_steps += new_constraint_satisfied != constraint_satisfied
            LM_optimization_total_steps += 1

            if new_constraint_satisfied == constraint_satisfied:
                consecutive_constraint_not_satisfied += 1
            else:
                consecutive_constraint_not_satisfied = 0

            # print(new_constraint_satisfied, constraint_satisfied, new_constraint_satisfied != constraint_satisfied, LM_optimization_steps, consecutive_constraint_not_satisfied)
            if noise_position == 'local':
                agent.lagrange_multiplier += (consecutive_constraint_not_satisfied +1 ) *  LM_optimization_learning_rate/LM_optimization_steps * (- (agent.average_constraint +  local_noise[list_agents.index(agent)]) + average_policy_cost)
            else:
                agent.lagrange_multiplier += (consecutive_constraint_not_satisfied +1 ) *  LM_optimization_learning_rate/LM_optimization_steps * (- agent.average_constraint + average_policy_cost)
            agent.lagrange_multiplier = max(0, agent.lagrange_multiplier)
            # evaluate if we have reached the convergence ofr the value of LM
            if abs(agent.lagrange_multiplier - old_LM) < 0.001:
                # print(abs(agent.lagrange_multiplier - old_LM), agent.lagrange_multiplier, old_LM)
                LM_agent_optimized = False
            
            old_LM = agent.lagrange_multiplier
            constraint_satisfied = new_constraint_satisfied

            # print("Agent {}. Constraint = {}. Lagrange multiplier = {}. Discounted reward = {}. Average cost = {}".format(list_agents.index(agent)+1, agent.average_constraint, agent.lagrange_multiplier, discounted_policy_reward, average_policy_cost))
        # print("Agent {}'s optimal policy found. Discounted reward = {}".format(list_agents.index(agent)+1, discounted_policy_reward[0]))
        policy_list.append(best_policy)

        agent.sum_average_constraints = agent_local_sum_average_constraints_copy

        # now we evaluate the final policy
        final_discounted_reward, final_discounted_cost, final_average_reward, final_average_cost = evaluate_single_agent_policy(agent, agent_optimal_solution[1], num_episodes = 100)

        local_optimized_discounted_rewards[list_agents.index(agent)] = final_discounted_reward
        local_optimized_discounted_costs[list_agents.index(agent)] = final_discounted_cost

        local_optimized_average_rewards[list_agents.index(agent)] = final_average_reward
        local_optimized_average_costs[list_agents.index(agent)] = final_average_cost
        print("Agent {}. Lagrange multiplier = {}. Discounted reward = {}. Discounted penalized reward = {}. Average cost = {}".format(list_agents.index(agent)+1, agent.lagrange_multiplier, best_discounted_policy_reward, best_discounted_penalized_reward, best_average_policy_cost))

    error_constraint = 0
    for i in range(n_players):
        error_constraint += np.abs(local_optimized_average_costs[i]/list_agents[i].average_constraint - 1) / n_players

    print("Local optimized average costs = ", local_optimized_average_costs)
    print("Error constraint = ", error_constraint)


    return local_optimized_discounted_rewards, policy_list, local_optimized_discounted_costs, local_optimized_discounted_rewards, local_optimized_average_costs

def global_evaluation(list_agents, policy_list, n_evaluation_episodes = 100, algorithm = None):
    # function that evaluates the whole environment (i.e. takes into account the interactions between the agents)

    reward_evaluations = np.zeros((n_evaluation_episodes,))
    # cost_evaluations = np.zeros((len(list_agents),))

    agents_action_interaction = np.zeros((n_players, ))
    for episode_index in range(n_evaluation_episodes):
        state = []
        for agent in list_agents:
            state.append(agent.reset(training = False))
        done = False
        reward_sum = 0
        cost_sum = 0
        for t in range(list_agents[0].max_episodes_steps):
            next_state = []
            global_action = np.zeros((n_players, ))
            for agent in list_agents:
                if state[list_agents.index(agent)]['e'] < 0:
                    action = 0
                else:
                    if algorithm == 'LP':
                    # if len(policy_list[list_agents.index(agent)])[state[list_agents.index(agent)]['x']-1, state[list_agents.index(agent)]['e']] >1 :
                        print("State: ", state[list_agents.index(agent)]['x']-1, state[list_agents.index(agent)]['e'])
                        print('Policy: ', policy_list[list_agents.index(agent)][state[list_agents.index(agent)]['x']-1, state[list_agents.index(agent)]['e']])
                        action = np.random.choice([0, 1, 2], p = policy_list[list_agents.index(agent)][state[list_agents.index(agent)]['x']-1, state[list_agents.index(agent)]['e']])
                    else:
                        action = policy_list[list_agents.index(agent)][state[list_agents.index(agent)]['x']-1, state[list_agents.index(agent)]['e']]
                next_state_agent, _, _, info = agent.step(state[list_agents.index(agent)], action, training = False)
                global_action[list_agents.index(agent)] = action
                if action == 2:
                    agents_action_interaction[list_agents.index(agent)] += 1
                next_state.append(next_state_agent)
            global_reward = agent.reward_with_interactions(state, global_action)
            # cost_sum += env.gamma**t * info['cost']
            # print(state, global_action, global_reward)
            state = copy.copy(next_state)
            reward_sum += list_agents[0].gamma**t * global_reward  
        reward_evaluations[episode_index] = reward_sum
        # cost_evaluations[episode_index] = cost_sum 

    # print('Mean interaction action: ', agents_action_interaction/(n_evaluation_episodes * list_agents[0].max_episodes_steps), np.sum(agents_action_interaction)/(n_evaluation_episodes * list_agents[0].max_episodes_steps))

    return np.mean(reward_evaluations)# , np.mean(cost_evaluations)

def global_evaluation_A2C(centralized_env, model, n_evaluation_episodes = 1000, verbose = False):
    # function that evaluates the whole environment (i.e. takes into account the interactions between the agents)

    reward_evaluations = np.zeros((n_evaluation_episodes,))
    # cost_evaluations = np.zeros((len(list_agents),))

    agents_action_interaction = np.zeros((n_players, ))
    frequency_unfeasible_action = 0
    frequency_negative_energy = 0


    for episode_index in range(n_evaluation_episodes):
        centralized_env.reset()
        done = False
        reward_sum = 0
        cost_sum = 0
        for t in range(list_agents[0].max_episodes_steps):
            action = model.predict(centralized_env._get_obs())[0]
            for i in range(n_players):
                if centralized_env._agents_energy[i]<0:
                    frequency_negative_energy += 1
                    if action[i] > 0:
                        frequency_unfeasible_action += 1
                    action[i] = 0
            next_state, global_reward, _, _, info = centralized_env.step(action)
            # print(state, global_action, global_reward)
            # state = copy.copy(next_state)
            reward_sum += list_agents[0].gamma**t * global_reward  
        reward_evaluations[episode_index] = reward_sum
        # cost_evaluations[episode_index] = cost_sum 
    if verbose:
        print("Frequency of unfeasible actions: ", frequency_unfeasible_action / frequency_negative_energy)
    # print("Frequency of negative energy: ", frequency_negative_energy, frequency_unfeasible_action)
    return -np.mean(reward_evaluations)# , np.mean(cost_evaluations)



def single_experiment_wrapper(list_agents):

    if list_agents is None:
        raise ValueError("Missing list agents")
    if cAlg_flag + IQL_flag + A2C_flag + PPO_flag != 1:
        raise ValueError("Too many algorithms")

    print(cAlg_flag, IQL_flag, A2C_flag)

    TOT_STEPS = 1.6e6
    n_episodes = int(TOT_STEPS / (list_agents[0].max_episodes_steps)) + 1
    evaluations_interval = 1000

    if cAlg_flag:
        single_agents_reward_RL, global_reward_RL, final_average_constraint, average_offloading_action_RL = single_experiment_constraint_optimization_with_RL(list_agents, max_iterations = 5)
    else:
        single_agents_reward_RL, global_reward_RL, final_average_constraint, average_offloading_action_RL = [0, 0, 0, 0]

    if IQL_flag:
        single_agents_reward_IQL_individual, global_reward_IQL_individual, average_offloading_action_IQL = single_experiment_IQL_individual_reward(list_agents, n_episodes = n_episodes, episodes_between_evaluation = evaluations_interval, verbose=True)
        if only_individual_IQL:
            single_agents_reward_IQL_common, global_reward_IQL_common = [0, 0] # single_experiment_IQL_common_reward(list_agents, n_episodes = n_episodes, episodes_between_evaluation = evaluations_interval, verbose=True)
        else:
            single_agents_reward_IQL_common, global_reward_IQL_common = single_experiment_IQL_common_reward(list_agents, n_episodes = n_episodes, episodes_between_evaluation = evaluations_interval, verbose=True)

    else:    
        single_agents_reward_IQL_individual, global_reward_IQL_individual, average_offloading_action_IQL = [0, 0, 0] 
        single_agents_reward_IQL_common, global_reward_IQL_common = [0, 0] 
    
    if A2C_flag or PPO_flag:
        global_reward_centralized = single_experiment_centralized_optimization(list_agents, tot_steps = TOT_STEPS/3, episodes_between_evaluation = evaluations_interval)
    else:
        global_reward_centralized = 0

    single_agent_rewards_list = [single_agents_reward_RL, single_agents_reward_IQL_individual, single_agents_reward_IQL_common]
    centralized_list = [global_reward_RL, global_reward_IQL_individual, global_reward_IQL_common, global_reward_centralized]

    return single_agent_rewards_list, centralized_list, final_average_constraint, average_offloading_action_RL, average_offloading_action_IQL


# def single_experiment_constraint_optimization_with_RL(list_agents, constraint_optimization_learning_rate = .25, max_iterations = 10):
def single_experiment_constraint_optimization_with_RL(list_agents, constraint_optimization_learning_rate = .25, max_iterations = 10):

    # initial constraint: [0, ..., 0]
    average_constraint_list = np.zeros((n_players, ))
    # discounted_constraint_list = [average_constraint_list[list_agents.index(agent)] * (1 - agent.gamma**(agent.max_episodes_steps+1))/(1 - agent.gamma) for agent in list_agents]


    constraints_optimized = False
    constraint_optimization_steps = np.ones((n_players, )) 
    positive_gradient = np.zeros((n_players, ))

    new_optimization_steps = 0
    remaining_iterations = max_iterations

    list_objective_function = []
    list_average_offloading_action = []
    list_sum_local_rewards = []

    while not constraints_optimized:
        
        print([agent.average_constraint for agent in list_agents], np.sum([agent.average_constraint for agent in list_agents]))
        # first we have to compute the sum of the average constraint
        # sum_average_constraints = 0
        # for agent in list_agents:
        #     agent.average_constraint = average_constraint_list[list_agents.index(agent)]
        #     sum_average_constraints += agent.average_constraint 
        # for agent in list_agents:
        #     agent.sum_average_constraints = sum_average_constraints


        # 1) COMPUTATION OF THE OPTIMAL LOCAL REWARD FOR EACH AGENT, GIVEN THE CONSTRAINT THETA    
        discounted_reward_no_noise, policy_list, discounted_cost_no_noise, average_reward_no_noise, average_cost_no_noise = single_agent_policy_learning(list_agents = list_agents, algorithm = 'RL', noise_position = 'none')
        print('Local reward with no noise: {}. Sum = {}'.format(discounted_reward_no_noise, np.sum(discounted_reward_no_noise)))
        # evaluation of the global objective function (before we just had the evaluation of the approximated reward)
        global_objective_function = global_evaluation(list_agents, policy_list, algorithm = 'RL')
        print('Global objective function = {}'.format(global_objective_function))
        # average cost for each agent in the evaluation of theta withoute noise
        list_average_offloading_action.append(average_cost_no_noise)
        # 2) DEFINITION OF THE RANDOM NOISE
        if remaining_iterations > 0:
            remaining_iterations -= 1
            variance_local_noise = 0.05
            random_local_noise = np.random.normal(0, variance_local_noise/constraint_optimization_steps, (n_players, ))
            random_other_noise = np.random.normal(0, variance_local_noise * n_players /constraint_optimization_steps, (n_players, ))

            # we want to make sure that the absolute value of random_other_noise is larger than 0.25
            for i in range(n_players):
                if random_other_noise[i] < 0:
                    random_local_noise[i] = -max(0.01, math.fabs(random_local_noise[i]))
                    random_other_noise[i] = -max(0.1, math.fabs(random_other_noise[i]))
                else:
                    random_local_noise[i] = max(0.01, math.fabs(random_local_noise[i]))
                    random_other_noise[i] = max(0.1, math.fabs(random_other_noise[i]))
            print('Random other noise = ', random_other_noise)

            # 3) DEFINITION OF THE NOISY CONSTRAINTS
            # for i in range(n_players):
            #     if average_constraint_list[i] <= 0:
            #         random_noise[i] = 1 * math.fabs(random_noise[i])         
            new_average_constraint_list = copy.copy([agent.average_constraint for agent in list_agents])   
            # new_average_constraint_list = [min(.9995, max(0, list_agents[i].average_constraint + random_noise[i])) for i in range(n_players)]
            new_discounted_constraint_list = [new_average_constraint_list[i] * (1 - list_agents[i].gamma**(list_agents[i].max_episodes_steps+1))/(1 - list_agents[i].gamma) for i in range(n_players)]

            # 4) COMPUTING THE NOISY OPTIMAL REWARDS
            # 4.1) COMPUTATION OF THE OPTIMAL LOCAL REWARD FOR EACH AGENT, GIVEN THE NOISY LOCAL CONSTRAINT THETA (local constraint: theta_i + noise_i, other constraints: theta_j, j != i)
            # discounted_reward_noisy_local_value, _, _ , average_reward_noisy_local_value, _ = single_agent_policy_learning(list_agents = list_agents, algorithm= 'RL', noisy_average_constraint_list=new_average_constraint_list, noisy_discounted_constraint_list=new_discounted_constraint_list, noise_position = 'local')
            discounted_reward_noisy_local_value, _, _ , average_reward_noisy_local_value, average_cost_local_noise = single_agent_policy_learning(list_agents = list_agents, algorithm= 'RL', local_noise=random_local_noise, other_noise=0, noise_position = 'local')
            list_average_offloading_action.append(average_cost_local_noise)
            
            # 4.2) COMPUTATION OF THE OPTIMAL LOCAL REWARD FOR EACH AGENT, GIVEN THE NOISY GLOBAL CONSTRAINT THETA (local constraint: theta_i, other constraints: theta_j + noise_j, j != i)
            # discounted_reward_noisy_other_values, _, _, average_reward_noisy_other_values, _ = single_agent_policy_learning(list_agents= list_agents, algorithm= 'RL', noisy_average_constraint_list=new_average_constraint_list, noisy_discounted_constraint_list=new_discounted_constraint_list, noise_position = 'other')
            discounted_reward_noisy_other_values, _, _, average_reward_noisy_other_values, average_cost_other_noise = single_agent_policy_learning(list_agents= list_agents, algorithm= 'RL', local_noise=0, other_noise=random_other_noise, noise_position = 'other')
            list_average_offloading_action.append(average_cost_other_noise)


            # 5) COMPUTATION OF THE GRADIENT
            derivative_local_values = np.zeros((n_players, ))
            derivative_other_values = np.zeros((n_players, ))
            for i in range(n_players):
                # derivative_local_values[i] = (discounted_reward_noisy_local_value[i] - discounted_reward_no_noise[i]) / random_local_noise[i]
                if average_cost_local_noise[i] == average_cost_no_noise[i]:
                    derivative_local_values[i] = random_local_noise[i]
                else:
                    derivative_local_values[i] = (discounted_reward_noisy_local_value[i] - discounted_reward_no_noise[i]) / (average_cost_local_noise[i] - average_cost_no_noise[i])

                # total_noise_other_values = np.sum([math.fabs(random_noise[k]) for k in range(n_players) if k != i])
                # total_noise_other_values = np.sum([random_noise[k] for k in range(n_players) if k != i])
                # derivative_other_values[i] = (discounted_reward_noisy_other_values[i] - discounted_reward_no_noise[i]) / random_other_noise[i]
                if average_cost_other_noise[i] == average_cost_no_noise[i]:
                    derivative_other_values[i] = random_other_noise[i]
                else:
                    derivative_other_values[i] = (discounted_reward_noisy_other_values[i] - discounted_reward_no_noise[i]) / (average_cost_other_noise[i] - average_cost_no_noise[i])
            print(random_local_noise, random_other_noise)
            print('Derivative local values = ', derivative_local_values)
            # these values have to be negative
            # derivative_local_values  = np.clip(derivative_local_values, 0, math.inf)
            print('Derivative other values = ', derivative_other_values)
            # these have to be positive
            # derivative_other_values = np.clip(derivative_other_values, -math.inf, 0)
            # we clip the values to avoid too large gradients
            derivative_local_values = np.clip(derivative_local_values, -50, 50)
            derivative_other_values = np.clip(derivative_other_values, -50/n_players, 50/n_players)
            gradient = np.zeros((n_players, ))
            for i in range(n_players):
                gradient[i] = derivative_local_values[i] + np.sum([derivative_other_values[j] for j in range(n_players) if j != i])
            # print("Gradient = {}".format(gradient))
            
            # first we ahve to increment (if necessary) the amount of steps for the optimization
            for i in range(n_players):
                if ( (gradient[i] > 0) + (positive_gradient[i]) ) == 1:
                    constraint_optimization_steps[i] += 1

            # gradient normalization
            gradient = gradient / np.linalg.norm(gradient)
            # print('Normalized gradient = ', gradient)   

            # NOTE: our goal is to minimize the global objective function, so we have to move the constraints in the opposite direction of the gradient
            average_constraint_list = [min(1, max(0, list_agents[i].average_constraint - (constraint_optimization_learning_rate/constraint_optimization_steps[i]**.51) * gradient[i])) for i in range(n_players)]
            for i in range(n_players):
                if gradient[i] < 0:
                    constraint_optimization_steps[i] += 1

            # average_constraint_list = average_constraint_list / np.sum(average_constraint_list)
            # average_constraint_list = average_constraint_list / np.sum(average_constraint_list)
            discounted_constraint_list = [average_constraint_list[i] * (1 - list_agents[i].gamma**(list_agents[i].max_episodes_steps+1))/(1 - list_agents[i].gamma) for i in range(n_players)]
            # print('New constraints: ', average_constraint_list, discounted_constraint_list)

            # now we update the constraint of each agent
            for agent in list_agents:
                agent.average_constraint = average_constraint_list[list_agents.index(agent)]
                agent.discounted_constraint = discounted_constraint_list[list_agents.index(agent)]
                agent.sum_average_constraints = np.sum(average_constraint_list)


            # 5) we check if we have to continue the optimization process (evaluation through simulation)
            #  an idea could be that we stop after K ocnsecutive failures to improve the global objective function

        list_objective_function.append(global_objective_function)
        list_sum_local_rewards.append(np.sum(discounted_reward_no_noise))

        # constraint_optimization_steps += 1

        new_optimization_steps += 1

        if new_optimization_steps > max_iterations:
            constraints_optimized = True
        print('RL----------------------------------------------', new_optimization_steps)

        
    return list_sum_local_rewards, list_objective_function, average_constraint_list, list_average_offloading_action 


def single_experiment_constraint_optimization_with_linear_program(list_agents, constraint_optimization_learning_rate = .5, max_iterations = 10):

    # initial constraint: [0, ..., 0]
    # average_constraint_list = np.zeros((n_players, ))

    # discounted_constraint_list = [average_constraint_list[list_agents.index(agent)] * (1 - agent.gamma**(agent.max_episodes_steps+1))/(1 - agent.gamma) for agent in list_agents]

    constraints_optimized = False
    constraint_optimization_steps = np.ones((n_players, )) 
    positive_gradient = np.zeros((n_players, ))

    new_optimization_steps = 0

    list_objective_function = []
    list_sum_local_rewards = []

    while not constraints_optimized:

        # first we have to compute the sum of the average constraint
        # sum_average_constraints = 0
        # for agent in list_agents:
        #     agent.average_constraint = average_constraint_list[list_agents.index(agent)]
        #     sum_average_constraints += agent.average_constraint 
        # for agent in list_agents:
        #     agent.sum_average_constraints = sum_average_constraints
        # optimized_agent_reward, policy_list, optimized_agent_cost = single_agent_policy_learning(list_agents = list_agents, algorithm='LP', noise_position = 'none')
        
        policy_list = compute_optimal_solution(list_agents)
        
        discounted_reward, discounted_approximated_reward, discounted_cost, average_reward, average_approximated_reward, average_cost = stochastic_policy_evaluation(list_agents, policy_list)
        
        # once we have the optimal local reward for each agent given a certain value of the constraints, it remains to optimize the value of the constraints
        # we still perform the operations in any of the stochastic approximation methods to optimize the global objective function
        # print(optimized_agent_reward)
        # print('Optimized agent cost = ',optimized_agent_cost)
        global_objective_function = global_evaluation(list_agents, policy_list, algorithm = 'LP')
        optimized_average_cost = np.zeros((n_players, ))
        for i in range(n_players):
            optimized_average_cost[i] = optimized_agent_cost[i] / ((1 - list_agents[i].gamma**list_agents[i].max_episodes_steps)/(1 - list_agents[i].gamma))
        # print('Expected mean interaction action: ', optimized_average_cost, np.sum(optimized_average_cost), sum_average_constraints)
        # print('Global objective function = {}. Sum of local rewards = {}'.format(global_objective_function, np.sum(optimized_agent_reward)))

        # we could have a stochastic evaluation of the environment (that considers the real reward obtained by each agent, not the "average" one)

        # 1) we have to define a random noise foe each agent
        random_local_noise = np.random.normal(0, 0.05/constraint_optimization_steps, (n_players, ))
        random_other_noise = np.random.normal(0, 0.05/constraint_optimization_steps, (n_players, ))

        # 2) we have to find the new optimal policy for each agent
        for i in range(n_players):
            if average_constraint_list[i] <= 0:
                random_noise[i] = 1 * math.fabs(random_noise[i])
        # print("Random noise = ", random_noise)
            
        new_average_constraint_list = [min(.9995, max(0, average_constraint_list[i] + random_noise[i])) for i in range(n_players)]
        new_discounted_constraint_list = [new_average_constraint_list[i] * (1 - agent.gamma**(agent.max_episodes_steps+1))/(1 - agent.gamma) for i in range(n_players)]

        # print(new_average_constraint_list)
        # print(new_discounted_constraint_list)

        optimized_agent_reward_noisy_local_value, _, _, average_reward_noisy_local_value, _ = single_agent_policy_learning(list_agents = list_agents, algorithm='LP', local_noise=random_local_noise, other_noise=0, noise_position = 'local')
        optimized_agent_reward_noisy_other_values, _, _, average_reward_noisy_local_value, _ = single_agent_policy_learning(list_agents = list_agents, algorithm='LP', local_noise=o, other_noise=random_other_noise, noise_position = 'other')
        
        # print(optimized_agent_reward_noisy_local_value)
        # print(optimized_agent_reward_noisy_other_values)

        # raise ValueError("Not implemented yet")
        # print("Difference = ", optimized_agent_reward_noisy - optimized_agent_reward)

        gradient_local_values = np.zeros((n_players, ))
        gradient_other_values = np.zeros((n_players, ))
        for i in range(n_players):
            gradient_local_values[i] = (optimized_agent_reward_noisy_local_value[i] - discounted_reward_no_noise[i]) / random_local_noise[i]
            # total_noise_other_values = np.sum([math.fabs(random_noise[i]) for i in range(n_players) if i != list_agents.index(agent)])
            # gradient_other_values[i] = (optimized_agent_reward_noisy_other_values[i] - optimized_agent_reward[i]) / total_noise_other_values
            derivative_other_values[i] = (discounted_reward_noisy_other_values[i] - discounted_reward_no_noise[i]) / random_other_noise[i]

        # print('Gradient local value = ', gradient_local_values)
        # print('Gradient other values = ', gradient_other_values)

        derivative_local_values = np.clip(derivative_local_values, -50, 50)
        derivative_other_values = np.clip(derivative_other_values, -10, 10)

        # computation of the actual gradient
        gradient = np.zeros((n_players, ))
        for i in range(n_players):
            gradient[i] = gradient_local_values[i] + np.sum([gradient_other_values[j] for j in range(n_players) if j != i])
        # print("Gradient = {}".format(gradient))
        
        # 3) we update the values of the constraints

        # first we ahve to increment (if necessary) the amount of steps for the optimization
        for i in range(n_players):
            if ( (gradient[i] > 0) + (positive_gradient[i]) ) == 1:
                constraint_optimization_steps[i] += 1

        gradient = gradient / np.linalg.norm(gradient)
        # print('Normalized gradient = ', gradient) 

        # NOTE: our goal is to minimize the global objective function, so we have to move the constraints in the opposite direction of the gradient
        average_constraint_list = [min(1, max(0, average_constraint_list[i] - (constraint_optimization_learning_rate/constraint_optimization_steps[i]**.51) * gradient[i])) for i in range(n_players)]
        for i in range(n_players):
            if gradient[i] < 0:
                constraint_optimization_steps[i] += 1
        discounted_constraint_list = [average_constraint_list[i] * (1 - agent.gamma**(agent.max_episodes_steps+1))/(1 - agent.gamma) for i in range(n_players)]
        # print('New constraints: ', average_constraint_list, discounted_constraint_list)

        for agent in list_agents:
            agent.average_constraint = average_constraint_list[list_agents.index(agent)]
            agent.discounted_constraint = discounted_constraint_list[list_agents.index(agent)]
            agent.sum_average_constraints = np.sum(average_constraint_list)


        # 5) we check if we have to continue the optimization process (evaluation through simulation)
        #  an idea could be that we stop after K ocnsecutive failures to improve the global objective function

        list_objective_function.append(global_objective_function)
        list_sum_local_rewards.append(np.sum(optimized_agent_reward))

        # constraint_optimization_steps += 1

        new_optimization_steps += 1

        if new_optimization_steps > max_iterations:
            constraints_optimized = True
        print('VI----------------------------------------------', new_optimization_steps)


    return list_sum_local_rewards, list_objective_function


def single_experiment_centralized_optimization(list_agents, tot_steps = 100000, episodes_between_evaluation = 100):

    centralized_env = CentralizedSystem(list_agents)
    if A2C_flag:
        model = A2C("MlpPolicy", centralized_env, verbose=0, device = 'cpu', gamma = list_agents[0].gamma)
    elif PPO_flag:
        model = PPO("MlpPolicy", centralized_env, verbose=0, device = 'cpu', gamma = list_agents[0].gamma)
    # model = PPO("MlpPolicy", centralized_env, verbose=0, device = 'cpu', gamma = list_agents[0].gamma) #``, device="cpu", policy_kwargs=dict(net_arch=[32, 32, 32, 32]), batch_size=64, learning_rate=1e-3, n_epochs=50, n_steps=2048, gamma = .95)
    list_objective_function = []
    n_iterations = int(tot_steps / (list_agents[0].max_episodes_steps * episodes_between_evaluation)) + 1
    for k in range(1, n_iterations + 1):
        # print('Training done')
        # we can now evaluate the model
        model.learn(total_timesteps=list_agents[0].max_episodes_steps * episodes_between_evaluation, reset_num_timesteps=True)
        evaluation = global_evaluation_A2C(centralized_env, model, n_evaluation_episodes = 100)
        print("Reward after {} steps: {}".format(k * list_agents[0].max_episodes_steps * episodes_between_evaluation, evaluation))
        if A2C_flag:
            print("A2C----------------------------------------------------------------", k)
        elif PPO_flag:
            print("PPO----------------------------------------------------------------", k)
        list_objective_function.append(evaluation)


    return list_objective_function

def single_experiment_IQL_individual_reward(list_agents, n_episodes = 1000, episodes_between_evaluation = 100, verbose = False):

    if list_agents is None:
        raise ValueError("list_agents cannot be None")
    if n_episodes <= 0:
        raise ValueError("n_episodes must be a positive integer")

    global_reward_evolution = []
    learning_rate = .05 # 0.5
    exploration_rate = .05
    exploration_decay = 0.995

    q_table = np.zeros((len(list_agents), list_agents[0].num_states, list_agents[0].action_space_size))
    state_visits = np.zeros((len(list_agents), list_agents[0].num_states))

    # edit the table to make sure we never choose the forbidden actions
    for i in range(len(list_agents)):
        for x in range(1, list_agents[i].M+1):
            for e in range(1, -list_agents[i].min_energy):
                index = list_agents[i].compute_state_index(x, -e)
                q_table[i, index, 1] = math.inf
                q_table[i, index, 2] = math.inf

        for x in range(1, list_agents[i].M+1):
            e = list_agents[i].B
            index = list_agents[i].compute_state_index(x, e)
            q_table[i, index, 0] = math.inf

        for e in range(1, list_agents[i].B+1):
            x = list_agents[i].M
            index = list_agents[i].compute_state_index(x, e)
            q_table[i, index, 0] = math.inf

    # Initialize the rewards list
    discounted_rewards = []
    discounted_penalized_rewards = []
    average_rewards = []
    discounted_costs = []
    episode_average_cost = np.zeros((n_episodes, len(list_agents)))

    for episode in range(n_episodes):
        # if episode % 20 == 0 and verbose:
        #     state_visits = np.zeros((len(list_agents), list_agents[0].num_states))

        if episode % episodes_between_evaluation == 0 and episode > 0:
            # we evaluate the performance of the current policies

            policy_list = []
            for i in range(len(list_agents)):
                policy = np.zeros((list_agents[i].M, list_agents[i].B+1))
                for x in range(1, list_agents[i].M+1):
                    for e in range(1, list_agents[i].B+1):
                        state_index = list_agents[i].compute_state_index(x, e)
                        # we choose the action with the minimum Q-value
                        best_action = np.argmin(q_table[i, state_index])
                        policy[x-1, e] = best_action
                policy_list.append(policy)

            # print("Policies computed for all agents.")

            # now we evaluate a global reward with the computed policies
            global_reward = global_evaluation(list_agents, policy_list, n_evaluation_episodes = 50)
            if verbose:
                print("IQL with individual reward: global reward after {} episoded: {}. Average use of offloading: {}".format(episode, global_reward, episode_average_cost[episode-1, :]))
            global_reward_evolution.append(global_reward)

        
        state = []
        for i in range(len(list_agents)):
            state_i = list_agents[i].reset()
            state.append(state_i)
            state_visits[i, state_i['index']] += 1

        # we need to compute the index of each state
        done = False
        total_reward = 0

        episode_discounted_reward = 0
        episode_discounted_penalized_reward = 0
        episode_discounted_cost = 0
        episode_average_reward = 0

        exploration_rate *= exploration_decay


        for t in range(list_agents[0].max_episodes_steps):
            global_action = np.zeros((len(list_agents),))
            for i in range(len(list_agents)):
                if random.uniform(0, 1) < exploration_rate/state_visits[i, state[i]['index']]**.75:
                    # depending on the state, the possible actions are limited
                    if state[i]['e']<0:
                        global_action[i] = 0
                    elif state[i]['e'] == list_agents[i].B or state[i]['x'] == list_agents[i].M:
                        global_action[i] = np.random.choice([1, 2])
                    else:
                        global_action[i] = np.random.randint(0, list_agents[i].action_space_size)
                else:
                    if state[i]['e'] < 0:
                        global_action[i] = 0
                    else:
                        global_action[i] = np.argmin(q_table[i, state[i]['index']])

            global_reward, reward_vector = list_agents[0].reward_with_interactions(state, global_action, return_vector_reward = True)
            next_state = []
            for i in range(len(list_agents)):
                next_state_i, _, _, info = list_agents[i].step(state[i], global_action[i], training=True)
                next_state.append(next_state_i)
            # print("Global action: ", global_action, "Global reward: ", global_reward, "Reward vector: ", reward_vector, "Sum of rewards: ", np.sum(reward_vector))
            # Update the Q-value using the Bellman equation

            # update counter offloading action
            episode_average_cost[episode, :] += (global_action == 2)


            for i in range(len(list_agents)):
                if state_visits[i, state[i]['index']] <= 0:
                    # we completely substitute with the new value
                    q_table[i, state[i]['index'], int(global_action[i])] = reward_vector[i] + list_agents[i].gamma * np.min(q_table[i, next_state[i]['index']])
                else:
                    q_table[i, state[i]['index'], int(global_action[i])] += (learning_rate/state_visits[i, state[i]['index']]**.75) * (reward_vector[i] + list_agents[i].gamma * np.min(q_table[i, next_state[i]['index']]) - q_table[i, state[i]['index'], int(global_action[i])])
            state_visits[i, state[i]['index']] += 1   

            state = next_state

            episode_discounted_reward += (list_agents[0].gamma ** t) * global_reward
            

        # rewards.append(total_reward)
        exploration_rate *= exploration_decay
        discounted_rewards.append(episode_discounted_reward)
        discounted_penalized_rewards.append(episode_discounted_penalized_reward)
        discounted_costs.append(episode_discounted_cost)
        average_rewards.append(episode_average_reward / (t + 1))
        episode_average_cost[episode, :] /= (t + 1)

    if verbose:
        print('Mean discounted reward: ', np.mean(discounted_rewards))
    # global_reward = global_evaluation(list_agents, policy_list, n_evaluation_episodes = 100)
    # print("Global reward for all agents: {}".format(global_reward))


    # now we need to compute the policy for each agent
    policy_list = []
    for i in range(len(list_agents)):
        policy = np.zeros((list_agents[i].M, list_agents[i].B+1))
        for x in range(1, list_agents[i].M+1):
            for e in range(1, list_agents[i].B+1):
                state_index = list_agents[i].compute_state_index(x, e)
                # we choose the action with the minimum Q-value
                best_action = np.argmin(q_table[i, state_index])
                policy[x-1, e] = best_action
        policy_list.append(policy)
    if verbose:
        print("Policies computed for all agents.")

    # now we evaluate a global reward with the computed policies
    global_reward = global_evaluation(list_agents, policy_list, n_evaluation_episodes = 100)
    if verbose:
        print("IQL with individual reward: global reward = {}. Average cost: {}".format(global_reward, np.mean(episode_average_cost[episode, :])))
    global_reward_evolution.append(global_reward)

    return 0, global_reward_evolution, episode_average_cost


def single_experiment_IQL_common_reward(list_agents = None, n_episodes = 1000, episodes_between_evaluation = 100, verbose =False):
    if list_agents is None:
        raise ValueError("list_agents cannot be None")
    if n_episodes <= 0:
        raise ValueError("n_episodes must be a positive integer")


    learning_rate = .05
    global_reward_evolution = []
    exploration_rate = .05
    exploration_decay = 0.995

    q_table = np.zeros((len(list_agents), list_agents[0].num_states, list_agents[0].action_space_size))
    state_visits = np.ones((len(list_agents), list_agents[0].num_states))

    # edit the table to make sure we never choose the forbidden actions
    for i in range(len(list_agents)):
        for x in range(1, list_agents[i].M+1):
            for e in range(1, -list_agents[i].min_energy):
                index = list_agents[i].compute_state_index(x, -e)
                q_table[i, index, 1] = math.inf
                q_table[i, index, 2] = math.inf

        for x in range(1, list_agents[i].M+1):
            e = list_agents[i].B
            index = list_agents[i].compute_state_index(x, e)
            q_table[i, index, 0] = math.inf

        for e in range(1, list_agents[i].B+1):
            x = list_agents[i].M
            index = list_agents[i].compute_state_index(x, e)
            q_table[i, index, 0] = math.inf

    # Initialize the rewards list
    discounted_rewards = []
    discounted_penalized_rewards = []
    average_rewards = []
    discounted_costs = []
    average_costs = []

    for episode in range(n_episodes):

        if episode % episodes_between_evaluation == 0 and episode > 0:
            # we evaluate the performance of the current policies

            policy_list = []
            for i in range(len(list_agents)):
                policy = np.zeros((list_agents[i].M, list_agents[i].B+1))
                for x in range(1, list_agents[i].M+1):
                    for e in range(1, list_agents[i].B+1):
                        state_index = list_agents[i].compute_state_index(x, e)
                        # we choose the action with the minimum Q-value
                        best_action = np.argmin(q_table[i, state_index])
                        policy[x-1, e] = best_action
                policy_list.append(policy)

            # print("Policies computed for all agents.")

            # now we evaluate a global reward with the computed policies
            global_reward = global_evaluation(list_agents, policy_list, n_evaluation_episodes = 50)
            if verbose:
                print("IQL with common reward: global reward after {} episoded: {}".format(episode, global_reward))
            global_reward_evolution.append(global_reward)
        
        state = []
        for i in range(len(list_agents)):
            state_i = list_agents[i].reset()
            state.append(state_i)
            state_visits[i, state_i['index']] += 1

        # we need to compute the index of each state
        done = False
        total_reward = 0

        episode_discounted_reward = 0
        episode_discounted_penalized_reward = 0
        episode_discounted_cost = 0
        episode_average_reward = 0
        episode_average_cost = 0

        exploration_rate *= exploration_decay


        for t in range(list_agents[0].max_episodes_steps):
            global_action = np.zeros((len(list_agents),))
            for i in range(len(list_agents)):
                if random.uniform(0, 1) < exploration_rate/state_visits[i, state[i]['index']]**.75:
                    # depending on the state, the possible actions are limited
                    if state[i]['e']<0:
                        global_action[i] = 0
                    elif state[i]['e'] == list_agents[i].B or state[i]['x'] == list_agents[i].M:
                        global_action[i] = np.random.choice([1, 2])
                    else:
                        global_action[i] = np.random.randint(0, list_agents[i].action_space_size)
                else:
                    if state[i]['e'] < 0:
                        global_action[i] = 0
                    else:
                        global_action[i] = np.argmin(q_table[i, state[i]['index']])

            global_reward = list_agents[0].reward_with_interactions(state, global_action)
            next_state = []
            for i in range(len(list_agents)):
                next_state_i, _, _, info = list_agents[i].step(state[i], global_action[i], training=True)
                next_state.append(next_state_i)

            # Update the Q-value using the Bellman equation
            for i in range(len(list_agents)):
                if state_visits[i, state[i]['index']] <= 0:
                    # we completely substitute with the new value
                    q_table[i, state[i]['index'], int(global_action[i])] = global_reward + list_agents[i].gamma * np.min(q_table[i, next_state[i]['index']])
                else:
                    q_table[i, state[i]['index'], int(global_action[i])] += (learning_rate/state_visits[i, state[i]['index']]**.75) * (global_reward + list_agents[i].gamma * np.min(q_table[i, next_state[i]['index']]) - q_table[i, state[i]['index'], int(global_action[i])])
            state_visits[i, state[i]['index']] += 1   

            state = next_state

            episode_discounted_reward += (list_agents[0].gamma ** t) * global_reward
            

        # rewards.append(total_reward)
        exploration_rate *= exploration_decay
        discounted_rewards.append(episode_discounted_reward)
        discounted_penalized_rewards.append(episode_discounted_penalized_reward)
        discounted_costs.append(episode_discounted_cost)
        average_rewards.append(episode_average_reward / (t + 1))
        average_costs.append(episode_average_cost / (t + 1))

    if verbose:
        print('Mean discounted reward: ', np.mean(discounted_rewards))
    policy_list = []
    for i in range(len(list_agents)):
        policy = np.zeros((list_agents[i].M, list_agents[i].B+1))
        for x in range(1, list_agents[i].M+1):
            for e in range(1, list_agents[i].B+1):
                state_index = list_agents[i].compute_state_index(x, e)
                # we choose the action with the minimum Q-value
                best_action = np.argmin(q_table[i, state_index])
                policy[x-1, e] = best_action
        policy_list.append(policy)

    if verbose:
        print("Policies computed for all agents.")

    # now we evaluate a global reward with the computed policies
    global_reward = global_evaluation(list_agents, policy_list, n_evaluation_episodes = 100)
    if verbose:
        print("IQL with common reward: global reward = {}".format(global_reward))

    global_reward_evolution.append(global_reward)
    
    return 0, global_reward_evolution

random.seed(0)
np.random.seed(0)


n_players = 50
n_experiments = 15
congestion_penalty_exponent = 1

experiments_environments = []

zero_initial_constraints = True

list_chosen_values_min_harvesting_rate = [] 
list_chosen_values_max_harvesting_rate = []
list_chosen_values_min_processing_cost = []
list_chosen_values_max_processing_cost = []


# GENERATION OF THE ENVIRONMENTS
for experiment in range(n_experiments):
    chosen_values_min_harvesting_rate = []
    chosen_values_max_harvesting_rate = []
    chosen_values_min_processing_cost = []
    chosen_values_max_processing_cost = []

    list_agents = []
    total_objective_function = 0
    constraints_optimized = False
    # definition of the distributions for the parameters of each agent
    M_values = [15]
    B_values = [15] 

    min_value_harvesting_list = [0, 1]
    max_value_harvesting_list = [1, 2, 3] #, 4, 5]

    min_value_cost_list = [1, 2]
    max_value_cost_list = [5, 7, 10]


    
    
    
    
    # if congestion_penalty_exponent != 1 and n_players != 10:
    #     raise ValueError("Error in the exponent and number of players")
    # definition of the players in the system
    for player in range(n_players):
        print('Creating agent {}'.format(player+1))
        M = random.choice(M_values)
        B = random.choice(B_values)

        agent = SingleDevice(M = M, B=B, total_players=n_players, congestion_penalty_multiplier = 1, congestion_penalty_exponent = congestion_penalty_exponent) 
        agent.gamma = .95

        harvesting_distribution = np.zeros((agent.B + 1, ))
        min_value_harvesting = random.choice(min_value_harvesting_list)
        chosen_values_min_harvesting_rate.append(min_value_harvesting)
        max_value_harvesting = random.choice(max_value_harvesting_list)
        chosen_values_max_harvesting_rate.append(max_value_harvesting)
        
        for i in range(min_value_harvesting, max_value_harvesting+1):
            harvesting_distribution[i] = 1
        harvesting_distribution = harvesting_distribution / np.sum(harvesting_distribution)
        agent.p_H = harvesting_distribution

        min_value_cost = random.choice(min_value_cost_list)
        chosen_values_min_processing_cost.append(min_value_cost)
        max_value_cost = random.choice(max_value_cost_list)
        chosen_values_max_processing_cost.append(max_value_cost)
        p_C=np.zeros((agent.B+1, ))              # vector of energy cost probabilities (0,1,...,B+h)
        for i in range(min_value_cost, max_value_cost+1):
            p_C[i] = 1
        p_C = p_C / np.sum(p_C)
        agent.p_C = p_C

        print(min_value_harvesting, max_value_harvesting, min_value_cost, max_value_cost)

        agent.average_constraint = 0 # 1/n_players
        agent.discounted_constraint = 0 # agent.average_constraint * (1 - agent.gamma**agent.max_episodes_steps)/(1 - agent.gamma)

        list_agents.append(agent)

    # now we define the constraints for the agents
    if zero_initial_constraints:
        # average_constraint_list = np.random.uniform(0, 1, (n_players, ))
        average_constraint_list = np.zeros((n_players, )) # all the agents can send as much as they want
        # average_constraint_list = 1 * average_constraint_list / np.sum(average_constraint_list)
        # average_constraint_list = .25 * np.ones((n_players, )) # all the agents can send as much as they want
        discounted_constraint_list = [average_constraint_list[i] * (1 - list_agents[i].gamma**(list_agents[i].max_episodes_steps+1))/(1 - list_agents[i].gamma) for i in range(n_players)]
    else:
        average_constraint_list = 0.25 * np.ones((n_players, ))
        discounted_constraint_list = [average_constraint_list[i] * (1 - list_agents[i].gamma**(list_agents[i].max_episodes_steps+1))/(1 - list_agents[i].gamma) for i in range(n_players)]

    for agent in list_agents:
        agent.average_constraint = average_constraint_list[list_agents.index(agent)]
        agent.discounted_constraint = discounted_constraint_list[list_agents.index(agent)]
        agent.sum_average_constraints = np.sum(average_constraint_list)

    experiments_environments.append((list_agents, ))


    list_chosen_values_min_harvesting_rate.append(chosen_values_min_harvesting_rate)
    list_chosen_values_max_harvesting_rate.append(chosen_values_max_harvesting_rate)
    list_chosen_values_min_processing_cost.append(chosen_values_min_processing_cost)
    list_chosen_values_max_processing_cost.append(chosen_values_max_processing_cost)

print(list_chosen_values_min_harvesting_rate)
print(list_chosen_values_max_harvesting_rate)
print(list_chosen_values_min_processing_cost)
print(list_chosen_values_max_processing_cost)

# raise ValueError("Stop here")

cAlg_flag = False
IQL_flag = True
only_individual_IQL = False
A2C_flag = False
PPO_flag = False
# list_inputs = [experiments_environments, cAlg, IQL, A2C]

pool = mp.Pool(15)
results = pool.starmap(single_experiment_wrapper, experiments_environments)
print(results)

print('obtained results')



if True and n_experiments > 1:
    # save results in a file

    # get current folder
    folder = os.getcwd() + '/results/'

    # create directory with the amount of agents
    if congestion_penalty_exponent == 1:
        folder = folder + str(n_players) + '_agents' + '/'
    elif congestion_penalty_exponent > 0:
        folder = folder + str(n_players) + '_agents_' + str(congestion_penalty_exponent) +  '_penalty_exponent/'

    if not os.path.exists(folder):
        os.makedirs(folder)

    if cAlg_flag:
        if zero_initial_constraints:
            with open(folder + 'results_RL_0_initial_constraint.pkl', 'wb') as f:
                pk.dump(results, f)
        else:
            with open(folder + 'results_RL_uniform_initial.pkl', 'wb') as f:
                pk.dump(results, f)
    elif IQL_flag and not only_individual_IQL:
        with open(folder + 'results_IQL_low_LR.pkl', 'wb') as f:
            pk.dump(results, f)
    elif IQL_flag and only_individual_IQL:
        with open(folder + 'results_IQL_individual_low_LR.pkl', 'wb') as f:
            pk.dump(results, f)
    elif A2C_flag:
        with open(folder + 'results_A2C.pkl', 'wb') as f:
            pk.dump(results, f)
    elif PPO_flag:
        with open(folder + 'results_PPO.pkl', 'wb') as f:
            pk.dump(results, f)
    else:
        raise ValueError("No flag found")

    print('saved results')

    plt.plot()
    plt.show()