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

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

# Define our assets
assets = ['Stocks', 'Bonds', 'Real Estate', 'Commodities']
n_assets = len(assets)

# Expected returns (annual)
expected_returns = np.array([0.10, 0.05, 0.08, 0.12])

# Volatilities (annual)
volatilities = np.array([0.20, 0.08, 0.15, 0.25])

# Create a correlation matrix
# This is a symmetric matrix with 1s on the diagonal
correlation_matrix = np.array([
    [1.00, 0.20, 0.50, 0.30],
    [0.20, 1.00, 0.30, 0.10],
    [0.50, 0.30, 1.00, 0.40],
    [0.30, 0.10, 0.40, 1.00]
])

# Calculate the covariance matrix
covariance_matrix = np.zeros((n_assets, n_assets))
for i in range(n_assets):
    for j in range(n_assets):
        covariance_matrix[i, j] = correlation_matrix[i, j] * volatilities[i] * volatilities[j]

# Visualize the correlation matrix
plt.figure(figsize=(10, 8))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', xticklabels=assets, yticklabels=assets)
plt.title('Asset Correlation Matrix')
plt.tight_layout()
plt.show()

# Risk-free rate
risk_free_rate = 0.02

# Function to calculate portfolio return and volatility
def portfolio_performance(weights):
    # Returns
    returns = np.sum(weights * expected_returns)
    
    # Volatility
    volatility = np.sqrt(np.dot(weights.T, np.dot(covariance_matrix, weights)))
    
    return returns, volatility

# Function to minimize (negative Sharpe ratio)
def negative_sharpe_ratio(weights):
    returns, volatility = portfolio_performance(weights)
    sharpe_ratio = (returns - risk_free_rate) / volatility
    return -sharpe_ratio

# Constraints: weights sum to 1
constraints = {'type': 'eq', 'fun': lambda x: np.sum(x) - 1}

# Bounds: weights between 0 and 1 (no short selling)
bounds = tuple((0, 1) for _ in range(n_assets))

# Initial guess: equal weights
initial_weights = np.array([1/n_assets] * n_assets)

# Optimize
result = minimize(negative_sharpe_ratio, initial_weights, method='SLSQP', 
                 bounds=bounds, constraints=constraints)

# Extract the optimal weights
optimal_weights = result['x']

# Calculate the performance of the optimal portfolio
optimal_returns, optimal_volatility = portfolio_performance(optimal_weights)
optimal_sharpe = (optimal_returns - risk_free_rate) / optimal_volatility

print("Optimal Portfolio Weights:")
for asset, weight in zip(assets, optimal_weights):
    print(f"{asset}: {weight:.4f} ({weight*100:.2f}%)")
    
print(f"\nOptimal Portfolio Performance:")
print(f"Expected Annual Return: {optimal_returns:.4f} ({optimal_returns*100:.2f}%)")
print(f"Annual Volatility: {optimal_volatility:.4f} ({optimal_volatility*100:.2f}%)")
print(f"Sharpe Ratio: {optimal_sharpe:.4f}")

# Calculate total investment amounts (assuming $10,000 total)
total_investment = 10000
investment_amounts = optimal_weights * total_investment

print(f"\nCapital Allocation for ${total_investment}:")
for asset, amount in zip(assets, investment_amounts):
    print(f"{asset}: ${amount:.2f}")

# Monte Carlo simulation to generate random portfolios
num_portfolios = 10000
results = np.zeros((3, num_portfolios))
all_weights = np.zeros((num_portfolios, n_assets))

for i in range(num_portfolios):
    # Generate random weights
    weights = np.random.random(n_assets)
    weights = weights / np.sum(weights)
    all_weights[i, :] = weights
    
    # Calculate portfolio return and volatility
    portfolio_return, portfolio_volatility = portfolio_performance(weights)
    
    # Store results
    results[0, i] = portfolio_return
    results[1, i] = portfolio_volatility
    results[2, i] = (portfolio_return - risk_free_rate) / portfolio_volatility

# Plot the efficient frontier
plt.figure(figsize=(12, 8))
plt.scatter(results[1, :], results[0, :], c=results[2, :], cmap='viridis', 
            s=10, alpha=0.3, marker='o')

# Plot the optimal portfolio
plt.scatter(optimal_volatility, optimal_returns, c='red', s=100, marker='*', 
            label=f'Optimal Portfolio (Sharpe: {optimal_sharpe:.4f})')

# Add labels and legend
plt.colorbar(label='Sharpe Ratio')
plt.xlabel('Volatility (Standard Deviation)')
plt.ylabel('Expected Return')
plt.title('Portfolio Optimization: Risk vs Return')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Plot the optimal portfolio allocation
plt.figure(figsize=(12, 6))
colors = get_cmap('tab10').colors[:n_assets]

# Bar chart
plt.subplot(1, 2, 1)
plt.bar(assets, optimal_weights, color=colors)
plt.title('Optimal Portfolio Weights')
plt.ylabel('Weight')
plt.xticks(rotation=45)

# Pie chart
plt.subplot(1, 2, 2)
plt.pie(optimal_weights, labels=assets, autopct='%1.1f%%', colors=colors)
plt.title('Optimal Portfolio Allocation')

plt.tight_layout()
plt.show()

# Test some alternative allocations
alternative_allocations = [
    {'name': 'Equal Weight', 'weights': np.array([0.25, 0.25, 0.25, 0.25])},
    {'name': 'Stocks Heavy', 'weights': np.array([0.70, 0.10, 0.10, 0.10])},
    {'name': 'Bonds Heavy', 'weights': np.array([0.10, 0.70, 0.10, 0.10])},
    {'name': 'Optimal', 'weights': optimal_weights}
]

# Calculate performance metrics for each allocation
performance_data = []
for alloc in alternative_allocations:
    returns, volatility = portfolio_performance(alloc['weights'])
    sharpe = (returns - risk_free_rate) / volatility
    performance_data.append({
        'Allocation': alloc['name'],
        'Return': returns * 100,  # Convert to percentage
        'Volatility': volatility * 100,  # Convert to percentage
        'Sharpe Ratio': sharpe
    })

# Convert to DataFrame for easier visualization
performance_df = pd.DataFrame(performance_data)

# Create a comparison chart
fig, axes = plt.subplots(1, 3, figsize=(18, 6))

# Return comparison
axes[0].bar(performance_df['Allocation'], performance_df['Return'], color='skyblue')
axes[0].set_title('Expected Annual Return (%)')
axes[0].set_ylabel('Return (%)')
axes[0].set_ylim(0, max(performance_df['Return']) * 1.2)
axes[0].grid(axis='y', alpha=0.3)
axes[0].set_xticklabels(performance_df['Allocation'], rotation=45)

# Volatility comparison
axes[1].bar(performance_df['Allocation'], performance_df['Volatility'], color='salmon')
axes[1].set_title('Annual Volatility (%)')
axes[1].set_ylabel('Volatility (%)')
axes[1].set_ylim(0, max(performance_df['Volatility']) * 1.2)
axes[1].grid(axis='y', alpha=0.3)
axes[1].set_xticklabels(performance_df['Allocation'], rotation=45)

# Sharpe Ratio comparison
axes[2].bar(performance_df['Allocation'], performance_df['Sharpe Ratio'], color='lightgreen')
axes[2].set_title('Sharpe Ratio')
axes[2].set_ylabel('Sharpe Ratio')
axes[2].set_ylim(0, max(performance_df['Sharpe Ratio']) * 1.2)
axes[2].grid(axis='y', alpha=0.3)
axes[2].set_xticklabels(performance_df['Allocation'], rotation=45)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linear_sum_assignment

# Define the cost matrix
cost_matrix = np.array([
    [9, 2, 7],
    [6, 4, 3],
    [5, 8, 1]
])

# Solve the assignment problem
row_ind, col_ind = linear_sum_assignment(cost_matrix)

# Output the assignments and total cost
total_cost = cost_matrix[row_ind, col_ind].sum()

# Print results
print("Optimal assignment (Machine -> Task):")
for r, c in zip(row_ind, col_ind):
    print(f"Machine {r} -> Task {c} (Cost: {cost_matrix[r][c]})")
print(f"\nTotal Minimum Cost: {total_cost}")

import seaborn as sns

# Create a heatmap to show costs
plt.figure(figsize=(6, 5))
sns.heatmap(cost_matrix, annot=True, fmt="d", cmap="YlGnBu", cbar=False)

# Highlight the chosen assignments
for i, j in zip(row_ind, col_ind):
    plt.text(j + 0.2, i + 0.4, '✓', fontsize=18, color='red')

plt.title("Machine Assignment Cost Matrix")
plt.xlabel("Tasks")
plt.ylabel("Machines")
plt.xticks(np.arange(3) + 0.5, [f"Task {i}" for i in range(3)])
plt.yticks(np.arange(3) + 0.5, [f"Machine {i}" for i in range(3)])
plt.grid(False)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import random
from itertools import permutations
from math import sqrt

class VRPSolver:
    def __init__(self, depot, customers, demands, vehicle_capacity, num_vehicles):
        """
        Initialize the VRP solver with problem parameters
        
        Parameters:
        depot: tuple (x, y) - coordinates of the depot
        customers: list of tuples [(x1, y1), (x2, y2), ...] - coordinates of customers
        demands: list [d1, d2, ...] - demand of each customer
        vehicle_capacity: int - capacity of each vehicle
        num_vehicles: int - number of available vehicles
        """
        self.depot = depot
        self.customers = customers
        self.demands = demands
        self.vehicle_capacity = vehicle_capacity
        self.num_vehicles = num_vehicles
        self.distance_matrix = self._compute_distance_matrix()
        
    def _compute_distance_matrix(self):
        """Compute the distance matrix between all locations (depot + customers)"""
        all_locations = [self.depot] + self.customers
        n = len(all_locations)
        dist_matrix = np.zeros((n, n))
        
        for i in range(n):
            for j in range(n):
                if i != j:
                    x1, y1 = all_locations[i]
                    x2, y2 = all_locations[j]
                    dist_matrix[i, j] = sqrt((x2 - x1)**2 + (y2 - y1)**2)
                    
        return dist_matrix
    
    def _calculate_route_cost(self, route):
        """Calculate the total distance of a route"""
        total_distance = 0
        prev_location = 0  # Start at depot (index 0)
        
        for customer in route:
            total_distance += self.distance_matrix[prev_location, customer]
            prev_location = customer
            
        # Return to depot
        total_distance += self.distance_matrix[prev_location, 0]
        
        return total_distance
    
    def _is_route_feasible(self, route):
        """Check if a route satisfies capacity constraints"""
        total_demand = sum(self.demands[customer-1] for customer in route)
        return total_demand <= self.vehicle_capacity
    
    def solve_greedy(self):
        """Solve VRP using a greedy algorithm"""
        unassigned = list(range(1, len(self.customers) + 1))
        routes = [[] for _ in range(self.num_vehicles)]
        
        # Assign customers to vehicles until all are assigned
        for vehicle_idx in range(self.num_vehicles):
            current_location = 0  # Depot
            remaining_capacity = self.vehicle_capacity
            
            while unassigned:
                best_next = None
                best_distance = float('inf')
                
                # Find the nearest unassigned customer that fits in the vehicle
                for customer in unassigned:
                    if self.demands[customer-1] <= remaining_capacity:
                        distance = self.distance_matrix[current_location, customer]
                        if distance < best_distance:
                            best_distance = distance
                            best_next = customer
                
                if best_next is None:
                    # No more customers can fit in this vehicle
                    break
                    
                # Assign the customer to this vehicle
                routes[vehicle_idx].append(best_next)
                remaining_capacity -= self.demands[best_next-1]
                current_location = best_next
                unassigned.remove(best_next)
                
            if not unassigned:
                break
        
        return routes
    
    def solve_heuristic(self):
        """Solve VRP using a simple heuristic with local search"""
        # Start with a greedy solution
        routes = self.solve_greedy()
        
        # Try to improve each route using 2-opt local search
        for i in range(len(routes)):
            if len(routes[i]) >= 2:
                routes[i] = self._improve_route_2opt(routes[i])
        
        return routes
    
    def _improve_route_2opt(self, route):
        """Improve a route using 2-opt local search"""
        improved = True
        best_route = route.copy()
        best_cost = self._calculate_route_cost(best_route)
        
        while improved:
            improved = False
            
            for i in range(len(route)):
                for j in range(i + 2, len(route)):
                    # Skip if indices would make an invalid swap
                    if j - i == 1:
                        continue
                        
                    # Create a new route with a 2-opt swap
                    new_route = route[:i + 1] + route[i + 1:j + 1][::-1] + route[j + 1:]
                    new_cost = self._calculate_route_cost(new_route)
                    
                    if new_cost < best_cost:
                        best_route = new_route
                        best_cost = new_cost
                        improved = True
                        
            route = best_route.copy()
        
        return best_route
    
    def visualize_solution(self, routes):
        """Visualize the VRP solution"""
        plt.figure(figsize=(10, 8))
        
        # Plot depot
        plt.scatter([self.depot[0]], [self.depot[1]], c='red', s=100, marker='s', label='Depot')
        
        # Plot customers
        customer_x = [c[0] for c in self.customers]
        customer_y = [c[1] for c in self.customers]
        plt.scatter(customer_x, customer_y, c='blue', s=50, label='Customers')
        
        # Add customer numbers
        for i, (x, y) in enumerate(self.customers):
            plt.annotate(f"{i+1}\n(d:{self.demands[i]})", (x, y), xytext=(5, 5), 
                         textcoords='offset points')
        
        # Plot routes with different colors
        colors = ['green', 'purple', 'orange', 'brown', 'pink', 'gray', 'olive', 'cyan']
        for i, route in enumerate(routes):
            if not route:
                continue
                
            # Get coordinates for the route (including depot at start and end)
            route_x = [self.depot[0]]
            route_y = [self.depot[1]]
            
            for customer in route:
                route_x.append(self.customers[customer-1][0])
                route_y.append(self.customers[customer-1][1])
                
            # Return to depot
            route_x.append(self.depot[0])
            route_y.append(self.depot[1])
            
            plt.plot(route_x, route_y, c=colors[i % len(colors)], label=f'Vehicle {i+1}')
        
        plt.title("Vehicle Routing Problem Solution")
        plt.xlabel("X Coordinate")
        plt.ylabel("Y Coordinate")
        plt.legend()
        plt.grid(True)
        plt.axis('equal')
        plt.tight_layout()
        
        return plt

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

# Problem setup
depot = (0, 0)
num_customers = 10
customer_locations = [(random.uniform(-10, 10), random.uniform(-10, 10)) for _ in range(num_customers)]
customer_demands = [random.randint(1, 5) for _ in range(num_customers)]
vehicle_capacity = 10
num_vehicles = 3

# Create and solve the VRP
vrp = VRPSolver(depot, customer_locations, customer_demands, vehicle_capacity, num_vehicles)
solution = vrp.solve_heuristic()

# Print the solution
print("Solution:")
total_distance = 0
for i, route in enumerate(solution):
    route_distance = vrp._calculate_route_cost(route)
    total_distance += route_distance
    route_demand = sum(vrp.demands[customer-1] for customer in route)
    print(f"Vehicle {i+1}: {route} (Distance: {route_distance:.2f}, Load: {route_demand}/{vehicle_capacity})")
print(f"Total distance: {total_distance:.2f}")

# Visualize the solution
plt = vrp.visualize_solution(solution)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx

# Define the adjacency matrix
cost_matrix = np.array([
    [0, 3, 0, 0, 6, 5],
    [3, 0, 1, 0, 0, 4],
    [0, 1, 0, 6, 0, 4],
    [0, 0, 6, 0, 8, 5],
    [6, 0, 0, 8, 0, 2],
    [5, 4, 4, 5, 2, 0]
])

labels = ['A', 'B', 'C', 'D', 'E', 'F']
num_nodes = len(cost_matrix)

# Prim's Algorithm
def prim_mst(cost_matrix):
    selected_nodes = [False] * num_nodes
    selected_nodes[0] = True
    edges = []

    for _ in range(num_nodes - 1):
        min_cost = float('inf')
        a = b = -1
        for i in range(num_nodes):
            if selected_nodes[i]:
                for j in range(num_nodes):
                    if not selected_nodes[j] and cost_matrix[i][j]:
                        if cost_matrix[i][j] < min_cost:
                            min_cost = cost_matrix[i][j]
                            a, b = i, j
        if a != -1 and b != -1:
            edges.append((a, b, cost_matrix[a][b]))
            selected_nodes[b] = True

    return edges

# Compute MST
mst_edges = prim_mst(cost_matrix)

# Create graph and MST subgraph
G = nx.Graph()
MST = nx.Graph()

# Add all edges for base graph
for i in range(num_nodes):
    for j in range(i + 1, num_nodes):
        if cost_matrix[i][j] != 0:
            G.add_edge(labels[i], labels[j], weight=cost_matrix[i][j])

# Add only MST edges
for i, j, w in mst_edges:
    MST.add_edge(labels[i], labels[j], weight=w)

# Position nodes in a circular layout
pos = nx.circular_layout(G)

# Plot full graph with light edges
plt.figure(figsize=(10, 6))
nx.draw_networkx(G, pos, with_labels=True, node_color='lightblue', edge_color='gray', node_size=1200)
nx.draw_networkx_edges(MST, pos, edge_color='red', width=3)
edge_labels = nx.get_edge_attributes(MST, 'weight')
nx.draw_networkx_edge_labels(MST, pos, edge_labels=edge_labels)

plt.title("Minimum Spanning Tree (in Red) over Full Network")
plt.axis('off')
plt.show()

for u, v, w in mst_edges:
    print(f"{labels[u]} - {labels[v]} : {w}")

# Install PuLP if needed
!pip install pulp

import pulp
import itertools
import matplotlib.pyplot as plt

# Parameters
stock_length = 100
demand = {20: 48, 45: 35, 50: 24}
sizes = list(demand.keys())

# Generate all valid cutting patterns
def generate_patterns(sizes, stock_length):
    patterns = []
    for counts in itertools.product(*(range(stock_length // s + 1) for s in sizes)):
        total = sum(s * c for s, c in zip(sizes, counts))
        if 0 < total <= stock_length:
            patterns.append(counts)
    return patterns

patterns = generate_patterns(sizes, stock_length)

# Define the optimization problem
prob = pulp.LpProblem("CuttingStockProblem", pulp.LpMinimize)

# Define variables: how many times to use each pattern
x = [pulp.LpVariable(f"x_{i}", lowBound=0, cat='Integer') for i in range(len(patterns))]

# Objective: minimize number of stock rods used
prob += pulp.lpSum(x)

# Constraints: satisfy demand for each piece size
for i, size in enumerate(sizes):
    prob += pulp.lpSum(pattern[i] * x[j] for j, pattern in enumerate(patterns)) >= demand[size]

# Solve
prob.solve()

# Collect and display the solution
used_patterns = [(j, x[j].varValue, patterns[j]) for j in range(len(patterns)) if x[j].varValue > 0]
used_patterns = sorted(used_patterns, key=lambda tup: -tup[1])

print("Optimal number of stock rods used:", pulp.value(prob.objective))
print("\nCutting patterns used:")
for j, count, pattern in used_patterns:
    print(f"Use pattern {pattern} -> {int(count)} times")

# Visualization
colors = ['#ff9999','#66b3ff','#99ff99']
labels = [f"{s} unit" for s in sizes]

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

for i, (j, count, pattern) in enumerate(used_patterns):
    for k in range(int(count)):
        y = i + k * 0.1
        start = 0
        for size_idx, num in enumerate(pattern):
            for _ in range(num):
                ax.barh(y, sizes[size_idx], left=start, height=0.08, color=colors[size_idx])
                start += sizes[size_idx]

ax.set_title("Cutting Stock Solution Visualization")
ax.set_xlabel("Rod Length (100 units)")
ax.set_ylabel("Pattern Instances")
ax.set_xlim(0, stock_length)
ax.legend(labels, loc='upper right')
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def calculate_eoq(D, S, H):
    return np.sqrt((2 * D * S) / H)

def total_cost(D, S, H, Q):
    ordering_cost = (D / Q) * S
    holding_cost = (Q / 2) * H
    return ordering_cost + holding_cost

# Given data
D = 10000  # Annual demand
S = 50     # Ordering cost per order
H = 2      # Holding cost per unit per year

# Compute EOQ
eoq = calculate_eoq(D, S, H)
cost_at_eoq = total_cost(D, S, H, eoq)

print(f"Economic Order Quantity (EOQ): {eoq:.2f} units")
print(f"Total cost at EOQ: ${cost_at_eoq:.2f}")

# Visualizing the cost function
Q_values = np.linspace(50, 1000, 100)  # Order quantity range
costs = [total_cost(D, S, H, Q) for Q in Q_values]

plt.figure(figsize=(8, 5))
plt.plot(Q_values, costs, label='Total Cost', color='b')
plt.axvline(eoq, color='r', linestyle='--', label=f'EOQ: {eoq:.2f}')
plt.xlabel("Order Quantity (Q)")
plt.ylabel("Total Cost ($)")
plt.title("Total Inventory Cost vs Order Quantity")
plt.legend()
plt.grid()
plt.show()

import numpy as np
import pulp as pl
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

# Define the problem
def solve_staff_allocation(shift_requirements, shift_costs):
    """
    Solves the staff allocation problem
    
    Parameters:
    - shift_requirements: List of minimum staff required for each shift
    - shift_costs: List of cost per employee for each shift
    
    Returns:
    - model: The solved PuLP model
    - solution: Dictionary with solution details
    """
    # Create shift names
    shift_names = [f"Shift_{i+1}" for i in range(len(shift_requirements))]
    
    # Create the model
    model = pl.LpProblem("Staff_Allocation_Problem", pl.LpMinimize)
    
    # Define variables: number of staff for each shift (must be integers)
    staff_vars = {
        shift: pl.LpVariable(f"staff_{shift}", lowBound=0, cat=pl.LpInteger)
        for shift in shift_names
    }
    
    # Objective function: minimize total cost
    model += pl.lpSum([shift_costs[i] * staff_vars[shift] 
                      for i, shift in enumerate(shift_names)])
    
    # Constraints: meet minimum staffing requirements
    for i, shift in enumerate(shift_names):
        model += staff_vars[shift] >= shift_requirements[i], f"Min_Staff_{shift}"
    
    # Solve the model
    model.solve(pl.PULP_CBC_CMD(msg=False))
    
    # Extract solution details
    solution = {
        "status": pl.LpStatus[model.status],
        "total_cost": pl.value(model.objective),
        "staff_allocation": {shift: pl.value(var) for shift, var in staff_vars.items()}
    }
    
    return model, solution

# Example data
shift_requirements = [10, 8, 5]  # Minimum staff for morning, afternoon, night shifts
shift_costs = [100, 120, 140]    # Cost per employee for each shift

# Solve the problem
model, solution = solve_staff_allocation(shift_requirements, shift_costs)

# Print results
print("Solution Status:", solution["status"])
print("Total Cost:", solution["total_cost"])
print("Staff Allocation:")
for shift, staff in solution["staff_allocation"].items():
    print(f"  {shift}: {int(staff)} employees")

# Visualization
def visualize_solution(solution, shift_requirements, shift_costs):
    """
    Creates visualizations for the staff allocation solution
    """
    # Prepare data for visualization
    shifts = list(solution["staff_allocation"].keys())
    allocated_staff = [solution["staff_allocation"][shift] for shift in shifts]
    min_required = shift_requirements
    shift_labels = ["Morning (8am-4pm)", "Afternoon (4pm-12am)", "Night (12am-8am)"]
    
    # Create DataFrame for easier plotting
    df = pd.DataFrame({
        'Shift': shift_labels,
        'Allocated': allocated_staff,
        'Required': min_required,
        'Cost per Employee': shift_costs,
        'Total Cost': [allocated_staff[i] * shift_costs[i] for i in range(len(shifts))]
    })
    
    # Set up figure with multiple subplots
    fig, axs = plt.subplots(2, 1, figsize=(12, 12))
    
    # Staff allocation vs requirements
    ax1 = axs[0]
    bar_width = 0.35
    x = np.arange(len(shifts))
    
    bars1 = ax1.bar(x - bar_width/2, allocated_staff, bar_width, label='Allocated Staff', color='skyblue')
    bars2 = ax1.bar(x + bar_width/2, min_required, bar_width, label='Minimum Required', color='lightgreen')
    
    ax1.set_xlabel('Shift')
    ax1.set_ylabel('Number of Employees')
    ax1.set_title('Staff Allocation vs. Minimum Requirements')
    ax1.set_xticks(x)
    ax1.set_xticklabels(shift_labels)
    ax1.legend()
    
    # Add value labels on top of bars
    for bars in [bars1, bars2]:
        for bar in bars:
            height = bar.get_height()
            ax1.annotate(f'{int(height)}',
                        xy=(bar.get_x() + bar.get_width() / 2, height),
                        xytext=(0, 3),
                        textcoords="offset points",
                        ha='center', va='bottom')
    
    # Cost breakdown
    ax2 = axs[1]
    cost_data = df[['Shift', 'Total Cost']]
    sns.barplot(x='Shift', y='Total Cost', data=cost_data, ax=ax2, color='salmon')
    ax2.set_title('Cost Breakdown by Shift')
    ax2.set_ylabel('Total Cost ($)')
    
    # Add value labels on top of bars
    for i, p in enumerate(ax2.patches):
        height = p.get_height()
        ax2.annotate(f'${int(height)}',
                    xy=(p.get_x() + p.get_width() / 2, height),
                    xytext=(0, 3),
                    textcoords="offset points",
                    ha='center', va='bottom')
    
    plt.tight_layout()
    return fig

# Generate and display visualization
fig = visualize_solution(solution, shift_requirements, shift_costs)
plt.show()

# Let's try another scenario with different requirements and costs
print("\n--- Alternative Scenario ---")

# Different scenario with varying call volumes throughout the day
alt_requirements = [15, 12, 4]  # Higher morning/afternoon, lower night
alt_costs = [100, 110, 150]     # More expensive night shift

alt_model, alt_solution = solve_staff_allocation(alt_requirements, alt_costs)

# Print results for alternative scenario
print("Solution Status:", alt_solution["status"])
print("Total Cost:", alt_solution["total_cost"])
print("Staff Allocation:")
for shift, staff in alt_solution["staff_allocation"].items():
    print(f"  {shift}: {int(staff)} employees")

# Visualize alternative scenario
alt_fig = visualize_solution(alt_solution, alt_requirements, alt_costs)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
from matplotlib.colors import to_rgba

def greedy_set_cover(universe, sets, costs=None):
    """
    Greedy algorithm for the set cover problem.
    
    Args:
        universe: Set of elements to cover
        sets: Dictionary mapping set_id to the set of elements
        costs: Dictionary mapping set_id to the cost (default: all 1)
    
    Returns:
        List of set_ids that form the cover and the total cost
    """
    if costs is None:
        costs = {set_id: 1 for set_id in sets}
    
    # Make a copy of the universe and initialize the cover
    universe_copy = set(universe)
    cover = []
    total_cost = 0
    
    # While there are uncovered elements
    while universe_copy:
        # Find the set with the best cost-effectiveness
        best_set = None
        best_cost_effectiveness = 0
        
        for set_id, elements in sets.items():
            # Skip if already in the cover
            if set_id in cover:
                continue
            
            # Calculate new elements covered
            new_elements = len(universe_copy & elements)
            if new_elements == 0:
                continue
            
            # Calculate cost effectiveness (elements covered per unit cost)
            cost_effectiveness = new_elements / costs[set_id]
            
            if best_set is None or cost_effectiveness > best_cost_effectiveness:
                best_set = set_id
                best_cost_effectiveness = cost_effectiveness
        
        # If no set covers any remaining elements, break
        if best_set is None:
            break
        
        # Add the best set to the cover
        cover.append(best_set)
        total_cost += costs[best_set]
        
        # Remove covered elements
        universe_copy -= sets[best_set]
    
    return cover, total_cost

def visualize_set_cover(universe, sets, cover):
    """
    Visualize the set cover problem and solution using a bipartite graph.
    
    Args:
        universe: Set of elements
        sets: Dictionary mapping set_id to the set of elements
        cover: List of set_ids that form the cover
    """
    G = nx.Graph()
    
    # Add nodes
    for element in universe:
        G.add_node(element, bipartite=0)  # 0 for elements
    
    for set_id in sets:
        G.add_node(set_id, bipartite=1)  # 1 for sets
    
    # Add edges
    for set_id, elements in sets.items():
        for element in elements:
            G.add_edge(set_id, element)
    
    # Create layout
    pos = nx.bipartite_layout(G, [node for node, attr in G.nodes(data=True) if attr['bipartite'] == 0])
    
    # Set up colors
    cover_color = 'red'
    non_cover_color = 'lightblue'
    element_color = 'lightgreen'
    
    # Create node colors
    node_colors = []
    for node in G.nodes():
        if node in cover:
            node_colors.append(cover_color)
        elif node in universe:
            node_colors.append(element_color)
        else:
            node_colors.append(non_cover_color)
    
    # Create edge colors
    edge_colors = []
    for u, v in G.edges():
        if u in cover or v in cover:
            edge_colors.append(cover_color)
        else:
            edge_colors.append('gray')
    
    # Draw the graph
    plt.figure(figsize=(12, 8))
    nx.draw(G, pos, 
            with_labels=True, 
            node_color=node_colors, 
            edge_color=edge_colors, 
            node_size=700,
            font_weight='bold',
            width=1.5)
    
    # Create a custom legend
    from matplotlib.patches import Patch
    legend_elements = [
        Patch(facecolor=cover_color, edgecolor='black', label='Sets in Cover'),
        Patch(facecolor=non_cover_color, edgecolor='black', label='Sets not in Cover'),
        Patch(facecolor=element_color, edgecolor='black', label='Elements to Cover')
    ]
    plt.legend(handles=legend_elements, loc='upper right')
    
    plt.title('Set Cover Visualization')
    plt.tight_layout()
    plt.show()

def visualize_set_coverage_steps(universe, sets, cover):
    """
    Visualize the progressive coverage of elements by the sets in the cover.
    
    Args:
        universe: Set of elements
        sets: Dictionary mapping set_id to the set of elements
        cover: List of set_ids that form the cover (ordered by selection)
    """
    # Create a new figure
    plt.figure(figsize=(12, 6))
    
    # Convert universe to a list for consistent ordering
    universe_list = list(universe)
    universe_list.sort()  # Sort for better visualization
    
    # Create a matrix to show coverage
    # Rows are steps (adding each set)
    # Columns are elements
    coverage_matrix = np.zeros((len(cover) + 1, len(universe_list)))
    
    # First row shows no coverage
    covered_elements = set()
    
    # For each step (adding a set to the cover)
    for i, set_id in enumerate(cover):
        # Add newly covered elements
        covered_elements.update(sets[set_id])
        # Mark covered elements in the matrix
        for j, element in enumerate(universe_list):
            if element in covered_elements:
                coverage_matrix[i+1, j] = 1
    
    # Plot as a heatmap
    plt.imshow(coverage_matrix, aspect='auto', cmap='Blues')
    
    # Add labels
    plt.yticks(np.arange(len(cover) + 1), 
               ['Initial'] + [f'Add {set_id}' for set_id in cover])
    plt.xticks(np.arange(len(universe_list)), universe_list, rotation=45)
    
    # Add grid
    plt.grid(False)
    for i in range(len(cover) + 1):
        plt.axhline(i - 0.5, color='gray', linestyle='-', linewidth=0.5)
    for j in range(len(universe_list)):
        plt.axvline(j - 0.5, color='gray', linestyle='-', linewidth=0.5)
    
    # Add coverage percentage for each step
    for i in range(len(cover) + 1):
        coverage_pct = np.sum(coverage_matrix[i]) / len(universe_list) * 100
        plt.text(len(universe_list) - 0.5, i, f'{coverage_pct:.1f}%', 
                 verticalalignment='center')
    
    plt.colorbar(label='Covered (1) / Not Covered (0)')
    plt.title('Progressive Element Coverage')
    plt.tight_layout()
    plt.show()

# Define our problem
universe = {'City1', 'City2', 'City3', 'City4', 'City5', 
            'City6', 'City7', 'City8', 'City9', 'City10'}

sets = {
    'Station1': {'City1', 'City2', 'City3'},
    'Station2': {'City2', 'City4', 'City5', 'City6'},
    'Station3': {'City3', 'City6', 'City7'},
    'Station4': {'City4', 'City8', 'City9'},
    'Station5': {'City5', 'City7', 'City9', 'City10'},
    'Station6': {'City1', 'City10'}
}

# Solve using the greedy algorithm
cover, total_cost = greedy_set_cover(universe, sets)
print(f"Optimal cover: {cover}")
print(f"Total cost: {total_cost}")
print(f"Number of sets used: {len(cover)}")

# Verify the cover is valid
covered = set()
for set_id in cover:
    covered.update(sets[set_id])

print(f"All elements covered: {covered == universe}")

# Visualize the solution
visualize_set_cover(universe, sets, cover)
visualize_set_coverage_steps(universe, sets, cover)

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import random
from matplotlib.colors import ListedColormap
from scipy.spatial import Voronoi

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

# Generate synthetic data for our state
num_precincts = 100
num_districts = 5

# Create precinct locations and populations
x_coords = np.random.rand(num_precincts) * 100
y_coords = np.random.rand(num_precincts) * 100
positions = np.column_stack((x_coords, y_coords))

# Generate population data with some variation
base_population = 10000
population_variation = 3000
populations = np.random.randint(
    base_population - population_variation,
    base_population + population_variation,
    num_precincts
)

# Create a graph where precincts are nodes connected to nearby precincts
G = nx.Graph()

# Add nodes with attributes
for i in range(num_precincts):
    G.add_node(i, pos=(x_coords[i], y_coords[i]), population=populations[i])

# Connect precincts that are close to each other
distance_threshold = 15  # Adjust based on your coordinate scale
for i in range(num_precincts):
    for j in range(i+1, num_precincts):
        dist = np.sqrt((x_coords[i] - x_coords[j])**2 + (y_coords[i] - y_coords[j])**2)
        if dist < distance_threshold:
            G.add_edge(i, j)

# Check if the graph is connected
if not nx.is_connected(G):
    largest_cc = max(nx.connected_components(G), key=len)
    if len(largest_cc) < num_precincts:
        print(f"Warning: Generated graph is not connected. Using largest component with {len(largest_cc)} nodes.")
        G = G.subgraph(largest_cc).copy()
        # Adjust our precinct count
        num_precincts = len(G.nodes)

# Helper function to check if a district is contiguous
def is_contiguous(G, district_nodes):
    if not district_nodes:
        return True
    subgraph = G.subgraph(district_nodes)
    return nx.is_connected(subgraph)

# Helper function to calculate district population
def district_population(district, populations):
    return sum(populations[node] for node in district)

# Initialize districts using region growing approach
def initialize_districts(G, num_districts, populations):
    nodes = list(G.nodes())
    
    # Start with seed nodes far apart
    seed_indices = np.linspace(0, len(nodes)-1, num_districts, dtype=int)
    seeds = [nodes[i] for i in seed_indices]
    
    # Initialize districts with seed nodes
    districts = [set([seed]) for seed in seeds]
    unassigned = set(nodes) - set(seeds)
    
    # Grow districts by adding adjacent nodes one at a time
    while unassigned:
        for d_idx, district in enumerate(districts):
            if not unassigned:
                break
                
            # Find unassigned nodes adjacent to this district
            frontier = set()
            for node in district:
                frontier.update(n for n in G.neighbors(node) if n in unassigned)
            
            if frontier:
                # Add one node to the district
                new_node = random.choice(list(frontier))
                district.add(new_node)
                unassigned.remove(new_node)
    
    return districts

# Score function to evaluate a districting plan
def score_plan(districts, populations):
    district_pops = [district_population(d, populations) for d in districts]
    
    # Population balance score (lower is better)
    avg_pop = sum(district_pops) / len(district_pops)
    pop_deviation = max(abs(p - avg_pop) / avg_pop for p in district_pops)
    
    return pop_deviation

# Local search algorithm to improve the districting plan
def optimize_districts(G, initial_districts, populations, iterations=1000):
    current_districts = [set(d) for d in initial_districts]
    current_score = score_plan(current_districts, populations)
    
    best_districts = [set(d) for d in current_districts]
    best_score = current_score
    
    for i in range(iterations):
        # Select a random border precinct
        moved = False
        
        # Try to move a random node from one district to another
        for source_idx in random.sample(range(len(current_districts)), len(current_districts)):
            if moved:
                break
                
            source_district = current_districts[source_idx]
            
            # Find border nodes (nodes with neighbors in other districts)
            border_nodes = set()
            for node in source_district:
                for neighbor in G.neighbors(node):
                    if any(neighbor in d for d in current_districts if d != source_district):
                        border_nodes.add(node)
                        break
            
            if not border_nodes:
                continue
                
            # Try to move a random border node
            node_to_move = random.choice(list(border_nodes))
            
            # Find neighboring districts
            target_districts = []
            for d_idx, district in enumerate(current_districts):
                if d_idx != source_idx:
                    for neighbor in G.neighbors(node_to_move):
                        if neighbor in district:
                            target_districts.append(d_idx)
                            break
            
            if not target_districts:
                continue
                
            target_idx = random.choice(target_districts)
            
            # Check if removing the node keeps the source district contiguous
            temp_source = source_district.copy()
            temp_source.remove(node_to_move)
            
            if not is_contiguous(G, temp_source):
                continue
                
            # Make the move
            new_districts = [set(d) for d in current_districts]
            new_districts[source_idx].remove(node_to_move)
            new_districts[target_idx].add(node_to_move)
            
            # Evaluate the new plan
            new_score = score_plan(new_districts, populations)
            
            # Accept if it's better
            if new_score < current_score:
                current_districts = new_districts
                current_score = new_score
                moved = True
                
                # Update best if it's the best so far
                if new_score < best_score:
                    best_districts = [set(d) for d in new_districts]
                    best_score = new_score
                    print(f"Iteration {i}: New best score = {best_score:.4f}")
    
    return best_districts, best_score

# Run the algorithm
print("Initializing districts...")
initial_districts = initialize_districts(G, num_districts, populations)
print("Optimizing districts...")
optimized_districts, final_score = optimize_districts(G, initial_districts, populations)

# Calculate district populations
district_pops = [district_population(d, populations) for d in optimized_districts]
total_pop = sum(populations)
district_pop_percentages = [p/total_pop*100 for p in district_pops]

# Prepare data for visualization
node_colors = np.zeros(num_precincts)
for i, district in enumerate(optimized_districts):
    for node in district:
        node_colors[node] = i

# Get positions for drawing
pos = nx.get_node_attributes(G, 'pos')

# Set up a colorful visualization
plt.figure(figsize=(15, 10))

# Plot the graph with district colors
cmap = plt.cm.get_cmap('tab10', num_districts)
nx.draw(G, pos, node_color=node_colors, cmap=cmap, node_size=100, with_labels=False, width=0.5)

# Add a title and legend
plt.title('Optimized Redistricting Plan', fontsize=16)
patches = [plt.plot([],[], marker="o", ms=10, ls="", mec=None, color=cmap(i), 
                    label=f"District {i+1}: {district_pops[i]} people ({district_pop_percentages[i]:.1f}%)"
                   )[0] for i in range(num_districts)]
plt.legend(handles=patches, loc='upper right', fontsize=12)

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

# Create a population distribution bar chart
plt.figure(figsize=(12, 6))
bars = plt.bar(range(1, num_districts+1), district_pops, color=cmap(range(num_districts)))
plt.axhline(y=total_pop/num_districts, color='r', linestyle='--', label='Ideal Population')

# Add labels and annotations
plt.xlabel('District', fontsize=14)
plt.ylabel('Population', fontsize=14)
plt.title('Population Distribution by District', fontsize=16)
plt.xticks(range(1, num_districts+1))

# Add population values on top of bars
for bar in bars:
    height = bar.get_height()
    plt.text(bar.get_x() + bar.get_width()/2., height + 5000,
            f'{int(height)}',
            ha='center', va='bottom', fontsize=12)

plt.legend()
plt.grid(axis='y', alpha=0.3)
plt.savefig('district_populations.png', dpi=300, bbox_inches='tight')
plt.show()

# Calculate and display metrics
print("\nRedistricting Results:")
print(f"Number of districts: {num_districts}")
print(f"Total population: {total_pop}")
print(f"Ideal district population: {total_pop / num_districts:.1f}")
print("\nDistrict populations:")
for i, pop in enumerate(district_pops):
    deviation = (pop - total_pop/num_districts) / (total_pop/num_districts) * 100
    print(f"District {i+1}: {pop} people (deviation: {deviation:+.2f}%)")

print(f"\nPopulation deviation score: {final_score:.4f}")
print(f"Maximum population deviation: {max(abs(p - total_pop/num_districts) / (total_pop/num_districts) * 100 for p in district_pops):.2f}%")