geneticAlgorithm.py

import random
import numpy as np
import matplotlib.pyplot as plt

n = 24  # Number of firms
alpha = 1023  # Demand intercept
c = 0  # Marginal cost
l = 10  # Length of binary string (chromosomes)
T = 1000  # Max iterations
Y = 1  # Second to last digit of ID
Z = 3  # Last digit of ID

p_cross_bench = 0.200 + (0.01 * Y) + (0.001 * Z)
p_mut_bench = 0.00200 + (0.0001 * Y) + (0.00001 * Z)

# Theoretical Values
Q_nash_total = (n / (n + 1)) * (alpha - c)
q_nash_indiv = Q_nash_total / n


def BinaryToDecimal(binary_array):
    value = 0
    for bit in binary_array:
        value = (value * 2) + bit
    return value


def BinaryToOutput(binary_array):
    x_it = BinaryToDecimal(binary_array)
    q_i = alpha * (
        x_it / (2**l)
    )  # Scale the integer [0, 1023] to output (Note: 2^10 = 1024. Max x_it is 1023.)
    return q_i


def CalculatePrice(total_Q):
    price = alpha - total_Q
    return price


def CalculateProfit(q_i, price):
    return (price * q_i) - (c * q_i)


def Variance(values):
    if len(values) == 0:
        return 0
    mean = sum(values) / len(values)
    variance = sum((x - mean) ** 2 for x in values) / len(values)
    return variance


def run_simulation(p_cross, p_mut, scenario_name):
    print(f"\n--- Running Scenario: {scenario_name} ---")
    print(f"p_cross: {p_cross:.5f}, p_mut: {p_mut:.5f}")

    random.seed(42)
    np.random.seed(42)

    Population = []
    for i in range(n):
        # Generate random binary string of length l
        genes = [random.randint(0, 1) for _ in range(l)]
        Population.append(genes)

    # Storage for statistics
    total_Q_history = []
    variance_history = []

    for t in range(1, T + 1):
        outputs_list = []
        total_Q = 0

        # Calculate outputs for current population
        for firm in Population:
            q_i = BinaryToOutput(firm)
            outputs_list.append(q_i)
            total_Q += q_i

        current_price = CalculatePrice(total_Q)

        # Store stats
        total_Q_history.append(total_Q)
        variance_history.append(Variance(outputs_list))

        # A. RECOMBINATION (Roulette Wheel Selection)
        fitness_scores = []
        for firm in Population:
            q_i = BinaryToOutput(firm)
            profit = CalculateProfit(q_i, current_price)
            fitness_scores.append(profit)

        # Handle negative fitness for probability calculation
        min_fitness = min(fitness_scores)
        if min_fitness < 0:
            # Shift all fitness values to be positive
            offset = -min_fitness
            # Add a small epsilon to ensure non-zero probabilities if all are equal
            fitness_scores = [f + offset + 1e-6 for f in fitness_scores]
            # But this solution can left useless in a situation where all firms have zero fitness score

        total_fitness = sum(fitness_scores)

        if total_fitness == 0:
            probabilities = [
                1.0 / n
            ] * n  # If all firms have zero fitness score, each firm has equal probability of being selected
        else:
            probabilities = [f / total_fitness for f in fitness_scores]

        # Create mating pool
        mating_pool = []
        indices = np.random.choice(len(Population), size=n, p=probabilities)
        for idx in indices:
            mating_pool.append(Population[idx].copy())

        # B. CROSSOVER
        offspring_pool = []
        for i in range(0, n, 2):
            Parent1 = mating_pool[i].copy()
            Parent2 = mating_pool[i + 1].copy()

            Child1 = Parent1.copy()
            Child2 = Parent2.copy()

            if random.random() < p_cross:
                k = random.randint(1, l - 1)
                Child1 = Parent1[:k] + Parent2[k:]
                Child2 = Parent2[:k] + Parent1[k:]

            offspring_pool.append(Child1)
            offspring_pool.append(Child2)

        # C. MUTATION
        for child in offspring_pool:
            for bit_idx in range(l):
                if random.random() < p_mut:
                    child[bit_idx] = 1 - child[bit_idx]

        # D. ELECTION
        Next_Generation = []

        for j in range(n):
            Parent = Population[j]
            Child = offspring_pool[j]

            q_parent = BinaryToOutput(Parent)
            q_child = BinaryToOutput(Child)

            # Calculate Q_others (aggregate of all others in current generation)
            Q_others = total_Q - q_parent

            # Parent Profit (Actual)
            profit_parent = CalculateProfit(q_parent, current_price)

            # Child Profit (Hypothetical)
            # Evaluate Child against the *current* behavior of others
            new_total_Q = Q_others + q_child
            new_price = CalculatePrice(new_total_Q)
            profit_child = CalculateProfit(q_child, new_price)

            if profit_child > profit_parent:
                Next_Generation.append(Child)
            else:
                Next_Generation.append(Parent)

        Population = Next_Generation

    # RESULTS & PLOTTING
    print(f"Simulation completed!")
    print(f"Final Total Q: {total_Q_history[-1]:.2f}")
    print(f"Target Nash Total Q: {Q_nash_total:.2f}")
    print(f"Final Variance: {variance_history[-1]:.2f}")

    # Plotting
    plt.figure(figsize=(12, 5))

    # Plot 1: Total Quantity
    plt.subplot(1, 2, 1)
    plt.plot(total_Q_history, label="Actual Total Q", color="blue")
    plt.axhline(y=Q_nash_total, color="r", linestyle="--", label="Nash Equilibrium")
    plt.xlabel("Iteration")
    plt.ylabel("Total Quantity (Q)")
    plt.title(f"Industry Output Convergence\n({scenario_name})")
    plt.legend()
    plt.grid(True, alpha=0.3)

    # Plot 2: Variance
    plt.subplot(1, 2, 2)
    plt.plot(variance_history, color="green")
    plt.xlabel("Iteration")
    plt.ylabel("Variance")
    plt.title(f"Population Variance\n({scenario_name})")
    plt.grid(True, alpha=0.3)

    plt.tight_layout()
    filename = f"results_{scenario_name.replace(' ', '_')}.png"
    plt.savefig(f"results/{filename}")
    print(f"Saved plot to results/{filename}")
    plt.close()


if __name__ == "__main__":
    # (ii) Benchmark Scenario
    run_simulation(p_cross_bench, p_mut_bench, "Benchmark")

    # (iii) High Crossover (+0.250)
    run_simulation(p_cross_bench + 0.250, p_mut_bench, "High_Crossover")

    # (iv) High Mutation (+0.00250)
    run_simulation(p_cross_bench, p_mut_bench + 0.00250, "High_Mutation")