# Import necessary libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from pgmpy.models import DiscreteBayesianNetwork  # Updated import
from pgmpy.factors.discrete import TabularCPD
from pgmpy.inference import VariableElimination
import networkx as nx
from matplotlib.colors import LinearSegmentedColormap

# Set style for plots
plt.style.use('ggplot')
sns.set_palette("muted")
sns.set_context("notebook", font_scale=1.2)

# Define the Bayesian Network structure
# Edges represent dependencies between variables
model = DiscreteBayesianNetwork([  # Updated class
    ('SupplierQuality', 'DefectRate'),
    ('MachineCalibration', 'DefectRate'),
    ('WorkerExperience', 'DefectRate'),
    ('DefectRate', 'QCDecision'),
    ('QCDecision', 'Cost')
])

# Define the conditional probability distributions (CPDs)

# 1. Supplier Quality (Prior)
cpd_supplier = TabularCPD(
    variable='SupplierQuality',
    variable_card=2,
    values=[[0.7], [0.3]],  # [High, Low]
    state_names={'SupplierQuality': ['High', 'Low']}
)

# 2. Machine Calibration (Prior)
cpd_machine = TabularCPD(
    variable='MachineCalibration',
    variable_card=2,
    values=[[0.8], [0.2]],  # [Good, Poor]
    state_names={'MachineCalibration': ['Good', 'Poor']}
)

# 3. Worker Experience (Prior)
cpd_worker = TabularCPD(
    variable='WorkerExperience',
    variable_card=2,
    values=[[0.6], [0.4]],  # [Experienced, Inexperienced]
    state_names={'WorkerExperience': ['Experienced', 'Inexperienced']}
)

# 4. Defect Rate (Conditional on Supplier, Machine, and Worker)
# This is a more complex CPD with multiple parents
cpd_defect = TabularCPD(
    variable='DefectRate',
    variable_card=2,
    values=[
        # Low defect rate probabilities
        [0.95, 0.80, 0.70, 0.50, 0.65, 0.40, 0.30, 0.10],
        # High defect rate probabilities
        [0.05, 0.20, 0.30, 0.50, 0.35, 0.60, 0.70, 0.90]
    ],
    evidence=['SupplierQuality', 'MachineCalibration', 'WorkerExperience'],
    evidence_card=[2, 2, 2],
    state_names={
        'DefectRate': ['Low', 'High'],
        'SupplierQuality': ['High', 'Low'],
        'MachineCalibration': ['Good', 'Poor'],
        'WorkerExperience': ['Experienced', 'Inexperienced']
    }
)

# 5. Quality Control Decision (Conditional on Defect Rate)
cpd_qc = TabularCPD(
    variable='QCDecision',
    variable_card=2,
    values=[
        [0.7, 0.2],  # Basic inspection probabilities
        [0.3, 0.8]   # Detailed inspection probabilities
    ],
    evidence=['DefectRate'],
    evidence_card=[2],
    state_names={
        'QCDecision': ['BasicInspection', 'DetailedInspection'],
        'DefectRate': ['Low', 'High']
    }
)

# 6. Cost (Conditional on QC Decision)
# The values represent expected costs in monetary units
cpd_cost = TabularCPD(
    variable='Cost',
    variable_card=3,
    values=[
        [0.8, 0.1],  # Low cost probabilities
        [0.15, 0.3], # Medium cost probabilities
        [0.05, 0.6]  # High cost probabilities
    ],
    evidence=['QCDecision'],
    evidence_card=[2],
    state_names={
        'Cost': ['Low', 'Medium', 'High'],
        'QCDecision': ['BasicInspection', 'DetailedInspection']
    }
)

# Add CPDs to the model
model.add_cpds(cpd_supplier, cpd_machine, cpd_worker, cpd_defect, cpd_qc, cpd_cost)

# Check if the model is valid
print(f"Is the model valid? {model.check_model()}")

# Create an inference object
inference = VariableElimination(model)

# Function to visualize the Bayesian network structure
def visualize_network(model):
    plt.figure(figsize=(12, 8))
    
    # Create a directed graph
    G = nx.DiGraph()
    
    # Add nodes and edges from the Bayesian network
    for node in model.nodes():
        G.add_node(node)
    
    for edge in model.edges():
        G.add_edge(edge[0], edge[1])
    
    # Define node positions using a hierarchical layout
    pos = {
        'SupplierQuality': (0, 2),
        'MachineCalibration': (1, 2),
        'WorkerExperience': (2, 2),
        'DefectRate': (1, 1),
        'QCDecision': (1, 0),
        'Cost': (1, -1)
    }
    
    # Define node colors
    node_colors = {
        'SupplierQuality': 'skyblue',
        'MachineCalibration': 'skyblue',
        'WorkerExperience': 'skyblue',
        'DefectRate': 'lightgreen',
        'QCDecision': 'salmon',
        'Cost': 'gold'
    }
    colors = [node_colors[node] for node in G.nodes()]
    
    # Draw the graph
    nx.draw(G, pos, with_labels=True, node_color=colors, node_size=3000, 
            font_size=12, font_weight='bold', arrowsize=20, alpha=0.8)
    
    plt.title("Bayesian Network for Manufacturing Quality Control", fontsize=16)
    plt.tight_layout()
    return plt

# Visualize the network
network_plot = visualize_network(model)
network_plot.savefig('bayesian_network.png', dpi=300, bbox_inches='tight')
network_plot.show()

# Perform various queries to support decision making

# 1. Query the probability of defect rates
defect_prob = inference.query(variables=['DefectRate'])
print("\nProbability of Defect Rates:")
print(defect_prob)

# 2. Query the probability of defect rates given supplier quality is high
defect_given_high_quality = inference.query(
    variables=['DefectRate'],
    evidence={'SupplierQuality': 'High'}
)
print("\nProbability of Defect Rates given High Supplier Quality:")
print(defect_given_high_quality)

# 3. Query the probability of defect rates given poor machine calibration
defect_given_poor_calibration = inference.query(
    variables=['DefectRate'],
    evidence={'MachineCalibration': 'Poor'}
)
print("\nProbability of Defect Rates given Poor Machine Calibration:")
print(defect_given_poor_calibration)

# 4. Query the probability of quality control decisions
qc_prob = inference.query(variables=['QCDecision'])
print("\nProbability of Quality Control Decisions:")
print(qc_prob)

# 5. Query the probability of costs
cost_prob = inference.query(variables=['Cost'])
print("\nProbability of Costs:")
print(cost_prob)

# 6. Query the probability of costs given high defect rate
cost_given_high_defect = inference.query(
    variables=['Cost'],
    evidence={'DefectRate': 'High'}
)
print("\nProbability of Costs given High Defect Rate:")
print(cost_given_high_defect)

# 7. What if we know the worker is inexperienced, what's the probability of high defect rate?
defect_given_inexperienced = inference.query(
    variables=['DefectRate'],
    evidence={'WorkerExperience': 'Inexperienced'}
)
print("\nProbability of Defect Rates given Inexperienced Worker:")
print(defect_given_inexperienced)

# 8. Complex query: What is the probability distribution of costs given high supplier quality and poor machine calibration?
cost_complex_evidence = inference.query(
    variables=['Cost'],
    evidence={'SupplierQuality': 'High', 'MachineCalibration': 'Poor'}
)
print("\nProbability of Costs given High Supplier Quality and Poor Machine Calibration:")
print(cost_complex_evidence)

# Visualize the probability distributions
# Extract values from query results
defect_values = defect_prob.values
defect_states = defect_prob.state_names['DefectRate']

defect_high_quality_values = defect_given_high_quality.values
defect_poor_calib_values = defect_given_poor_calibration.values
defect_inexperienced_values = defect_given_inexperienced.values

defect_distributions = [
    (defect_states, defect_values, "Overall"),
    (defect_states, defect_high_quality_values, "High Supplier Quality"),
    (defect_states, defect_poor_calib_values, "Poor Machine Calibration"),
    (defect_states, defect_inexperienced_values, "Inexperienced Worker")
]

# Plot defect rate comparisons
plt.figure(figsize=(14, 8))

x = np.arange(len(defect_states))
width = 0.2
multiplier = 0

for attribute, values, label in defect_distributions:
    offset = width * multiplier
    rects = plt.bar(x + offset, values, width, label=label)
    multiplier += 1

plt.ylabel('Probability', fontsize=14)
plt.title('Defect Rate Probability Under Different Conditions', fontsize=16)
plt.xticks(x + width, defect_states, fontsize=12)
plt.legend(loc='upper left', fontsize=12)
plt.grid(axis='y', linestyle='--', alpha=0.7)

plt.tight_layout()
plt.savefig('defect_rate_comparison.png', dpi=300, bbox_inches='tight')
plt.show()

# Function to create a heatmap for sensitivity analysis
def create_sensitivity_heatmap():
    # Define the factors to vary
    supplier_quality = ['High', 'Low']
    machine_calibration = ['Good', 'Poor']
    worker_experience = ['Experienced', 'Inexperienced']
    
    # Create a DataFrame to store results
    results = pd.DataFrame(columns=['SupplierQuality', 'MachineCalibration', 
                                   'WorkerExperience', 'HighDefectProbability'])
    
    # Run queries for all combinations
    for sq in supplier_quality:
        for mc in machine_calibration:
            for we in worker_experience:
                evidence = {
                    'SupplierQuality': sq,
                    'MachineCalibration': mc,
                    'WorkerExperience': we
                }
                query_result = inference.query(variables=['DefectRate'], evidence=evidence)
                high_defect_prob = query_result.values[1]  # Index 1 corresponds to 'High'
                
                # Append to results
                new_row = pd.DataFrame({
                    'SupplierQuality': [sq],
                    'MachineCalibration': [mc],
                    'WorkerExperience': [we],
                    'HighDefectProbability': [high_defect_prob]
                })
                results = pd.concat([results, new_row], ignore_index=True)
    
    # Create a pivot table for the heatmap
    heatmap_data = results.pivot_table(
        index=['SupplierQuality', 'MachineCalibration'],
        columns='WorkerExperience',
        values='HighDefectProbability'
    )
    
    # Create the heatmap
    plt.figure(figsize=(12, 8))
    custom_cmap = LinearSegmentedColormap.from_list('custom_red_green', ['green', 'yellow', 'red'])
    
    ax = sns.heatmap(heatmap_data, annot=True, cmap=custom_cmap, fmt='.2f', 
                     linewidths=.5, cbar_kws={'label': 'Probability of High Defect Rate'})
    
    plt.title('Sensitivity Analysis: Probability of High Defect Rate', fontsize=16)
    plt.tight_layout()
    
    return plt, results

# Run sensitivity analysis and display heatmap
sensitivity_plot, sensitivity_data = create_sensitivity_heatmap()
sensitivity_plot.savefig('sensitivity_analysis.png', dpi=300, bbox_inches='tight')
sensitivity_plot.show()

# Calculate expected costs for decision analysis
def calculate_expected_costs():
    # Define utility values for each cost level
    cost_utilities = {
        'Low': 100,      # $100 for low cost
        'Medium': 500,   # $500 for medium cost
        'High': 2000     # $2000 for high cost
    }
    
    # Define different evidence scenarios
    scenarios = [
        ('No Evidence', {}),
        ('High Supplier Quality', {'SupplierQuality': 'High'}),
        ('Poor Machine Calibration', {'MachineCalibration': 'Poor'}),
        ('Inexperienced Worker', {'WorkerExperience': 'Inexperienced'}),
        ('High Defect Rate', {'DefectRate': 'High'}),
        ('Low Defect Rate', {'DefectRate': 'Low'})
    ]
    
    # Calculate expected costs for each scenario and decision
    results = []
    
    for scenario_name, evidence in scenarios:
        # Calculate probability distribution for costs given the evidence
        query_result = inference.query(variables=['Cost'], evidence=evidence)
        
        # Calculate expected cost
        expected_cost = 0
        for i, cost_level in enumerate(query_result.state_names['Cost']):
            expected_cost += query_result.values[i] * cost_utilities[cost_level]
        
        results.append((scenario_name, expected_cost))
    
    # Convert to DataFrame for plotting
    df = pd.DataFrame(results, columns=['Scenario', 'Expected Cost'])
    
    # Create the plot
    plt.figure(figsize=(14, 8))
    bars = plt.bar(df['Scenario'], df['Expected Cost'], color='skyblue', alpha=0.7)
    
    # Add data labels
    for bar in bars:
        height = bar.get_height()
        plt.text(bar.get_x() + bar.get_width()/2., height + 20,
                f'${height:.2f}', ha='center', va='bottom', fontsize=12)
    
    plt.xlabel('Scenario', fontsize=14)
    plt.ylabel('Expected Cost ($)', fontsize=14)
    plt.title('Expected Costs Under Different Scenarios', fontsize=16)
    plt.xticks(rotation=45, ha='right', fontsize=12)
    plt.grid(axis='y', linestyle='--', alpha=0.7)
    
    plt.tight_layout()
    return plt, df

# Calculate and display expected costs
cost_analysis_plot, cost_data = calculate_expected_costs()
cost_analysis_plot.savefig('expected_costs.png', dpi=300, bbox_inches='tight')
cost_analysis_plot.show()

# Decision analysis: Should we invest in better machine calibration?

# Current scenario (baseline)
baseline_query = inference.query(variables=['Cost'])
cost_utilities = {'Low': 100, 'Medium': 500, 'High': 2000}
baseline_expected_cost = sum(baseline_query.values[i] * cost_utilities[cost] 
                           for i, cost in enumerate(baseline_query.state_names['Cost']))

# Scenario with perfect machine calibration
perfect_calibration_query = inference.query(
    variables=['Cost'],
    evidence={'MachineCalibration': 'Good'}
)
perfect_calibration_cost = sum(perfect_calibration_query.values[i] * cost_utilities[cost] 
                             for i, cost in enumerate(perfect_calibration_query.state_names['Cost']))

# Cost savings from perfect calibration
cost_savings = baseline_expected_cost - perfect_calibration_cost

# Calculate the maximum amount worth spending on calibration
max_investment = cost_savings

print(f"\nBaseline expected cost: ${baseline_expected_cost:.2f}")
print(f"Expected cost with perfect calibration: ${perfect_calibration_cost:.2f}")
print(f"Expected cost savings: ${cost_savings:.2f}")
print(f"Maximum worth investing in perfect calibration: ${max_investment:.2f}")

# Value of Information Analysis
# How much would it be worth to know the defect rate with certainty?

# Expected cost with no information (baseline)
no_info_cost = baseline_expected_cost

# Expected cost with perfect information about defect rate
# We need to calculate the expected cost for each possible defect rate value,
# weighted by the probability of that value

# Get probabilities of defect rates
defect_probs = inference.query(variables=['DefectRate'])

# Calculate expected cost for each defect rate value
expected_costs_with_info = 0
for i, defect_state in enumerate(defect_probs.state_names['DefectRate']):
    # Query cost given this defect rate
    defect_specific_query = inference.query(
        variables=['Cost'],
        evidence={'DefectRate': defect_state}
    )
    
    # Calculate expected cost for this defect rate
    defect_specific_cost = sum(defect_specific_query.values[j] * cost_utilities[cost] 
                             for j, cost in enumerate(defect_specific_query.state_names['Cost']))
    
    # Weight by probability of this defect rate
    expected_costs_with_info += defect_probs.values[i] * defect_specific_cost

# Value of perfect information
vopi = no_info_cost - expected_costs_with_info

print(f"\nExpected cost with no information: ${no_info_cost:.2f}")
print(f"Expected cost with perfect information about defect rate: ${expected_costs_with_info:.2f}")
print(f"Value of Perfect Information (VOPI): ${vopi:.2f}")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import minimize
import networkx as nx
from matplotlib.colors import LinearSegmentedColormap
import matplotlib.ticker as ticker

# Set random seed for reproducibility
np.random.seed(42)

# Define the supply chain network
num_suppliers = 3
num_factories = 2
num_distribution_centers = 3

# Define the parameters
# Cost and carbon emissions for transportation from suppliers to factories
supplier_factory_cost = np.array([
    [100, 120],  # Supplier 1 to Factory 1 and 2
    [90, 130],   # Supplier 2 to Factory 1 and 2
    [110, 100]   # Supplier 3 to Factory 1 and 2
])

supplier_factory_carbon = np.array([
    [50, 60],    # Supplier 1 to Factory 1 and 2
    [70, 40],    # Supplier 2 to Factory 1 and 2
    [30, 80]     # Supplier 3 to Factory 1 and 2
])

# Cost and carbon emissions for transportation from factories to distribution centers
factory_dc_cost = np.array([
    [80, 90, 110],  # Factory 1 to DC 1, 2, and 3
    [100, 70, 90]   # Factory 2 to DC 1, 2, and 3
])

factory_dc_carbon = np.array([
    [40, 50, 30],   # Factory 1 to DC 1, 2, and 3
    [60, 45, 55]    # Factory 2 to DC 1, 2, and 3
])

# Production capacity at each factory
factory_capacity = np.array([1000, 800])

# Maximum supply capacity from each supplier
supplier_capacity = np.array([600, 700, 800])

# Demand at each distribution center
dc_demand = np.array([400, 500, 600])

# Define the decision variables and constraints for the optimization model
def supply_chain_cost(x):
    """Calculate the total cost of the supply chain."""
    x = x.reshape((num_suppliers, num_factories, num_distribution_centers))
    
    # Calculate cost from suppliers to factories
    sf_cost = 0
    for s in range(num_suppliers):
        for f in range(num_factories):
            sf_cost += supplier_factory_cost[s, f] * np.sum(x[s, f, :])
    
    # Calculate cost from factories to distribution centers
    fd_cost = 0
    for f in range(num_factories):
        for d in range(num_distribution_centers):
            fd_cost += factory_dc_cost[f, d] * np.sum(x[:, f, d])
    
    return sf_cost + fd_cost

def supply_chain_carbon(x):
    """Calculate the total carbon emissions of the supply chain."""
    x = x.reshape((num_suppliers, num_factories, num_distribution_centers))
    
    # Calculate carbon emissions from suppliers to factories
    sf_carbon = 0
    for s in range(num_suppliers):
        for f in range(num_factories):
            sf_carbon += supplier_factory_carbon[s, f] * np.sum(x[s, f, :])
    
    # Calculate carbon emissions from factories to distribution centers
    fd_carbon = 0
    for f in range(num_factories):
        for d in range(num_distribution_centers):
            fd_carbon += factory_dc_carbon[f, d] * np.sum(x[:, f, d])
    
    return sf_carbon + fd_carbon

def objective_function(x, weight_cost=0.5, weight_carbon=0.5):
    """
    Multi-objective function to minimize both cost and carbon emissions.
    The weights determine the importance of each objective.
    """
    # Normalize objectives to similar scales
    cost = supply_chain_cost(x) / 100000
    carbon = supply_chain_carbon(x) / 50000
    
    return weight_cost * cost + weight_carbon * carbon

def constraint_demand(x):
    """Ensure demand at each distribution center is met."""
    x = x.reshape((num_suppliers, num_factories, num_distribution_centers))
    result = []
    
    for d in range(num_distribution_centers):
        # Sum of flow to distribution center d should equal its demand
        flow_to_dc = np.sum(x[:, :, d])
        result.append(flow_to_dc - dc_demand[d])
    
    return np.array(result)

def constraint_factory_capacity(x):
    """Ensure factory capacity constraints are not violated."""
    x = x.reshape((num_suppliers, num_factories, num_distribution_centers))
    result = []
    
    for f in range(num_factories):
        # Sum of flow through factory f should not exceed its capacity
        flow_through_factory = np.sum(x[:, f, :])
        result.append(factory_capacity[f] - flow_through_factory)
    
    return np.array(result)

def constraint_supplier_capacity(x):
    """Ensure supplier capacity constraints are not violated."""
    x = x.reshape((num_suppliers, num_factories, num_distribution_centers))
    result = []
    
    for s in range(num_suppliers):
        # Sum of flow from supplier s should not exceed its capacity
        flow_from_supplier = np.sum(x[s, :, :])
        result.append(supplier_capacity[s] - flow_from_supplier)
    
    return np.array(result)

# Run the optimization for different weight combinations
weights = np.linspace(0, 1, 11)  # 0, 0.1, 0.2, ..., 1.0
results = []

for weight_cost in weights:
    weight_carbon = 1 - weight_cost
    
    # Initial solution: equal distribution
    initial_flow = np.ones((num_suppliers, num_factories, num_distribution_centers)) * 50
    initial_flow = initial_flow.flatten()
    
    # Define constraints
    constraints = [
        {'type': 'eq', 'fun': constraint_demand},
        {'type': 'ineq', 'fun': constraint_factory_capacity},
        {'type': 'ineq', 'fun': constraint_supplier_capacity}
    ]
    
    # Bounds: flow must be non-negative
    bounds = [(0, None) for _ in range(len(initial_flow))]
    
    # Solve the optimization problem
    solution = minimize(
        fun=lambda x: objective_function(x, weight_cost, weight_carbon),
        x0=initial_flow,
        method='SLSQP',
        bounds=bounds,
        constraints=constraints,
        options={'maxiter': 1000}
    )
    
    if solution.success:
        optimal_flow = solution.x.reshape((num_suppliers, num_factories, num_distribution_centers))
        cost = supply_chain_cost(solution.x)
        carbon = supply_chain_carbon(solution.x)
        results.append({
            'weight_cost': weight_cost,
            'weight_carbon': weight_carbon,
            'total_cost': cost,
            'total_carbon': carbon,
            'flow': optimal_flow
        })
        print(f"Weight cost: {weight_cost:.1f}, Weight carbon: {weight_carbon:.1f}, Cost: {cost:.2f}, Carbon: {carbon:.2f}")
    else:
        print(f"Optimization failed for weights: {weight_cost:.1f}, {weight_carbon:.1f}")

# Visualize the Pareto front
pareto_costs = [result['total_cost'] for result in results]
pareto_carbon = [result['total_carbon'] for result in results]

plt.figure(figsize=(10, 6))
plt.scatter(pareto_costs, pareto_carbon, c=weights, cmap='viridis', s=100)
cbar = plt.colorbar()
cbar.set_label('Weight on Cost Objective')
plt.xlabel('Total Cost ($)')
plt.ylabel('Total Carbon Emissions (kg CO₂)')
plt.title('Pareto Front: Trade-off between Cost and Carbon Emissions')
plt.grid(True, alpha=0.3)

# Add annotations for some key points
for i in [0, 5, 10]:
    if i < len(results):
        plt.annotate(
            f"Cost weight: {results[i]['weight_cost']:.1f}",
            (results[i]['total_cost'], results[i]['total_carbon']),
            textcoords="offset points",
            xytext=(0, 10),
            ha='center'
        )

plt.tight_layout()
plt.savefig('pareto_front.png', dpi=300)
plt.show()

# Visualize the optimal flow network for a selected solution (balanced weights)
# Let's use the middle solution (weight_cost = 0.5, weight_carbon = 0.5)
balanced_solution = results[5]  # Index 5 corresponds to weights 0.5, 0.5
optimal_flow = balanced_solution['flow']

# Create a directed graph
G = nx.DiGraph()

# Add nodes
suppliers = [f"S{i+1}" for i in range(num_suppliers)]
factories = [f"F{i+1}" for i in range(num_factories)]
dcs = [f"DC{i+1}" for i in range(num_distribution_centers)]

# Add all nodes
for node in suppliers + factories + dcs:
    G.add_node(node)

# Set node positions
pos = {}
# Suppliers on the left
for i, s in enumerate(suppliers):
    pos[s] = (0, (num_suppliers-1)/2 - i)
# Factories in the middle
for i, f in enumerate(factories):
    pos[f] = (1, (num_factories-1)/2 - i)
# Distribution centers on the right
for i, d in enumerate(dcs):
    pos[d] = (2, (num_distribution_centers-1)/2 - i)

# Add edges with weights based on the optimal flow
edge_widths = []
edge_labels = {}

# Supplier to factory connections
for s in range(num_suppliers):
    for f in range(num_factories):
        flow_amount = np.sum(optimal_flow[s, f, :])
        if flow_amount > 0:
            G.add_edge(suppliers[s], factories[f], weight=flow_amount)
            edge_widths.append(flow_amount / 100)  # Scale down for visualization
            edge_labels[(suppliers[s], factories[f])] = f"{int(flow_amount)}"

# Factory to distribution center connections
for f in range(num_factories):
    for d in range(num_distribution_centers):
        flow_amount = np.sum(optimal_flow[:, f, d])
        if flow_amount > 0:
            G.add_edge(factories[f], dcs[d], weight=flow_amount)
            edge_widths.append(flow_amount / 100)  # Scale down for visualization
            edge_labels[(factories[f], dcs[d])] = f"{int(flow_amount)}"

# Plot the graph
plt.figure(figsize=(12, 8))
node_colors = ['skyblue'] * num_suppliers + ['lightgreen'] * num_factories + ['lightcoral'] * num_distribution_centers
node_sizes = [1500] * (num_suppliers + num_factories + num_distribution_centers)

nx.draw_networkx(
    G, pos, 
    node_color=node_colors,
    node_size=node_sizes,
    with_labels=True,
    font_size=10,
    font_weight='bold',
    arrowsize=15,
    width=edge_widths,
    edge_color='gray',
    alpha=0.7
)

nx.draw_networkx_edge_labels(
    G, pos,
    edge_labels=edge_labels,
    font_size=8
)

plt.title(f"Optimal Supply Chain Flow (Cost Weight: {balanced_solution['weight_cost']}, Carbon Weight: {balanced_solution['weight_carbon']})")
plt.axis('off')
plt.tight_layout()
plt.savefig('optimal_flow_network.png', dpi=300)
plt.show()

# Create a heatmap visualization of the flow
fig, axes = plt.subplots(1, 2, figsize=(15, 6))

# Aggregate flows for Supplier to Factory
sf_flow = np.zeros((num_suppliers, num_factories))
for s in range(num_suppliers):
    for f in range(num_factories):
        sf_flow[s, f] = np.sum(optimal_flow[s, f, :])

# Aggregate flows for Factory to DC
fd_flow = np.zeros((num_factories, num_distribution_centers))
for f in range(num_factories):
    for d in range(num_distribution_centers):
        fd_flow[f, d] = np.sum(optimal_flow[:, f, d])

# Create custom color map (white to green)
green_cmap = LinearSegmentedColormap.from_list('GreenMap', ['white', 'forestgreen'])

# Plot heatmap for Supplier to Factory flows
sns.heatmap(sf_flow, annot=True, fmt='.0f', cmap=green_cmap, ax=axes[0], 
            xticklabels=[f'Factory {i+1}' for i in range(num_factories)],
            yticklabels=[f'Supplier {i+1}' for i in range(num_suppliers)])
axes[0].set_title('Material Flow: Suppliers to Factories')

# Plot heatmap for Factory to DC flows
sns.heatmap(fd_flow, annot=True, fmt='.0f', cmap=green_cmap, ax=axes[1],
            xticklabels=[f'DC {i+1}' for i in range(num_distribution_centers)],
            yticklabels=[f'Factory {i+1}' for i in range(num_factories)])
axes[1].set_title('Product Flow: Factories to Distribution Centers')

plt.tight_layout()
plt.savefig('flow_heatmaps.png', dpi=300)
plt.show()

# Visualize the trade-off details
costs = [result['total_cost'] for result in results]
carbon = [result['total_carbon'] for result in results]
weights_cost = [result['weight_cost'] for result in results]

fig, ax1 = plt.subplots(figsize=(10, 6))

color = 'tab:blue'
ax1.set_xlabel('Cost Weight')
ax1.set_ylabel('Total Cost ($)', color=color)
line1 = ax1.plot(weights_cost, costs, 'o-', color=color, label='Total Cost')
ax1.tick_params(axis='y', labelcolor=color)
ax1.grid(True, alpha=0.3)

ax2 = ax1.twinx()
color = 'tab:green'
ax2.set_ylabel('Total Carbon Emissions (kg CO₂)', color=color)
line2 = ax2.plot(weights_cost, carbon, 's-', color=color, label='Carbon Emissions')
ax2.tick_params(axis='y', labelcolor=color)

# Add combined legend
lines = line1 + line2
labels = [l.get_label() for l in lines]
ax1.legend(lines, labels, loc='upper center')

plt.title('Effect of Cost Weight on Total Cost and Carbon Emissions')
plt.tight_layout()
plt.savefig('weight_effects.png', dpi=300)
plt.show()

# Let's calculate and display the percentage changes
min_cost = min(costs)
max_cost = max(costs)
min_carbon = min(carbon)
max_carbon = max(carbon)

cost_range = max_cost - min_cost
carbon_range = max_carbon - min_carbon

# Create a table of results
results_df = pd.DataFrame({
    'Cost Weight': weights_cost,
    'Carbon Weight': [1-w for w in weights_cost],
    'Total Cost': costs,
    'Total Carbon': carbon,
    'Cost % of Max': [(c - min_cost) / cost_range * 100 for c in costs],
    'Carbon % of Max': [(c - min_carbon) / carbon_range * 100 for c in carbon]
})

# Display the formatted table
pd.set_option('display.float_format', '{:.2f}'.format)
print("\nTrade-off Analysis Table:")
print(results_df)

# Calculate the mathematical model explicitly
# Let's express the mathematical optimization model in LaTeX format
print("\nMathematical Model:")
print("Minimize: $\\alpha \\sum_{s=1}^{S} \\sum_{f=1}^{F} \\sum_{d=1}^{D} (c_{sf} + c_{fd}) x_{sfd} + (1-\\alpha) \\sum_{s=1}^{S} \\sum_{f=1}^{F} \\sum_{d=1}^{D} (e_{sf} + e_{fd}) x_{sfd}$")
print("\nSubject to:")
print("$\\sum_{s=1}^{S} \\sum_{f=1}^{F} x_{sfd} = d_d, \\forall d \\in D$ (Demand Constraints)")
print("$\\sum_{s=1}^{S} \\sum_{d=1}^{D} x_{sfd} \\leq cap_f, \\forall f \\in F$ (Factory Capacity Constraints)")
print("$\\sum_{f=1}^{F} \\sum_{d=1}^{D} x_{sfd} \\leq sup_s, \\forall s \\in S$ (Supplier Capacity Constraints)")
print("$x_{sfd} \\geq 0, \\forall s \\in S, f \\in F, d \\in D$ (Non-negativity Constraints)")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pulp
from matplotlib.ticker import MaxNLocator
import seaborn as sns

# Set seed for reproducibility
np.random.seed(42)

# Define the problem parameters
num_hours = 24
num_generators = 4

# Generator types and their properties
generator_types = ['Coal', 'Natural Gas', 'Solar', 'Wind']
generator_costs = {'Coal': 100, 'Natural Gas': 80, 'Solar': 5, 'Wind': 10}  # Cost per MWh
generator_min_output = {'Coal': 50, 'Natural Gas': 20, 'Solar': 0, 'Wind': 0}  # Min output in MW
generator_max_output = {'Coal': 200, 'Natural Gas': 150, 'Solar': 100, 'Wind': 120}  # Max output in MW

# Set up time-dependent availability for renewables (solar and wind)
solar_availability = np.zeros(num_hours)
wind_availability = np.zeros(num_hours)

# Solar availability follows a bell curve (daylight hours)
for t in range(num_hours):
    if 6 <= t < 18:  # Daylight hours (6 AM to 6 PM)
        solar_availability[t] = 0.8 * np.sin(np.pi * (t - 6) / 12) + 0.1
    else:
        solar_availability[t] = 0.1  # Minimal solar at night
        
# Wind availability has some variability throughout the day
wind_availability = 0.4 + 0.2 * np.sin(np.pi * np.arange(num_hours) / 12) + 0.1 * np.random.rand(num_hours)

# Demand profile throughout the day (in MW)
demand = 200 + 100 * np.sin(np.pi * (np.arange(num_hours) - 4) / 12) + 50 * np.sin(np.pi * (np.arange(num_hours) - 8) / 8)
demand = demand.astype(int)

# Storage parameters
storage_capacity = 150  # Maximum storage capacity in MWh
storage_cost = 15  # Cost per MWh of charging/discharging
initial_storage = 50  # Initial storage level in MWh

# Create the optimization model
model = pulp.LpProblem("Smart_Grid_Optimization", pulp.LpMinimize)

# Decision variables
# Power output for each generator at each hour
power_output = {(g, t): pulp.LpVariable(f"Power_{g}_{t}", 
                                         lowBound=generator_min_output[g], 
                                         upBound=generator_max_output[g] * (1 if g not in ['Solar', 'Wind'] else 
                                                                           (solar_availability[t] if g == 'Solar' else wind_availability[t])), 
                                         cat='Continuous')
                for g in generator_types for t in range(num_hours)}

# Storage level at each hour
storage = {t: pulp.LpVariable(f"Storage_{t}", lowBound=0, upBound=storage_capacity, cat='Continuous') 
          for t in range(num_hours)}

# Variables to track charging and discharging (for cost calculation)
charging = {t: pulp.LpVariable(f"Charging_{t}", lowBound=0, cat='Continuous') 
          for t in range(num_hours)}
discharging = {t: pulp.LpVariable(f"Discharging_{t}", lowBound=0, cat='Continuous') 
          for t in range(num_hours)}

# Objective function: minimize total cost
model += pulp.lpSum([generator_costs[g] * power_output[g, t] for g in generator_types for t in range(num_hours)]) + \
         pulp.lpSum([storage_cost * (charging[t] + discharging[t]) for t in range(num_hours)])

# Constraints
# Initial storage level
model += storage[0] == initial_storage

# Storage dynamics
for t in range(1, num_hours):
    model += storage[t] - storage[t-1] == charging[t] - discharging[t]
    
# Power balance: supply = demand at each hour
for t in range(num_hours):
    model += pulp.lpSum([power_output[g, t] for g in generator_types]) + discharging[t] - charging[t] == demand[t]

# Solve the model
model.solve(pulp.PULP_CBC_CMD(msg=False))

print(f"Status: {pulp.LpStatus[model.status]}")
print(f"Total Cost: ${pulp.value(model.objective):.2f}")

# Extract results
results = pd.DataFrame(index=range(num_hours), columns=generator_types + ['Storage', 'Demand'])

for t in range(num_hours):
    for g in generator_types:
        results.loc[t, g] = power_output[g, t].value()
    results.loc[t, 'Storage'] = storage[t].value()
    results.loc[t, 'Demand'] = demand[t]

# Plot the results
plt.figure(figsize=(15, 10))

# Plot 1: Generator outputs and demand
ax1 = plt.subplot(2, 1, 1)
bottom = np.zeros(num_hours)

for gen in generator_types:
    plt.bar(range(num_hours), results[gen], bottom=bottom, label=gen, alpha=0.7)
    bottom += results[gen].values

plt.plot(range(num_hours), results['Demand'], 'r-', linewidth=2, label='Demand')
plt.title('Hourly Power Generation by Source vs Demand', fontsize=14)
plt.xlabel('Hour of Day')
plt.ylabel('Power (MW)')
plt.legend()
plt.xticks(range(0, num_hours, 2))
plt.grid(True, alpha=0.3)

# Plot 2: Storage level
ax2 = plt.subplot(2, 1, 2)
plt.plot(range(num_hours), results['Storage'], 'g-', linewidth=2, marker='o')
plt.fill_between(range(num_hours), results['Storage'], alpha=0.3, color='green')
plt.title('Energy Storage Level Over Time', fontsize=14)
plt.xlabel('Hour of Day')
plt.ylabel('Storage Level (MWh)')
plt.xticks(range(0, num_hours, 2))
plt.grid(True, alpha=0.3)
plt.tight_layout()

# Additional plot: Cost breakdown
costs = pd.DataFrame(index=generator_types + ['Storage'])
for g in generator_types:
    costs.loc[g, 'Cost'] = sum(generator_costs[g] * power_output[g, t].value() for t in range(num_hours))
    
storage_total_cost = 0
for t in range(num_hours):
    if t > 0:
        storage_total_cost += storage_cost * abs(storage[t].value() - storage[t-1].value())
        
costs.loc['Storage', 'Cost'] = storage_total_cost

plt.figure(figsize=(10, 6))
sns.barplot(x=costs.index, y='Cost', data=costs)
plt.title('Total Cost by Energy Source', fontsize=14)
plt.xlabel('Source')
plt.ylabel('Cost ($)')
plt.grid(True, alpha=0.3)
plt.tight_layout()

# Export results to CSV
results.to_csv('smart_grid_optimization_results.csv')

plt.show()

# Print hourly breakdown
print("\nHourly Power Generation and Storage:")
print(results.round(2))

# Calculate some key metrics
total_renewable = results['Solar'].sum() + results['Wind'].sum()
total_conventional = results['Coal'].sum() + results['Natural Gas'].sum()
total_generation = total_renewable + total_conventional

print(f"\nPercentage of Renewable Energy: {100 * total_renewable / total_generation:.2f}%")
print(f"Maximum Storage Used: {results['Storage'].max():.2f} MWh")
print(f"Storage Utilization: {100 * results['Storage'].max() / storage_capacity:.2f}%")

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import pandas as pd
import seaborn as sns
from scipy.stats import norm

# Parameters
S0 = 90  # Current PV of expected cash flows ($ million)
K = 100  # Investment cost ($ million)
T = 5    # Time to expiration (years)
r = 0.05 # Risk-free rate
sigma = 0.3 # Volatility of the underlying asset

# Number of steps in the binomial tree
n_steps = 5

# Calculate up and down factors
dt = T / n_steps
u = np.exp(sigma * np.sqrt(dt))  # Up factor
d = 1 / u  # Down factor
p = (np.exp(r * dt) - d) / (u - d)  # Risk-neutral probability

# Function to create the underlying asset price tree
def create_asset_tree(S0, u, d, n_steps):
    tree = np.zeros((n_steps + 1, n_steps + 1))
    for i in range(n_steps + 1):
        for j in range(i + 1):
            tree[j, i] = S0 * (u ** (i - j)) * (d ** j)
    return tree

# Function to calculate the option value tree
def calculate_option_tree(asset_tree, K, r, dt, p, n_steps):
    option_tree = np.zeros((n_steps + 1, n_steps + 1))
    
    # Terminal payoffs (at expiration)
    for i in range(n_steps + 1):
        option_tree[i, n_steps] = max(0, asset_tree[i, n_steps] - K)
    
    # Backward induction
    for i in range(n_steps - 1, -1, -1):
        for j in range(i + 1):
            # Expected value of continuation
            continuation_value = np.exp(-r * dt) * (p * option_tree[j, i + 1] + (1 - p) * option_tree[j + 1, i + 1])
            # Immediate exercise value
            exercise_value = max(0, asset_tree[j, i] - K)
            # Option value is maximum of continuation and exercise
            option_tree[j, i] = max(continuation_value, exercise_value)
    
    return option_tree

# Create the asset price tree
asset_tree = create_asset_tree(S0, u, d, n_steps)

# Calculate the option value tree
option_tree = calculate_option_tree(asset_tree, K, r, dt, p, n_steps)

# Print the option value (value of flexibility)
real_option_value = option_tree[0, 0]
print(f"Traditional NPV: ${S0 - K} million")
print(f"Real Option Value: ${real_option_value:.2f} million")
print(f"Value of flexibility: ${real_option_value - (S0 - K):.2f} million")

# Create more detailed visualization
def plot_binomial_tree(tree, title, cmap=cm.viridis):
    n_steps = tree.shape[1] - 1
    fig, ax = plt.subplots(figsize=(12, 8))
    
    # Create x and y coordinates for each node
    x_coords = []
    y_coords = []
    values = []
    
    for i in range(n_steps + 1):
        for j in range(i + 1):
            x_coords.append(i)
            y_coords.append(j)
            values.append(tree[j, i])
    
    # Normalize values for color mapping
    normalized_values = (values - np.min(values)) / (np.max(values) - np.min(values))
    colors = cmap(normalized_values)
    
    # Plot points
    scatter = ax.scatter(x_coords, y_coords, c=values, cmap=cmap, s=100, edgecolors='black')
    
    # Add value labels
    for i, (x, y, val) in enumerate(zip(x_coords, y_coords, values)):
        ax.annotate(f"{val:.1f}", (x, y), xytext=(0, 0), textcoords="offset points",
                   ha='center', va='center', color='white' if normalized_values[i] > 0.5 else 'black',
                   fontweight='bold')
    
    # Set axis labels and title
    ax.set_xlabel('Time Step')
    ax.set_ylabel('Down Movements')
    ax.set_title(title)
    
    # Add colorbar
    cbar = plt.colorbar(scatter)
    cbar.set_label('Value ($ million)')
    
    # Set ticks
    ax.set_xticks(range(n_steps + 1))
    ax.set_yticks(range(n_steps + 1))
    
    # Adjust layout
    plt.tight_layout()
    
    return fig, ax

# Plot the asset price tree
fig1, ax1 = plot_binomial_tree(asset_tree, 'Project Value Tree ($ million)')
plt.show()

# Plot the option value tree
fig2, ax2 = plot_binomial_tree(option_tree, 'Real Option Value Tree ($ million)')
plt.show()

# Now let's also implement the Black-Scholes model for continuous-time valuation
def black_scholes_call(S, K, T, r, sigma):
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)

# Calculate using Black-Scholes
bs_option_value = black_scholes_call(S0, K, T, r, sigma)
print(f"Black-Scholes Real Option Value: ${bs_option_value:.2f} million")

# Sensitivity analysis - varying volatility
volatilities = np.linspace(0.1, 0.6, 20)
bs_option_values = [black_scholes_call(S0, K, T, r, vol) for vol in volatilities]

plt.figure(figsize=(12, 6))
plt.plot(volatilities, bs_option_values, 'b-', linewidth=2)
plt.grid(True)
plt.xlabel('Volatility (σ)')
plt.ylabel('Real Option Value ($ million)')
plt.title('Sensitivity of Real Option Value to Volatility')
plt.show()

# Sensitivity analysis - varying underlying asset value
asset_values = np.linspace(50, 150, 20)
bs_option_values = [black_scholes_call(S, K, T, r, sigma) for S in asset_values]
npv_values = [S - K for S in asset_values]

plt.figure(figsize=(12, 6))
plt.plot(asset_values, bs_option_values, 'b-', linewidth=2, label='Real Option Value')
plt.plot(asset_values, npv_values, 'r--', linewidth=2, label='Traditional NPV')
plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
plt.axvline(x=K, color='k', linestyle='--', alpha=0.3, label='Investment Cost')
plt.grid(True)
plt.xlabel('Present Value of Cash Flows ($ million)')
plt.ylabel('Value ($ million)')
plt.title('Real Option Value vs. Traditional NPV')
plt.legend()
plt.show()

# Sensitivity to time to expiration
times = np.linspace(0.5, 10, 20)
bs_option_values = [black_scholes_call(S0, K, t, r, sigma) for t in times]

plt.figure(figsize=(12, 6))
plt.plot(times, bs_option_values, 'g-', linewidth=2)
plt.grid(True)
plt.xlabel('Time to Expiration (years)')
plt.ylabel('Real Option Value ($ million)')
plt.title('Impact of Option Expiration Time on Real Option Value')
plt.show()

# Create a 3D surface plot showing option value as a function of PV and volatility
PV_range = np.linspace(70, 130, 30)
vol_range = np.linspace(0.1, 0.5, 30)
PV_grid, vol_grid = np.meshgrid(PV_range, vol_range)
option_values = np.zeros_like(PV_grid)

for i in range(len(vol_range)):
    for j in range(len(PV_range)):
        option_values[i, j] = black_scholes_call(PV_range[j], K, T, r, vol_range[i])

fig = plt.figure(figsize=(14, 10))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(PV_grid, vol_grid, option_values, cmap=cm.viridis, 
                       linewidth=0, antialiased=True)
ax.set_xlabel('Present Value of Cash Flows ($ million)')
ax.set_ylabel('Volatility (σ)')
ax.set_zlabel('Real Option Value ($ million)')
ax.set_title('3D Surface of Real Option Value')
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)
plt.show()

# Create a heatmap
heatmap_data = pd.DataFrame(option_values, index=vol_range, columns=PV_range)
plt.figure(figsize=(14, 10))
sns.heatmap(heatmap_data, cmap='viridis', annot=False)
plt.xlabel('Present Value of Cash Flows ($ million)')
plt.ylabel('Volatility (σ)')
plt.title('Heatmap of Real Option Value')
plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import norm
from datetime import datetime, timedelta
import matplotlib.dates as mdates

# Set styling for plots
plt.style.use('ggplot')
sns.set_theme(style="whitegrid")

class AirlineRevenueSim:
    def __init__(self, total_seats=100, days_to_departure=60, 
                 base_price=300, demand_model='normal',
                 price_sensitivity=-0.01, random_seed=42):
        """
        Initialize the airline revenue management simulation.
        
        Parameters:
        -----------
        total_seats : int
            Total capacity of the flight
        days_to_departure : int
            Number of days before departure to start selling
        base_price : float
            Starting price point for the ticket
        demand_model : str
            Type of demand model to use ('normal' for normally distributed demand)
        price_sensitivity : float
            How sensitive demand is to price changes (negative value)
        random_seed : int
            Seed for reproducibility
        """
        np.random.seed(random_seed)
        
        self.total_seats = total_seats
        self.days_to_departure = days_to_departure
        self.base_price = base_price
        self.price_sensitivity = price_sensitivity
        
        # Create a timeline from sales start to departure
        self.timeline = pd.date_range(
            end=datetime.now() + timedelta(days=5),  # departure date
            periods=days_to_departure,
            freq='D'
        )
        
        # Generate the base demand curve - demand typically increases as departure approaches
        if demand_model == 'normal':
            # Base demand curve follows a normal distribution with peak near departure
            mu = 0.8 * days_to_departure  # Peak demand closer to departure
            sigma = days_to_departure / 4
            base_demand = norm.pdf(np.arange(days_to_departure), mu, sigma)
            self.base_demand = base_demand / np.max(base_demand) * 30  # Scale to reasonable daily demand
        
        # Initialize remaining inventory
        self.remaining_seats = total_seats
        
        # Storage for simulation results
        self.results = {
            'date': [],
            'price': [],
            'demand': [],
            'sales': [],
            'revenue': [],
            'remaining_seats': []
        }
    
    def calculate_demand(self, price, day_index):
        """Calculate expected demand based on price and day."""
        # Base demand for this day
        base_demand = self.base_demand[day_index]
        
        # Apply price sensitivity - higher prices reduce demand
        price_effect = np.exp(self.price_sensitivity * (price - self.base_price))
        
        # Add some randomness to simulate real-world variability
        randomness = np.random.normal(1, 0.2)
        
        # Calculate expected demand
        expected_demand = base_demand * price_effect * randomness
        
        # Ensure non-negative demand
        return max(0, expected_demand)
    
    def run_static_pricing(self):
        """Run simulation with static pricing (control case)."""
        price = self.base_price
        
        for day in range(self.days_to_departure):
            # Calculate demand at current price
            demand = self.calculate_demand(price, day)
            
            # Actual sales are limited by remaining inventory
            sales = min(demand, self.remaining_seats)
            
            # Update remaining inventory
            self.remaining_seats -= sales
            
            # Record results
            self.results['date'].append(self.timeline[day])
            self.results['price'].append(price)
            self.results['demand'].append(demand)
            self.results['sales'].append(sales)
            self.results['revenue'].append(sales * price)
            self.results['remaining_seats'].append(self.remaining_seats)
            
            # If we've sold out, stop simulation
            if self.remaining_seats <= 0:
                break
                
        return pd.DataFrame(self.results)
    
    def run_dynamic_pricing(self, load_factor_thresholds=[0.7, 0.5, 0.3], 
                           price_adjustments=[1.0, 1.3, 1.6, 2.0]):
        """
        Run simulation with dynamic pricing based on load factor.
        
        Parameters:
        -----------
        load_factor_thresholds : list
            Load factor levels at which to adjust prices
        price_adjustments : list
            Price multipliers to apply at each threshold (should be one more than thresholds)
        """
        for day in range(self.days_to_departure):
            # Calculate days remaining
            days_remaining = self.days_to_departure - day - 1
            
            # Calculate current load factor
            load_factor = 1 - (self.remaining_seats / self.total_seats)
            
            # Determine price based on target sell-out
            target_load_factor = (day + 1) / self.days_to_departure
            
            # Adjust price based on current load factor vs target
            if days_remaining > 0:  # Avoid division by zero
                price_multiplier = 1.0
                
                # If we're ahead of target, increase prices
                if load_factor > target_load_factor:
                    # Find which threshold we've crossed
                    for i, threshold in enumerate(load_factor_thresholds):
                        if load_factor > threshold:
                            price_multiplier = price_adjustments[i]
                            break
                    else:
                        price_multiplier = price_adjustments[-1]
                
                # Apply time-based adjustment: prices increase as departure approaches
                time_factor = 1 + (0.5 * (1 - days_remaining / self.days_to_departure))
                
                price = self.base_price * price_multiplier * time_factor
            else:
                # Last day - final price adjustment
                if load_factor > 0.9:
                    price = self.base_price * 2.2  # Premium for last-minute travelers
                else:
                    price = self.base_price * 0.9  # Discount to sell remaining seats
            
            # Calculate demand at current price
            demand = self.calculate_demand(price, day)
            
            # Actual sales are limited by remaining inventory
            sales = min(demand, self.remaining_seats)
            
            # Update remaining inventory
            self.remaining_seats -= sales
            
            # Record results
            self.results['date'].append(self.timeline[day])
            self.results['price'].append(price)
            self.results['demand'].append(demand)
            self.results['sales'].append(sales)
            self.results['revenue'].append(sales * price)
            self.results['remaining_seats'].append(self.remaining_seats)
            
            # If we've sold out, stop simulation
            if self.remaining_seats <= 0:
                break
                
        return pd.DataFrame(self.results)

# Run simulations for both pricing strategies
np.random.seed(42)  # For reproducibility

# Static pricing simulation
static_sim = AirlineRevenueSim(total_seats=150, days_to_departure=60, base_price=300)
static_results = static_sim.run_static_pricing()

# Reset and run dynamic pricing simulation
np.random.seed(42)  # Use same seed for fair comparison
dynamic_sim = AirlineRevenueSim(total_seats=150, days_to_departure=60, base_price=300)
dynamic_results = dynamic_sim.run_dynamic_pricing()

# Create comparison visualizations
fig, axs = plt.subplots(2, 2, figsize=(16, 12))

# Plot 1: Price over time
axs[0, 0].plot(static_results['date'], static_results['price'], label='Static Pricing')
axs[0, 0].plot(dynamic_results['date'], dynamic_results['price'], label='Dynamic Pricing')
axs[0, 0].set_title('Ticket Price Over Time', fontsize=14)
axs[0, 0].set_xlabel('Date')
axs[0, 0].set_ylabel('Price ($)')
axs[0, 0].legend()
axs[0, 0].xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))

# Plot 2: Remaining seats over time
axs[0, 1].plot(static_results['date'], static_results['remaining_seats'], label='Static Pricing')
axs[0, 1].plot(dynamic_results['date'], dynamic_results['remaining_seats'], label='Dynamic Pricing')
axs[0, 1].set_title('Remaining Seats Over Time', fontsize=14)
axs[0, 1].set_xlabel('Date')
axs[0, 1].set_ylabel('Seats')
axs[0, 1].legend()
axs[0, 1].xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))

# Plot 3: Daily sales
axs[1, 0].bar(np.arange(len(static_results)) - 0.2, static_results['sales'], width=0.4, label='Static Pricing')
axs[1, 0].bar(np.arange(len(dynamic_results)) + 0.2, dynamic_results['sales'], width=0.4, label='Dynamic Pricing')
axs[1, 0].set_title('Daily Ticket Sales', fontsize=14)
axs[1, 0].set_xlabel('Days Before Departure')
axs[1, 0].set_ylabel('Number of Tickets Sold')
axs[1, 0].set_xticks(np.arange(0, len(static_results), 5))
axs[1, 0].set_xticklabels([f"{60-i}" for i in range(0, len(static_results), 5)])
axs[1, 0].legend()

# Plot 4: Cumulative revenue
static_cumulative = static_results['revenue'].cumsum()
dynamic_cumulative = dynamic_results['revenue'].cumsum()

axs[1, 1].plot(static_results['date'], static_cumulative, label='Static Pricing')
axs[1, 1].plot(dynamic_results['date'], dynamic_cumulative, label='Dynamic Pricing')
axs[1, 1].set_title('Cumulative Revenue', fontsize=14)
axs[1, 1].set_xlabel('Date')
axs[1, 1].set_ylabel('Revenue ($)')
axs[1, 1].legend()
axs[1, 1].xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))

plt.tight_layout()
plt.show()

# Print summary statistics
print("Static Pricing Summary:")
print(f"Total Revenue: ${static_results['revenue'].sum():.2f}")
print(f"Average Ticket Price: ${static_results['price'].mean():.2f}")
print(f"Tickets Sold: {static_results['sales'].sum():.0f}")
print(f"Remaining Seats: {static_results['remaining_seats'].iloc[-1]:.0f}")
print(f"Load Factor: {100 * (1 - static_results['remaining_seats'].iloc[-1] / 150):.1f}%")

print("\nDynamic Pricing Summary:")
print(f"Total Revenue: ${dynamic_results['revenue'].sum():.2f}")
print(f"Average Ticket Price: ${dynamic_results['price'].mean():.2f}")
print(f"Tickets Sold: {dynamic_results['sales'].sum():.0f}")
print(f"Remaining Seats: {dynamic_results['remaining_seats'].iloc[-1]:.0f}")
print(f"Load Factor: {100 * (1 - dynamic_results['remaining_seats'].iloc[-1] / 150):.1f}%")
print(f"Revenue Improvement: {100 * (dynamic_results['revenue'].sum() / static_results['revenue'].sum() - 1):.1f}%")

# Alternative model: Optimization-based approach
# Create a function to find the optimal price for each day based on expected demand

def optimize_prices(sim, remaining_days, remaining_seats):
    """
    Find optimal price for each day to maximize expected revenue
    using a simple grid search approach.
    """
    prices = np.linspace(sim.base_price * 0.5, sim.base_price * 3, 50)
    optimal_prices = []
    
    for day in range(remaining_days):
        day_index = sim.days_to_departure - remaining_days + day
        expected_revenues = []
        
        # For each price point, calculate expected revenue
        for price in prices:
            expected_demand = sim.calculate_demand(price, day_index)
            expected_sales = min(expected_demand, remaining_seats)
            expected_revenue = expected_sales * price
            expected_revenues.append(expected_revenue)
        
        # Find price that maximizes revenue
        optimal_price = prices[np.argmax(expected_revenues)]
        optimal_prices.append(optimal_price)
        
        # Update remaining seats for next iteration
        expected_demand = sim.calculate_demand(optimal_price, day_index)
        expected_sales = min(expected_demand, remaining_seats)
        remaining_seats -= expected_sales
        
        if remaining_seats <= 0:
            break
    
    return optimal_prices

# Run optimization-based dynamic pricing
np.random.seed(42)  # Reset seed
optim_sim = AirlineRevenueSim(total_seats=150, days_to_departure=60, base_price=300)

# Calculate optimal prices for the entire booking period
optimal_prices = optimize_prices(optim_sim, optim_sim.days_to_departure, optim_sim.total_seats)

# Now simulate sales with these optimal prices
optim_results = {
    'date': [],
    'price': [],
    'demand': [],
    'sales': [],
    'revenue': [],
    'remaining_seats': []
}

remaining_seats = optim_sim.total_seats

for day in range(len(optimal_prices)):
    price = optimal_prices[day]
    demand = optim_sim.calculate_demand(price, day)
    sales = min(demand, remaining_seats)
    remaining_seats -= sales
    
    optim_results['date'].append(optim_sim.timeline[day])
    optim_results['price'].append(price)
    optim_results['demand'].append(demand)
    optim_results['sales'].append(sales)
    optim_results['revenue'].append(sales * price)
    optim_results['remaining_seats'].append(remaining_seats)
    
    if remaining_seats <= 0:
        break

optim_results = pd.DataFrame(optim_results)

# Compare all three approaches
fig, axs = plt.subplots(2, 2, figsize=(16, 12))

# Plot 1: Price over time
axs[0, 0].plot(static_results['date'], static_results['price'], label='Static Pricing')
axs[0, 0].plot(dynamic_results['date'], dynamic_results['price'], label='Rule-Based Dynamic')
axs[0, 0].plot(optim_results['date'], optim_results['price'], label='Optimization-Based')
axs[0, 0].set_title('Ticket Price Over Time', fontsize=14)
axs[0, 0].set_xlabel('Date')
axs[0, 0].set_ylabel('Price ($)')
axs[0, 0].legend()
axs[0, 0].xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))

# Plot 2: Remaining seats over time
axs[0, 1].plot(static_results['date'], static_results['remaining_seats'], label='Static Pricing')
axs[0, 1].plot(dynamic_results['date'], dynamic_results['remaining_seats'], label='Rule-Based Dynamic')
axs[0, 1].plot(optim_results['date'], optim_results['remaining_seats'], label='Optimization-Based')
axs[0, 1].set_title('Remaining Seats Over Time', fontsize=14)
axs[0, 1].set_xlabel('Date')
axs[0, 1].set_ylabel('Seats')
axs[0, 1].legend()
axs[0, 1].xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))

# Plot 3: Daily sales comparison
days_to_show = min(len(static_results), len(dynamic_results), len(optim_results))
axs[1, 0].plot(range(days_to_show), static_results['sales'][:days_to_show], label='Static Pricing')
axs[1, 0].plot(range(days_to_show), dynamic_results['sales'][:days_to_show], label='Rule-Based Dynamic')
axs[1, 0].plot(range(days_to_show), optim_results['sales'][:days_to_show], label='Optimization-Based')
axs[1, 0].set_title('Daily Ticket Sales', fontsize=14)
axs[1, 0].set_xlabel('Days Before Departure')
axs[1, 0].set_ylabel('Number of Tickets Sold')
axs[1, 0].set_xticks(np.arange(0, days_to_show, 5))
axs[1, 0].set_xticklabels([f"{60-i}" for i in range(0, days_to_show, 5)])
axs[1, 0].legend()

# Plot 4: Cumulative revenue
static_cumulative = static_results['revenue'].cumsum()
dynamic_cumulative = dynamic_results['revenue'].cumsum()
optim_cumulative = optim_results['revenue'].cumsum()

axs[1, 1].plot(static_results['date'][:days_to_show], static_cumulative[:days_to_show], label='Static Pricing')
axs[1, 1].plot(dynamic_results['date'][:days_to_show], dynamic_cumulative[:days_to_show], label='Rule-Based Dynamic')
axs[1, 1].plot(optim_results['date'][:days_to_show], optim_cumulative[:days_to_show], label='Optimization-Based')
axs[1, 1].set_title('Cumulative Revenue', fontsize=14)
axs[1, 1].set_xlabel('Date')
axs[1, 1].set_ylabel('Revenue ($)')
axs[1, 1].legend()
axs[1, 1].xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))

plt.tight_layout()
plt.show()

# Print comparison of all three approaches
print("Static Pricing Summary:")
print(f"Total Revenue: ${static_results['revenue'].sum():.2f}")
print(f"Average Ticket Price: ${static_results['price'].mean():.2f}")
print(f"Tickets Sold: {static_results['sales'].sum():.0f}")
print(f"Load Factor: {100 * (1 - static_results['remaining_seats'].iloc[-1] / 150):.1f}%")

print("\nRule-Based Dynamic Pricing Summary:")
print(f"Total Revenue: ${dynamic_results['revenue'].sum():.2f}")
print(f"Average Ticket Price: ${dynamic_results['price'].mean():.2f}")
print(f"Tickets Sold: {dynamic_results['sales'].sum():.0f}")
print(f"Load Factor: {100 * (1 - dynamic_results['remaining_seats'].iloc[-1] / 150):.1f}%")
print(f"Revenue Improvement: {100 * (dynamic_results['revenue'].sum() / static_results['revenue'].sum() - 1):.1f}%")

print("\nOptimization-Based Pricing Summary:")
print(f"Total Revenue: ${optim_results['revenue'].sum():.2f}")
print(f"Average Ticket Price: ${optim_results['price'].mean():.2f}")
print(f"Tickets Sold: {optim_results['sales'].sum():.0f}")
print(f"Load Factor: {100 * (1 - optim_results['remaining_seats'].iloc[-1] / 150):.1f}%")
print(f"Revenue Improvement: {100 * (optim_results['revenue'].sum() / static_results['revenue'].sum() - 1):.1f}%")