import networkx as nx
import matplotlib.pyplot as plt

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

# Add edges and their capacities
edges = [
    (0, 1, 16),
    (0, 2, 13),
    (1, 2, 10),
    (1, 3, 12),
    (2, 1, 4),
    (2, 4, 14),
    (3, 2, 9),
    (3, 5, 20),
    (4, 3, 7),
    (4, 5, 4)
]

for u, v, c in edges:
    G.add_edge(u, v, capacity=c)

# Compute the maximum flow using Edmonds-Karp algorithm
flow_value, flow_dict = nx.maximum_flow(G, 0, 5, flow_func=nx.algorithms.flow.edmonds_karp)

print("Maximum flow:", flow_value)
print("Flow per edge:")
for u in flow_dict:
    for v in flow_dict[u]:
        if flow_dict[u][v] > 0:
            print(f"  {u} -> {v}: {flow_dict[u][v]}")

# Draw the graph with flow values
pos = nx.spring_layout(G, seed=42)  # layout for consistent node positioning
edge_labels = {
    (u, v): f"{flow_dict[u][v]}/{G[u][v]['capacity']}"
    for u, v in G.edges()
}

plt.figure(figsize=(10, 6))
nx.draw_networkx(G, pos, with_labels=True, node_size=800, node_color='lightblue')
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
nx.draw_networkx_edges(G, pos, width=2, edge_color='gray')
plt.title("Maximum Flow in the Network (flow/capacity)")
plt.axis('off')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Parameters
np.random.seed(42)
T = 50  # time steps
alpha = 0.3  # smoothing parameter

# Generate customer demand: small random walk
customer_demand = np.round(np.cumsum(np.random.normal(0, 1, T)) + 20).astype(int)
customer_demand = np.clip(customer_demand, 15, 30)

# Initialize arrays
forecast_demand = np.zeros(T)
retailer_orders = np.zeros(T)
wholesaler_orders = np.zeros(T)
manufacturer_orders = np.zeros(T)

# Initial forecast
forecast_demand[0] = customer_demand[0]

# Forecast using exponential smoothing
for t in range(1, T):
    forecast_demand[t] = alpha * customer_demand[t-1] + (1 - alpha) * forecast_demand[t-1]

# Order policies
retailer_orders = forecast_demand + 2  # +2 as safety stock
wholesaler_orders = np.roll(retailer_orders, 1)
manufacturer_orders = np.roll(wholesaler_orders, 1)
wholesaler_orders[0] = retailer_orders[0]
manufacturer_orders[0] = wholesaler_orders[0]

# Plot
plt.figure(figsize=(12, 6))
plt.plot(customer_demand, label='Customer Demand', marker='o')
plt.plot(retailer_orders, label='Retailer Orders', marker='o')
plt.plot(wholesaler_orders, label='Wholesaler Orders', marker='o')
plt.plot(manufacturer_orders, label='Manufacturer Orders', marker='o')
plt.title('Bullwhip Effect Simulation with Smoothing')
plt.xlabel('Time')
plt.ylabel('Order Quantity')
plt.legend()
plt.grid(True)
plt.show()

customer_demand = np.round(np.cumsum(np.random.normal(0, 1, T)) + 20).astype(int)
customer_demand = np.clip(customer_demand, 15, 30)

forecast_demand[t] = alpha * customer_demand[t-1] + (1 - alpha) * forecast_demand[t-1]

retailer_orders = forecast_demand + 2
wholesaler_orders = np.roll(retailer_orders, 1)
manufacturer_orders = np.roll(wholesaler_orders, 1)

# Importing necessary libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import minimize
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# Set the aesthetic style of the plots
plt.style.use('ggplot')
sns.set_palette("colorblind")

# Define the demand function for multiple products
def demand_function(prices, base_demand, price_sensitivity, cross_elasticity_matrix):
    """
    Calculate the demand for each product given their prices and elasticities.
    
    Parameters:
    - prices: Array of prices for each product
    - base_demand: Base demand when price is at reference level
    - price_sensitivity: Own-price elasticity for each product
    - cross_elasticity_matrix: Matrix of cross-price elasticities
    
    Returns:
    - Array of demand values for each product
    """
    # Initialize demand with base demand
    demand = base_demand.copy()
    
    # Apply own-price effects
    for i in range(len(prices)):
        demand[i] *= np.exp(-price_sensitivity[i] * prices[i])
    
    # Apply cross-price effects
    for i in range(len(prices)):
        for j in range(len(prices)):
            if i != j:
                # Positive cross-elasticity means products are substitutes
                # Negative cross-elasticity means products are complements
                demand[i] *= np.exp(cross_elasticity_matrix[i][j] * prices[j])
    
    return demand

# Define the profit function to maximize
def profit_function(prices, base_demand, price_sensitivity, cross_elasticity_matrix, costs):
    """
    Calculate the total profit across all products.
    
    Parameters:
    - prices: Array of prices for each product
    - base_demand: Base demand when price is at reference level
    - price_sensitivity: Own-price elasticity for each product
    - cross_elasticity_matrix: Matrix of cross-price elasticities
    - costs: Variable costs for each product
    
    Returns:
    - Total profit (negative for minimization algorithm)
    """
    demand = demand_function(prices, base_demand, price_sensitivity, cross_elasticity_matrix)
    profit = np.sum((prices - costs) * demand)
    # Return negative profit for minimization (scipy.optimize.minimize minimizes the function)
    return -profit

# Define a function to analyze price elasticity across a range
def price_elasticity_analysis(base_price, product_index, price_range_percentage, 
                            base_demand, price_sensitivity, cross_elasticity_matrix, costs):
    """
    Analyze how profit and demand change as one product's price changes.
    
    Parameters:
    - base_price: Base prices for all products
    - product_index: Index of the product to vary
    - price_range_percentage: Percentage range around base price to analyze
    - Other parameters as defined in previous functions
    
    Returns:
    - DataFrame with price, demand, and profit data
    """
    min_price = base_price[product_index] * (1 - price_range_percentage/100)
    max_price = base_price[product_index] * (1 + price_range_percentage/100)
    price_points = np.linspace(min_price, max_price, 100)
    
    results = []
    for price in price_points:
        # Create a copy of the base price and modify the selected product's price
        current_prices = base_price.copy()
        current_prices[product_index] = price
        
        # Calculate demand and profit
        demand = demand_function(current_prices, base_demand, price_sensitivity, cross_elasticity_matrix)
        profit = -profit_function(current_prices, base_demand, price_sensitivity, 
                                cross_elasticity_matrix, costs)
        
        # Store results
        results.append({
            'Price': price,
            'Demand': demand[product_index],
            'Total Demand': np.sum(demand),
            'Profit': profit
        })
    
    return pd.DataFrame(results)

# Set up an example for a company offering three subscription tiers
# Product 0: Basic tier
# Product 1: Premium tier
# Product 2: Enterprise tier

# Initial parameters
num_products = 3
initial_prices = np.array([10.0, 25.0, 50.0])  # Initial prices
base_demand = np.array([5000, 2000, 1000])     # Base demand at reference prices
costs = np.array([2.0, 5.0, 10.0])             # Variable costs per unit

# Price sensitivity (how demand responds to product's own price)
# Higher values mean demand drops more as price increases
price_sensitivity = np.array([0.15, 0.1, 0.05])

# Cross-price elasticity matrix
# Positive values: substitutes (if price of j increases, demand for i increases)
# Negative values: complements (if price of j increases, demand for i decreases)
cross_elasticity_matrix = np.array([
    [0.0,  0.08, 0.04],  # Impact of products on Basic tier
    [0.03, 0.0,  0.06],  # Impact of products on Premium tier
    [0.01, 0.05, 0.0]    # Impact of products on Enterprise tier
])

# Run the optimization
result = minimize(
    profit_function,
    initial_prices,
    args=(base_demand, price_sensitivity, cross_elasticity_matrix, costs),
    method='SLSQP',
    bounds=[(cost*1.1, cost*10) for cost in costs],  # Price must be at least 10% above cost
    options={'disp': True}
)

# Extract optimized prices
optimal_prices = result.x
print(f"Optimal prices: {optimal_prices}")

# Calculate demand and profit at optimal prices
optimal_demand = demand_function(optimal_prices, base_demand, price_sensitivity, cross_elasticity_matrix)
optimal_profit = -profit_function(optimal_prices, base_demand, price_sensitivity, 
                               cross_elasticity_matrix, costs)

print(f"Demand at optimal prices: {optimal_demand}")
print(f"Total profit at optimal prices: ${optimal_profit:.2f}")

# Analyze how changing each product's price affects overall profit
price_range = 50  # Analyze prices ±50% around optimal
product_names = ["Basic Tier", "Premium Tier", "Enterprise Tier"]

# Create visualizations

# 1. Bar chart comparing initial vs. optimal prices
plt.figure(figsize=(12, 6))
width = 0.35
x = np.arange(num_products)
plt.bar(x - width/2, initial_prices, width, label='Initial Prices')
plt.bar(x + width/2, optimal_prices, width, label='Optimal Prices')
plt.xlabel('Product')
plt.ylabel('Price ($)')
plt.title('Comparison of Initial vs. Optimal Prices')
plt.xticks(x, product_names)
plt.legend()
plt.tight_layout()
plt.savefig('price_comparison.png', dpi=300)
plt.show()

# 2. Plot profit curves for each product around optimal price
plt.figure(figsize=(16, 6))

for i in range(num_products):
    plt.subplot(1, 3, i+1)
    analysis_df = price_elasticity_analysis(
        optimal_prices, i, price_range, base_demand, 
        price_sensitivity, cross_elasticity_matrix, costs
    )
    
    # Find the optimal price within our analyzed range
    max_profit_idx = analysis_df['Profit'].idxmax()
    max_profit_price = analysis_df.loc[max_profit_idx, 'Price']
    
    plt.plot(analysis_df['Price'], analysis_df['Profit'])
    plt.axvline(x=optimal_prices[i], color='red', linestyle='--', 
                label=f'Optimal: ${optimal_prices[i]:.2f}')
    plt.title(f'Profit Curve for {product_names[i]}')
    plt.xlabel('Price ($)')
    plt.ylabel('Total Profit ($)')
    plt.legend()
    plt.grid(True)

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

# 3. Plot demand curves
plt.figure(figsize=(16, 6))

for i in range(num_products):
    plt.subplot(1, 3, i+1)
    analysis_df = price_elasticity_analysis(
        optimal_prices, i, price_range, base_demand, 
        price_sensitivity, cross_elasticity_matrix, costs
    )
    
    plt.plot(analysis_df['Price'], analysis_df['Demand'])
    plt.axvline(x=optimal_prices[i], color='red', linestyle='--', 
                label=f'Optimal: ${optimal_prices[i]:.2f}')
    plt.title(f'Demand Curve for {product_names[i]}')
    plt.xlabel('Price ($)')
    plt.ylabel('Quantity Demanded')
    plt.legend()
    plt.grid(True)

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

# 4. Create a 3D surface plot for two products to visualize the profit landscape
# We'll vary prices for Basic and Premium tiers while keeping Enterprise at optimal
plt.figure(figsize=(12, 10))
ax = plt.axes(projection='3d')

# Create a grid of prices for products 0 and 1
price_range_percent = 30
p0_range = np.linspace(optimal_prices[0] * (1 - price_range_percent/100), 
                     optimal_prices[0] * (1 + price_range_percent/100), 20)
p1_range = np.linspace(optimal_prices[1] * (1 - price_range_percent/100), 
                     optimal_prices[1] * (1 + price_range_percent/100), 20)
P0, P1 = np.meshgrid(p0_range, p1_range)
profit_values = np.zeros(P0.shape)

# Calculate profit for each price combination
for i in range(len(p0_range)):
    for j in range(len(p1_range)):
        current_prices = optimal_prices.copy()
        current_prices[0] = P0[i, j]
        current_prices[1] = P1[i, j]
        profit_values[i, j] = -profit_function(current_prices, base_demand, 
                                             price_sensitivity, cross_elasticity_matrix, costs)

# Create the surface plot
surf = ax.plot_surface(P0, P1, profit_values, cmap=cm.coolwarm,
                       linewidth=0, antialiased=True, alpha=0.8)

# Mark the optimal point
ax.scatter([optimal_prices[0]], [optimal_prices[1]], 
           [-profit_function(optimal_prices, base_demand, price_sensitivity, 
                          cross_elasticity_matrix, costs)],
           color='black', s=100, label='Optimal Price Point')

ax.set_xlabel('Basic Tier Price ($)')
ax.set_ylabel('Premium Tier Price ($)')
ax.set_zlabel('Profit ($)')
ax.set_title('Profit Landscape for Basic and Premium Tiers')
plt.colorbar(surf, ax=ax, shrink=0.5, aspect=5, label='Profit ($)')
plt.savefig('profit_landscape_3d.png', dpi=300)
plt.show()

# 5. Create a heatmap for the same profit landscape
plt.figure(figsize=(10, 8))
profit_df = pd.DataFrame(profit_values, index=p0_range, columns=p1_range)
sns.heatmap(profit_df, cmap='viridis', annot=False)
plt.xlabel('Premium Tier Price ($)')
plt.ylabel('Basic Tier Price ($)')
plt.title('Profit Heatmap for Basic vs Premium Tier Pricing')
plt.tight_layout()
plt.savefig('profit_heatmap.png', dpi=300)
plt.show()

# 6. Create a summary table of results
summary_data = {
    'Product': product_names,
    'Initial Price ($)': initial_prices,
    'Optimal Price ($)': optimal_prices,
    'Price Change (%)': (optimal_prices - initial_prices) / initial_prices * 100,
    'Demand at Optimal': optimal_demand,
    'Unit Profit ($)': optimal_prices - costs,
    'Total Profit ($)': (optimal_prices - costs) * optimal_demand
}

summary_df = pd.DataFrame(summary_data)
summary_df['Price Change (%)'] = summary_df['Price Change (%)'].round(2)
summary_df['Unit Profit ($)'] = summary_df['Unit Profit ($)'].round(2)
summary_df['Total Profit ($)'] = summary_df['Total Profit ($)'].round(2)

print("\nSummary of Pricing Optimization Results:")
print(summary_df)

# Display the cross-elasticity matrix as a heatmap to show relationships between products
plt.figure(figsize=(8, 6))
sns.heatmap(cross_elasticity_matrix, annot=True, cmap='coolwarm', 
            xticklabels=product_names, yticklabels=product_names)
plt.title('Cross-Price Elasticity Matrix')
plt.xlabel('Price Change for Product')
plt.ylabel('Demand Impact on Product')
plt.tight_layout()
plt.savefig('cross_elasticity_heatmap.png', dpi=300)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun']
demand = np.array([80, 100, 60, 90, 110, 70])

!pip install pulp

from pulp import LpMinimize, LpProblem, LpVariable, lpSum, LpStatus

# Problem definition
model = LpProblem("Production_Planning", LpMinimize)

n_months = len(demand)
production = [LpVariable(f"x_{t}", lowBound=0, upBound=100) for t in range(n_months)]
inventory = [LpVariable(f"I_{t}", lowBound=0) for t in range(n_months)]

# Objective function: production cost + inventory cost
production_cost = lpSum(10 * production[t] for t in range(n_months))
inventory_cost = lpSum(2 * inventory[t] for t in range(n_months))
model += production_cost + inventory_cost

# Inventory balance constraints
for t in range(n_months):
    if t == 0:
        model += inventory[t] == production[t] - demand[t]
    else:
        model += inventory[t] == inventory[t-1] + production[t] - demand[t]

# Solve
model.solve()

# Output
print(f"Status: {LpStatus[model.status]}")
for t in range(n_months):
    print(f"{months[t]} - Production: {production[t].value():.0f}, Inventory: {inventory[t].value():.0f}")

prod_vals = [production[t].value() for t in range(n_months)]
inv_vals = [inventory[t].value() for t in range(n_months)]

plt.figure(figsize=(10, 5))
plt.plot(months, demand, label="Demand", marker='o')
plt.plot(months, prod_vals, label="Production", marker='s')
plt.plot(months, inv_vals, label="Inventory", marker='^')
plt.title("Optimal Production Plan")
plt.xlabel("Month")
plt.ylabel("Units")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Generate random delivery locations
np.random.seed(42)
num_points = 50
X = np.random.rand(num_points, 2) * 100  # Coordinates within 100x100 grid

# Plot delivery points
plt.figure(figsize=(6, 6))
plt.scatter(X[:, 0], X[:, 1], c='blue')
plt.title("Delivery Points")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.grid(True)
plt.show()

from sklearn.cluster import KMeans

num_clusters = 4
kmeans = KMeans(n_clusters=num_clusters, random_state=42)
labels = kmeans.fit_predict(X)

# Plot clustered points
plt.figure(figsize=(6, 6))
for i in range(num_clusters):
    cluster = X[labels == i]
    plt.scatter(cluster[:, 0], cluster[:, 1], label=f'Zone {i+1}')
plt.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1], 
            c='black', marker='x', s=100, label='Centers')
plt.title("Clustered Delivery Zones")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.legend()
plt.grid(True)
plt.show()

from ortools.constraint_solver import pywrapcp, routing_enums_pb2

def create_distance_matrix(locations):
    from scipy.spatial.distance import cdist
    return cdist(locations, locations).astype(int)

def solve_tsp(distance_matrix):
    manager = pywrapcp.RoutingIndexManager(len(distance_matrix), 1, 0)
    routing = pywrapcp.RoutingModel(manager)

    def distance_callback(from_idx, to_idx):
        return distance_matrix[manager.IndexToNode(from_idx)][manager.IndexToNode(to_idx)]

    transit_callback_index = routing.RegisterTransitCallback(distance_callback)
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    search_params = pywrapcp.DefaultRoutingSearchParameters()
    search_params.first_solution_strategy = routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC

    solution = routing.SolveWithParameters(search_params)
    
    if not solution:
        return None
    
    # Extract the route
    index = routing.Start(0)
    route = []
    while not routing.IsEnd(index):
        route.append(manager.IndexToNode(index))
        index = solution.Value(routing.NextVar(index))
    route.append(manager.IndexToNode(index))
    return route

zone_routes = []
plt.figure(figsize=(8, 8))

colors = ['r', 'g', 'b', 'orange']
for zone_id in range(num_clusters):
    zone_points = X[labels == zone_id]
    if len(zone_points) < 2:
        continue  # Skip small clusters

    distance_matrix = create_distance_matrix(zone_points)
    route = solve_tsp(distance_matrix)
    if route is None:
        print(f"No route found for zone {zone_id}")
        continue

    # Save and plot
    zone_routes.append(route)
    zone_route_coords = zone_points[route]

    plt.plot(zone_route_coords[:, 0], zone_route_coords[:, 1], 
             marker='o', color=colors[zone_id], label=f'Zone {zone_id+1}')

plt.title("Optimized Delivery Routes per Zone")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.grid(True)
plt.legend()
plt.show()

import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.patches import Rectangle
from matplotlib.colors import LinearSegmentedColormap

# Define the project activities with their durations and dependencies
activities = {
    'A': {'name': 'Requirements gathering', 'duration': 3, 'predecessors': []},
    'B': {'name': 'System design', 'duration': 4, 'predecessors': ['A']},
    'C': {'name': 'Database setup', 'duration': 2, 'predecessors': ['A']},
    'D': {'name': 'Frontend development', 'duration': 5, 'predecessors': ['B']},
    'E': {'name': 'Backend development', 'duration': 6, 'predecessors': ['B', 'C']},
    'F': {'name': 'Integration', 'duration': 3, 'predecessors': ['D', 'E']},
    'G': {'name': 'Testing', 'duration': 4, 'predecessors': ['F']},
    'H': {'name': 'Deployment', 'duration': 2, 'predecessors': ['G']}
}

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

# Add nodes (activities) to the graph
for activity_id, activity_info in activities.items():
    G.add_node(activity_id, **activity_info)

# Add edges (dependencies) to the graph
for activity_id, activity_info in activities.items():
    for predecessor in activity_info['predecessors']:
        G.add_edge(predecessor, activity_id)

# Calculate earliest start and finish times (forward pass)
earliest_start = {'A': 0}  # Activity A starts at time 0
earliest_finish = {}

# Topological sort ensures we process activities in the correct order
for activity_id in nx.topological_sort(G):
    activity_info = G.nodes[activity_id]
    
    # Set earliest start time
    if activity_id not in earliest_start:
        if not activity_info['predecessors']:
            earliest_start[activity_id] = 0
        else:
            earliest_start[activity_id] = max(earliest_finish.get(pred, 0) for pred in activity_info['predecessors'])
    
    # Calculate earliest finish time
    earliest_finish[activity_id] = earliest_start[activity_id] + activity_info['duration']

# Calculate latest start and finish times (backward pass)
project_duration = max(earliest_finish.values())
latest_finish = {activity_id: project_duration for activity_id in activities}
latest_start = {}

# Process activities in reverse topological order
for activity_id in reversed(list(nx.topological_sort(G))):
    # Find successors of the current activity
    successors = list(G.successors(activity_id))
    
    if successors:
        latest_finish[activity_id] = min(latest_start.get(succ, project_duration) for succ in successors)
    
    latest_start[activity_id] = latest_finish[activity_id] - G.nodes[activity_id]['duration']

# Calculate float (slack) for each activity
float_times = {}
for activity_id in activities:
    float_times[activity_id] = latest_start[activity_id] - earliest_start[activity_id]

# Identify the critical path (activities with zero float)
critical_path = [activity_id for activity_id, float_time in float_times.items() if float_time == 0]

# Print the results
print(f"Project Duration: {project_duration} days")
print(f"Critical Path: {' -> '.join(critical_path)}")
print("\nActivity Details:")
print(f"{'Activity':<10}{'Description':<25}{'Duration':<10}{'ES':<5}{'EF':<5}{'LS':<5}{'LF':<5}{'Float':<5}{'Critical':<8}")
print("-" * 80)

for activity_id, activity_info in sorted(activities.items()):
    is_critical = activity_id in critical_path
    print(f"{activity_id:<10}{activity_info['name']:<25}{activity_info['duration']:<10}{earliest_start[activity_id]:<5}"
          f"{earliest_finish[activity_id]:<5}{latest_start[activity_id]:<5}{latest_finish[activity_id]:<5}"
          f"{float_times[activity_id]:<5}{'Yes' if is_critical else 'No':<8}")

# Visualization of the network diagram
plt.figure(figsize=(14, 8))
pos = nx.spring_layout(G, seed=42)  # Layout for the graph

# Draw the network
node_colors = ['red' if node in critical_path else 'skyblue' for node in G.nodes()]
edge_colors = ['red' if u in critical_path and v in critical_path else 'black' for u, v in G.edges()]
edge_width = [2 if u in critical_path and v in critical_path else 1 for u, v in G.edges()]

nx.draw(G, pos, with_labels=True, node_color=node_colors, edge_color=edge_colors, 
        width=edge_width, node_size=700, font_size=10, font_weight='bold')

plt.title('Project Network Diagram with Critical Path (in red)', fontsize=16)
plt.savefig('network_diagram.png', dpi=300, bbox_inches='tight')
plt.close()

# Create a Gantt chart
plt.figure(figsize=(14, 8))

# Create custom colormap: blue for normal activities, red for critical path
cmap = LinearSegmentedColormap.from_list('custom_cmap', ['skyblue', 'red'])

# Sort activities by earliest start time
sorted_activities = sorted(activities.items(), key=lambda x: earliest_start[x[0]])

# Draw the Gantt chart
y_positions = np.arange(len(activities))
y_labels = [f"{activity_id}: {info['name']}" for activity_id, info in sorted_activities]

# Draw activity bars
for i, (activity_id, info) in enumerate(sorted_activities):
    is_critical = activity_id in critical_path
    color = 'red' if is_critical else 'skyblue'
    
    # Draw the main activity bar
    plt.barh(i, info['duration'], left=earliest_start[activity_id], color=color, 
             edgecolor='black', alpha=0.8)
    
    # Draw the float/slack time (if any)
    if float_times[activity_id] > 0:
        plt.barh(i, float_times[activity_id], left=earliest_finish[activity_id], 
                 color='lightgray', alpha=0.5, edgecolor='gray', hatch='/')
    
    # Add text labels on the bars
    plt.text(earliest_start[activity_id] + info['duration']/2, i, 
             f"{activity_id} ({info['duration']}d)", 
             ha='center', va='center', color='black', fontweight='bold')

# Draw time grid
for t in range(0, project_duration + 1):
    plt.axvline(x=t, color='gray', linestyle='--', alpha=0.3)

# Set chart properties
plt.yticks(y_positions, y_labels)
plt.xlabel('Time (days)')
plt.grid(axis='x', alpha=0.3)
plt.title('Project Gantt Chart with Critical Path', fontsize=16)

# Add legend
legend_elements = [
    Rectangle((0, 0), 1, 1, color='red', alpha=0.8, label='Critical Activity'),
    Rectangle((0, 0), 1, 1, color='skyblue', alpha=0.8, label='Non-Critical Activity'),
    Rectangle((0, 0), 1, 1, color='lightgray', alpha=0.5, hatch='/', label='Float Time')
]
plt.legend(handles=legend_elements, loc='upper right')

# Show the Gantt chart
plt.tight_layout()
plt.savefig('gantt_chart.png', dpi=300, bbox_inches='tight')
plt.show()

import numpy as np
import pandas as pd
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
import seaborn as sns

# Setting a consistent visual style
plt.style.use('seaborn-v0_8-darkgrid')
sns.set_palette("viridis")

# Define the problem parameters
n_apps = 5  # Number of applications
n_resources = 2  # Types of resources: CPU and Memory

# Define resource constraints (total available)
total_resources = np.array([100, 200])  # 100 CPU units, 200 Memory units

# Define resource requirements per unit of allocation for each app
# Each row represents an app, columns are [CPU, Memory]
resource_requirements = np.array([
    [2, 3],    # App 1: 2 CPU units, 3 Memory units per allocation unit
    [1, 2],    # App 2: 1 CPU unit, 2 Memory units per allocation unit
    [3, 1],    # App 3: 3 CPU units, 1 Memory unit per allocation unit
    [2, 2],    # App 4: 2 CPU units, 2 Memory units per allocation unit
    [1, 3],    # App 5: 1 CPU unit, 3 Memory units per allocation unit
])

# Utility functions: Diminishing returns modeled with logarithmic utility
# For each app, define parameters a and b for utility function a*log(1+b*x)
utility_params = np.array([
    [10, 0.5],  # App 1: high priority but moderate scaling
    [8, 0.8],   # App 2: medium priority with good scaling
    [12, 0.3],  # App 3: highest priority but scales poorly
    [7, 0.7],   # App 4: lower priority with decent scaling
    [9, 0.6]    # App 5: medium-high priority with moderate scaling
])

# Define the utility function for a given allocation
def calculate_utility(allocation):
    utility = np.sum([
        utility_params[i, 0] * np.log(1 + utility_params[i, 1] * allocation[i])
        for i in range(n_apps)
    ])
    return utility

# Define the objective function to maximize (we minimize the negative utility)
def objective(allocation):
    return -calculate_utility(allocation)

# Define the constraint function: total resources used must be <= total available
def resource_constraint(allocation):
    # Calculate total resources used
    total_used = np.zeros(n_resources)
    for i in range(n_apps):
        total_used += allocation[i] * resource_requirements[i]
    # Return the slack in each resource (must be >= 0)
    return total_resources - total_used

# Define the constraints for the optimizer
constraints = [{
    'type': 'ineq',
    'fun': lambda x: resource_constraint(x)[i]
} for i in range(n_resources)]

# Add non-negativity constraints for allocations
bounds = [(0, None) for _ in range(n_apps)]

# Initial guess: equal allocation to all apps
initial_allocation = np.ones(n_apps) * 5

# Solve the optimization problem
result = minimize(
    objective,
    initial_allocation,
    method='SLSQP',
    bounds=bounds,
    constraints=constraints,
    options={'disp': True}
)

# Get the optimal allocation
optimal_allocation = result.x
print("Optimal Allocation:", optimal_allocation)

# Calculate the utility of the optimal allocation
optimal_utility = -result.fun
print("Optimal Utility:", optimal_utility)

# Calculate resource usage
resource_usage = np.zeros(n_resources)
for i in range(n_apps):
    resource_usage += optimal_allocation[i] * resource_requirements[i]
print("Resource Usage (CPU, Memory):", resource_usage)
print("Resource Utilization Rate:", resource_usage / total_resources * 100, "%")

# Calculate individual app utilities
app_utilities = [
    utility_params[i, 0] * np.log(1 + utility_params[i, 1] * optimal_allocation[i])
    for i in range(n_apps)
]
print("App Utilities:", app_utilities)

# Now let's visualize the results
fig, axes = plt.subplots(2, 2, figsize=(16, 12))

# Plot 1: Optimal Allocation by App
ax1 = axes[0, 0]
bars = ax1.bar(range(1, n_apps+1), optimal_allocation, color=sns.color_palette("viridis", n_apps))
ax1.set_xlabel('Application', fontsize=12)
ax1.set_ylabel('Allocation Units', fontsize=12)
ax1.set_title('Optimal Resource Allocation by Application', fontsize=14)
ax1.xaxis.set_major_locator(MaxNLocator(integer=True))
for bar in bars:
    height = bar.get_height()
    ax1.text(bar.get_x() + bar.get_width()/2., height + 0.1,
            f'{height:.2f}', ha='center', va='bottom')

# Plot 2: Utility by App
ax2 = axes[0, 1]
bars = ax2.bar(range(1, n_apps+1), app_utilities, color=sns.color_palette("viridis", n_apps))
ax2.set_xlabel('Application', fontsize=12)
ax2.set_ylabel('Utility', fontsize=12)
ax2.set_title('Utility Generated by Each Application', fontsize=14)
ax2.xaxis.set_major_locator(MaxNLocator(integer=True))
for bar in bars:
    height = bar.get_height()
    ax2.text(bar.get_x() + bar.get_width()/2., height + 0.1,
            f'{height:.2f}', ha='center', va='bottom')

# Plot 3: Resource Usage
ax3 = axes[1, 0]
resource_names = ['CPU', 'Memory']
x = np.arange(len(resource_names))
width = 0.35
bars1 = ax3.bar(x - width/2, resource_usage, width, label='Used')
bars2 = ax3.bar(x + width/2, total_resources - resource_usage, width, label='Available')
ax3.set_ylabel('Resource Units', fontsize=12)
ax3.set_title('Resource Utilization', fontsize=14)
ax3.set_xticks(x)
ax3.set_xticklabels(resource_names)
ax3.legend()
for i, bar in enumerate(bars1):
    height = bar.get_height()
    usage_pct = resource_usage[i] / total_resources[i] * 100
    ax3.text(bar.get_x() + bar.get_width()/2., height/2,
            f'{height:.1f}\n({usage_pct:.1f}%)', ha='center', va='center')

# Plot 4: Utility Curve Visualization for each app
ax4 = axes[1, 1]
x = np.linspace(0, 20, 100)
for i in range(n_apps):
    y = [utility_params[i, 0] * np.log(1 + utility_params[i, 1] * allocation) for allocation in x]
    ax4.plot(x, y, label=f'App {i+1}')
    # Mark the optimal allocation point
    opt_utility = utility_params[i, 0] * np.log(1 + utility_params[i, 1] * optimal_allocation[i])
    ax4.scatter(optimal_allocation[i], opt_utility, marker='o')
    ax4.text(optimal_allocation[i], opt_utility, f'  {optimal_allocation[i]:.2f}', va='bottom')

ax4.set_xlabel('Resource Allocation', fontsize=12)
ax4.set_ylabel('Utility', fontsize=12)
ax4.set_title('Utility Functions by Application', fontsize=14)
ax4.legend()

plt.tight_layout()
plt.show()

# Let's simulate a dynamic scenario where demand changes over time
np.random.seed(42)
time_periods = 10
demand_fluctuation = np.random.uniform(0.7, 1.3, size=(time_periods, n_apps))

# Store results for each time period
allocations_over_time = []
utilities_over_time = []
resource_usage_over_time = []

for t in range(time_periods):
    # Adjust utility parameters based on demand fluctuation
    adjusted_utility_params = utility_params.copy()
    adjusted_utility_params[:, 0] *= demand_fluctuation[t]
    
    # Define the utility function for this time period
    def calculate_utility_t(allocation):
        utility = np.sum([
            adjusted_utility_params[i, 0] * np.log(1 + adjusted_utility_params[i, 1] * allocation[i])
            for i in range(n_apps)
        ])
        return utility
    
    # Define the objective function to maximize
    def objective_t(allocation):
        return -calculate_utility_t(allocation)
    
    # Solve the optimization problem for this time period
    result_t = minimize(
        objective_t,
        initial_allocation,
        method='SLSQP',
        bounds=bounds,
        constraints=constraints,
        options={'disp': False}
    )
    
    # Store results
    allocation_t = result_t.x
    allocations_over_time.append(allocation_t)
    utilities_over_time.append(-result_t.fun)
    
    # Calculate resource usage - FIX: proper broadcasting for each app
    resource_usage_t = np.zeros(n_resources)
    for i in range(n_apps):
        # Fix here: allocate to each app individually
        resource_usage_t += allocation_t[i] * resource_requirements[i]
    resource_usage_over_time.append(resource_usage_t)
    
    # Use this period's allocation as starting point for next period
    initial_allocation = allocation_t

# Convert results to arrays for easier plotting
allocations_over_time = np.array(allocations_over_time)
utilities_over_time = np.array(utilities_over_time)
resource_usage_over_time = np.array(resource_usage_over_time)

# Visualize dynamic allocation over time
fig, axes = plt.subplots(3, 1, figsize=(14, 18))

# Plot 1: Allocation by app over time
ax1 = axes[0]
for i in range(n_apps):
    ax1.plot(range(1, time_periods+1), allocations_over_time[:, i], 
             marker='o', label=f'App {i+1}')
ax1.set_xlabel('Time Period', fontsize=12)
ax1.set_ylabel('Allocation Units', fontsize=12)
ax1.set_title('Dynamic Resource Allocation Over Time', fontsize=14)
ax1.legend()
ax1.grid(True)

# Plot 2: Total utility over time
ax2 = axes[1]
ax2.plot(range(1, time_periods+1), utilities_over_time, marker='o', color='green', linewidth=2)
ax2.set_xlabel('Time Period', fontsize=12)
ax2.set_ylabel('Total Utility', fontsize=12)
ax2.set_title('System Utility Over Time', fontsize=14)
ax2.grid(True)

# Plot 3: Resource utilization over time
ax3 = axes[2]
for i, resource_name in enumerate(resource_names):
    ax3.plot(range(1, time_periods+1), resource_usage_over_time[:, i], 
             marker='o', label=f'{resource_name} Used')
    ax3.axhline(y=total_resources[i], linestyle='--', 
                label=f'{resource_name} Limit', color=f'C{i}', alpha=0.6)
ax3.set_xlabel('Time Period', fontsize=12)
ax3.set_ylabel('Resource Units', fontsize=12)
ax3.set_title('Resource Utilization Over Time', fontsize=14)
ax3.legend()
ax3.grid(True)

plt.tight_layout()
plt.show()

# Calculate and display average metrics across all time periods
print("\n----- Dynamic Allocation Summary -----")
print(f"Average Total Utility: {np.mean(utilities_over_time):.2f}")
print(f"Average CPU Utilization: {np.mean(resource_usage_over_time[:,0])/total_resources[0]*100:.2f}%")
print(f"Average Memory Utilization: {np.mean(resource_usage_over_time[:,1])/total_resources[1]*100:.2f}%")

# Create a DataFrame to show the allocation data over time
df_allocations = pd.DataFrame(allocations_over_time, 
                             columns=[f'App {i+1}' for i in range(n_apps)])
df_allocations.index = [f'Period {i+1}' for i in range(time_periods)]
print("\nResource Allocation Over Time:")
print(df_allocations)

# Show the correlation between app demands and allocations
correlation_matrix = np.zeros((n_apps, n_apps))
for i in range(n_apps):
    for j in range(n_apps):
        correlation_matrix[i,j] = np.corrcoef(demand_fluctuation[:,i], allocations_over_time[:,j])[0,1]

plt.figure(figsize=(10, 8))
sns.heatmap(correlation_matrix, 
            annot=True, 
            xticklabels=[f'Alloc {i+1}' for i in range(n_apps)],
            yticklabels=[f'Demand {i+1}' for i in range(n_apps)],
            cmap="coolwarm")
plt.title('Correlation Between App Demand and Resource Allocation')
plt.tight_layout()
plt.show()

# 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}")