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

# Define products
products = ['A', 'B', 'C', 'D']

# Transition probabilities matrix
# Rows: From product, Columns: To product
P = np.array([
    [0, 0.5, 0.3, 0.2],  # From A
    [0.4, 0, 0.4, 0.2],  # From B
    [0.3, 0.3, 0, 0.4],  # From C
    [0.2, 0.3, 0.5, 0]   # From D
])

# Purchase probabilities at each product location
S = np.array([0.6, 0.5, 0.4, 0.3])

# Generate all possible arrangements (permutations)
arrangements = list(itertools.permutations(products))

# Evaluate expected sales for each arrangement
results = []

for arrangement in arrangements:
    # Map product to index in new arrangement
    idx_map = {product: i for i, product in enumerate(arrangement)}
    
    expected_sales = 0
    for i, from_product in enumerate(products):
        for j, to_product in enumerate(products):
            if i != j:
                expected_sales += P[i, j] * S[idx_map[to_product]]
    
    results.append((arrangement, expected_sales))

# Find the best arrangement
best_arrangement, best_sales = max(results, key=lambda x: x[1])

print("Best Arrangement:", best_arrangement)
print("Expected Sales:", round(best_sales, 3))

# Plot all arrangements
arrangement_labels = ['-'.join(a) for a, _ in results]
sales_values = [s for _, s in results]

plt.figure(figsize=(12, 6))
plt.barh(arrangement_labels, sales_values, color='skyblue')
plt.xlabel('Expected Sales')
plt.title('Expected Sales for Different Store Layouts')
plt.gca().invert_yaxis()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import itertools
from matplotlib.patches import Circle
import pulp

# Set up the problem
np.random.seed(42)  # For reproducibility

# Generate demand points (randomly distributed in a 2D space)
n_points = 15
points = np.random.rand(n_points, 2) * 10  # Points in a 10x10 grid

# Generate potential facility locations
n_facilities = 10
facilities = np.random.rand(n_facilities, 2) * 10  # Facilities in the same 10x10 grid

# Define coverage radius
coverage_radius = 3.0

# Compute coverage matrix
# coverage[i, j] = 1 if facility j covers point i, 0 otherwise
coverage = np.zeros((n_points, n_facilities), dtype=int)
for i in range(n_points):
    for j in range(n_facilities):
        distance = np.sqrt(np.sum((points[i] - facilities[j])**2))
        if distance <= coverage_radius:
            coverage[i, j] = 1

# Create a PuLP model
model = pulp.LpProblem("Set_Cover_Problem", pulp.LpMinimize)

# Define decision variables: 1 if facility j is selected, 0 otherwise
x = [pulp.LpVariable(f"x_{j}", cat='Binary') for j in range(n_facilities)]

# Define objective function: minimize the number of facilities
model += pulp.lpSum(x)

# Add constraints: each point must be covered by at least one facility
for i in range(n_points):
    model += pulp.lpSum(coverage[i, j] * x[j] for j in range(n_facilities)) >= 1

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

# Extract the solution
selected_facilities = [j for j in range(n_facilities) if pulp.value(x[j]) > 0.5]
print(f"Selected facilities: {selected_facilities}")
print(f"Number of facilities: {len(selected_facilities)}")

# Check that all points are covered
covered = np.zeros(n_points, dtype=bool)
for j in selected_facilities:
    covered |= (coverage[:, j] == 1)
print(f"All points covered: {np.all(covered)}")

# Visualization
plt.figure(figsize=(10, 8))

# Plot demand points
plt.scatter(points[:, 0], points[:, 1], color='blue', label='Demand Points')

# Plot all facilities
plt.scatter(facilities[:, 0], facilities[:, 1], color='lightgray', 
            marker='s', s=100, label='Unselected Facilities')

# Overlay selected facilities
plt.scatter(facilities[selected_facilities, 0], facilities[selected_facilities, 1], 
            color='red', marker='s', s=100, label='Selected Facilities')

# Draw coverage circles for selected facilities
for j in selected_facilities:
    plt.gca().add_patch(Circle(facilities[j], coverage_radius, 
                              fill=False, color='red', alpha=0.3))

plt.grid(True)
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.title('Set Cover Problem Solution')
plt.legend()
plt.axis('equal')
plt.xlim(-1, 11)
plt.ylim(-1, 11)
plt.tight_layout()
plt.show()

# Bonus: Let's analyze how many points each selected facility covers
coverage_counts = []
for j in selected_facilities:
    count = np.sum(coverage[:, j])
    coverage_counts.append(count)
    print(f"Facility {j} covers {count} points")

# Visualize the coverage distribution
plt.figure(figsize=(8, 5))
plt.bar(range(len(selected_facilities)), coverage_counts, 
        tick_label=[f"Facility {j}" for j in selected_facilities])
plt.xlabel('Selected Facilities')
plt.ylabel('Number of Points Covered')
plt.title('Coverage Distribution Among Selected Facilities')
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

# Bonus: Let's check for redundancy in our solution
for j in selected_facilities:
    # Remove this facility temporarily
    temp_coverage = np.zeros(n_points, dtype=bool)
    for k in selected_facilities:
        if k != j:
            temp_coverage |= (coverage[:, k] == 1)
    
    if np.all(temp_coverage):
        print(f"Facility {j} is redundant - all points would still be covered without it")

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import factorial
import pandas as pd
import seaborn as sns

def erlang_c(s, lam, mu):
    """
    Calculate the probability of waiting based on Erlang C formula
    
    Parameters:
    s (int): Number of agents
    lam (float): Arrival rate (calls per minute)
    mu (float): Service rate (calls per agent per minute)
    
    Returns:
    float: Probability that a call will need to wait (P(W>0))
    """
    rho = lam / (s * mu)
    
    if rho >= 1:
        return 1.0  # System is unstable if rho >= 1
    
    # Calculate the probability of zero customers in the system (P0)
    sum_term = sum([(s * rho)**n / factorial(n) for n in range(s)])
    p0 = 1 / (sum_term + (s * rho)**s / (factorial(s) * (1 - rho)))
    
    # Calculate P(W>0) using the Erlang C formula
    pw = (s * rho)**s / (factorial(s) * (1 - rho)) * p0
    
    return pw

def wait_probability(s, lam, mu, t):
    """
    Calculate the probability that wait time exceeds t seconds
    
    Parameters:
    s (int): Number of agents
    lam (float): Arrival rate (calls per minute)
    mu (float): Service rate (calls per agent per minute)
    t (float): Time threshold in seconds
    
    Returns:
    float: Probability that a call will wait more than t seconds
    """
    # Convert t from seconds to minutes for consistent units
    t_minutes = t / 60
    
    rho = lam / (s * mu)
    
    if rho >= 1:
        return 1.0  # System is unstable
    
    pw = erlang_c(s, lam, mu)
    wait_prob = pw * np.exp(-s * mu * (1 - rho) * t_minutes)
    
    return wait_prob

def min_agents_for_service_level(lam, mu, t, service_level):
    """
    Find the minimum number of agents needed to meet the service level requirement
    
    Parameters:
    lam (float): Arrival rate (calls per minute)
    mu (float): Service rate (calls per agent per minute)
    t (float): Time threshold in seconds
    service_level (float): Required service level (e.g., 0.8 for 80%)
    
    Returns:
    int: Minimum number of agents needed
    """
    # Start with the minimum theoretical number of agents needed for stability
    s_min = int(np.ceil(lam / mu))
    
    # Increase the number of agents until the service level constraint is met
    s = s_min
    while True:
        # Calculate the probability that a customer waits less than t seconds
        prob_wait_less_than_t = 1 - wait_probability(s, lam, mu, t)
        
        if prob_wait_less_than_t >= service_level:
            return s
        
        s += 1

# Example parameters for a call center
arrival_rate = 5  # 5 calls per minute (300 calls per hour)
service_rate = 0.5  # 0.5 calls per agent per minute (30 calls per hour per agent)
time_threshold = 20  # 20 seconds
service_level_target = 0.8  # 80% of calls answered within the time threshold

# Find the minimum number of agents needed
min_agents = min_agents_for_service_level(arrival_rate, service_rate, time_threshold, service_level_target)

print(f"Minimum agents needed: {min_agents}")
print(f"Service level achieved: {(1 - wait_probability(min_agents, arrival_rate, service_rate, time_threshold)) * 100:.2f}%")

# Let's analyze the impact of varying the number of agents
agent_range = range(min_agents-3, min_agents+4)
service_levels = []

for s in agent_range:
    if s <= arrival_rate / service_rate:  # Skip unstable configurations
        service_levels.append(0)
    else:
        sl = (1 - wait_probability(s, arrival_rate, service_rate, time_threshold)) * 100
        service_levels.append(sl)

# Create a DataFrame for the results
results_df = pd.DataFrame({
    'Number of Agents': list(agent_range),
    'Service Level (%)': service_levels
})

# Plot the results
plt.figure(figsize=(10, 6))
sns.lineplot(data=results_df, x='Number of Agents', y='Service Level (%)', marker='o', linewidth=2)
plt.axhline(y=service_level_target*100, color='r', linestyle='--', label=f'Target ({service_level_target*100}%)')
plt.axvline(x=min_agents, color='g', linestyle='--', label=f'Minimum Agents ({min_agents})')
plt.title('Service Level vs. Number of Agents', fontsize=14)
plt.grid(True, alpha=0.3)
plt.legend()
plt.xticks(list(agent_range))
plt.ylim(0, 100)
plt.show()

# Let's also analyze how service level varies with different arrival rates
arrival_rates = np.linspace(3, 7, 9)  # From 3 to 7 calls per minute
agent_counts = []

for lam in arrival_rates:
    agents = min_agents_for_service_level(lam, service_rate, time_threshold, service_level_target)
    agent_counts.append(agents)

# Create a DataFrame for sensitivity analysis
sensitivity_df = pd.DataFrame({
    'Arrival Rate (calls/min)': arrival_rates,
    'Minimum Agents Required': agent_counts
})

# Plot sensitivity analysis
plt.figure(figsize=(10, 6))
sns.lineplot(data=sensitivity_df, x='Arrival Rate (calls/min)', y='Minimum Agents Required', marker='o', linewidth=2)
plt.title('Sensitivity Analysis: Impact of Call Volume on Staffing Requirements', fontsize=14)
plt.grid(True, alpha=0.3)
plt.xticks(arrival_rates)
plt.yticks(range(min(agent_counts)-1, max(agent_counts)+2))
plt.show()

# Let's create a heatmap showing how the service level varies with different combinations of agents and arrival rates
agent_range = range(8, 17)
arrival_range = [3, 4, 5, 6, 7]
heatmap_data = []

for agents in agent_range:
    row = []
    for lam in arrival_range:
        if agents <= lam / service_rate:  # Skip unstable configurations
            sl = 0
        else:
            sl = (1 - wait_probability(agents, lam, service_rate, time_threshold)) * 100
        row.append(sl)
    heatmap_data.append(row)

# Create a heatmap
plt.figure(figsize=(10, 8))
sns.heatmap(heatmap_data, annot=True, fmt=".1f", cmap="YlGnBu",
            xticklabels=arrival_range, yticklabels=agent_range,
            cbar_kws={'label': 'Service Level (%)'})
plt.xlabel('Arrival Rate (calls/min)')
plt.ylabel('Number of Agents')
plt.title('Service Level Heatmap', fontsize=14)
plt.tight_layout()
plt.show()

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()