import networkx as nx
import matplotlib.pyplot as plt

# Define the graph edges and their weights
edges = [
    ("A", "B", 7),
    ("A", "D", 5),
    ("B", "C", 8),
    ("B", "D", 9),
    ("B", "E", 7),
    ("C", "E", 5),
    ("D", "E", 15),
    ("D", "F", 6),
    ("E", "F", 8),
    ("E", "G", 9),
    ("F", "G", 11)
]

# Create an undirected graph
G = nx.Graph()
G.add_weighted_edges_from(edges)

# Use Kruskal's algorithm to compute the Minimum Spanning Tree
mst = nx.minimum_spanning_tree(G, algorithm="kruskal")

# Display MST edges and total cost
print("Edges in MST:")
total_cost = 0
for u, v, data in mst.edges(data=True):
    print(f"{u} - {v}: {data['weight']}")
    total_cost += data['weight']
print(f"\nTotal Cost of MST: {total_cost}")

# Position nodes using spring layout for clarity
pos = nx.spring_layout(G, seed=42)

# Draw original graph in light gray
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
nx.draw(G, pos, with_labels=True, node_color='lightblue', node_size=800, edge_color='gray')
nx.draw_networkx_edge_labels(G, pos, edge_labels={(u,v): d['weight'] for u,v,d in G.edges(data=True)})
plt.title("Original Graph")

# Draw MST in second subplot
plt.subplot(1, 2, 2)
nx.draw(G, pos, with_labels=True, node_color='lightblue', node_size=800, edge_color='lightgray')
nx.draw(mst, pos, with_labels=True, node_color='lightgreen', node_size=800, edge_color='green', width=2)
nx.draw_networkx_edge_labels(mst, pos, edge_labels={(u,v): d['weight'] for u,v,d in mst.edges(data=True)})
plt.title("Minimum Spanning Tree")

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

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

# Problem parameters
total_budget = 10  # Total budget in million dollars
periods = 5  # Number of investment periods (years)
discretization = 100  # How finely we discretize the budget
step_size = total_budget / discretization  # The step size for discretization

# Generate return functions for each period
# For this example, we'll use quadratic functions with diminishing returns
# f_t(x) = a_t * sqrt(x) - b_t * x
def create_return_functions():
    # Parameters for return functions (a_t, b_t)
    # These create unique return profiles for each period
    params = [
        (4.0, 0.2),  # Period 1: Higher initial returns but quicker diminishment
        (3.5, 0.1),  # Period 2: Moderate returns with slower diminishment
        (5.0, 0.3),  # Period 3: Highest initial returns but quick diminishment
        (3.0, 0.05), # Period 4: Lower returns but very slow diminishment
        (4.5, 0.15)  # Period 5: Good returns with moderate diminishment
    ]
    
    # Create return functions for each period
    return_functions = []
    for a, b in params:
        # Define the return function for this period
        def f(x, a=a, b=b):
            if x == 0:
                return 0
            return a * np.sqrt(x) - b * x
        return_functions.append(f)
    
    return return_functions

# Create the return functions
return_functions = create_return_functions()

# Function to plot the return functions
def plot_return_functions(return_functions):
    plt.figure(figsize=(10, 6))
    x = np.linspace(0, total_budget, 100)
    
    for t, f in enumerate(return_functions):
        returns = [f(xi) for xi in x]
        plt.plot(x, returns, label=f'Period {t+1}')
    
    plt.title('Return Functions for Each Period')
    plt.xlabel('Investment Amount ($ Million)')
    plt.ylabel('Return ($ Million)')
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

# Plot the return functions
plot_return_functions(return_functions)

# Dynamic Programming Solution
def solve_dp():
    # Initialize value table and decision table
    # Value table: V[t][b] = maximum return from periods t to T with budget b
    # Decision table: D[t][b] = optimal investment amount in period t with budget b
    budget_steps = int(total_budget / step_size) + 1
    V = np.zeros((periods + 1, budget_steps))
    D = np.zeros((periods, budget_steps))
    
    # DP backward recursion
    for t in range(periods - 1, -1, -1):
        for b in range(budget_steps):
            budget = b * step_size
            
            # Try different investment amounts for period t
            best_value = -float('inf')
            best_decision = 0
            
            for i in range(b + 1):
                investment = i * step_size
                remaining = budget - investment
                remaining_idx = int(remaining / step_size)
                
                current_return = return_functions[t](investment)
                future_return = V[t + 1][remaining_idx]
                total_return = current_return + future_return
                
                if total_return > best_value:
                    best_value = total_return
                    best_decision = investment
            
            V[t][b] = best_value
            D[t][b] = best_decision
    
    # Extract the optimal policy
    optimal_policy = []
    remaining_budget = total_budget
    
    for t in range(periods):
        b_idx = int(remaining_budget / step_size)
        investment = D[t][b_idx]
        optimal_policy.append(investment)
        remaining_budget -= investment
    
    return V, D, optimal_policy

# Solve the problem
V, D, optimal_policy = solve_dp()

# Print the results
print("Optimal Investment Policy:")
for t, investment in enumerate(optimal_policy):
    print(f"Period {t+1}: ${investment:.2f} million")

print("\nTotal Expected Return: ${:.2f} million".format(V[0][int(total_budget / step_size)]))

# Function to plot the optimal investment policy
def plot_optimal_policy(optimal_policy):
    plt.figure(figsize=(10, 6))
    periods_list = list(range(1, periods + 1))
    
    # Bar plot for investments
    bars = plt.bar(periods_list, optimal_policy, color='skyblue', alpha=0.7)
    
    # Calculate returns for each period
    returns = [return_functions[t](optimal_policy[t]) for t in range(periods)]
    
    # Line plot for returns
    plt.plot(periods_list, returns, 'ro-', label='Return', linewidth=2)
    
    # Adding value labels on top of the bars
    for i, (bar, investment, ret) in enumerate(zip(bars, optimal_policy, returns)):
        plt.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1, 
                 f"${investment:.2f}M", ha='center', va='bottom')
        plt.text(i+1, ret + 0.1, f"${ret:.2f}M", ha='center', va='bottom', color='red')
    
    plt.title('Optimal Investment Strategy and Returns by Period')
    plt.xlabel('Period')
    plt.ylabel('Amount ($ Million)')
    plt.legend(['Return', 'Investment'])
    plt.xticks(periods_list)
    plt.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.show()

# Function to visualize the value function
def plot_value_function(V):
    plt.figure(figsize=(10, 6))
    budget_values = np.linspace(0, total_budget, V.shape[1])
    
    for t in range(periods):
        plt.plot(budget_values, V[t], label=f'Period {t+1}')
    
    plt.title('Value Function by Period and Available Budget')
    plt.xlabel('Available Budget ($ Million)')
    plt.ylabel('Maximum Expected Return ($ Million)')
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

# Function to create a heatmap of the decision table
def plot_decision_heatmap(D):
    plt.figure(figsize=(12, 8))
    
    # We'll use a subset of the decision table for better visualization
    step = max(1, int(D.shape[1] / 20))
    subset_D = D[:, ::step]
    subset_budgets = np.linspace(0, total_budget, subset_D.shape[1])
    
    ax = sns.heatmap(subset_D, cmap='viridis', 
                    xticklabels=[f'{x:.1f}' for x in subset_budgets],
                    yticklabels=[f'Period {t+1}' for t in range(periods)])
    
    plt.title('Optimal Investment Decision by Period and Available Budget')
    plt.xlabel('Available Budget ($ Million)')
    plt.ylabel('Period')
    cbar = ax.collections[0].colorbar
    cbar.set_label('Investment Amount ($ Million)')
    plt.tight_layout()
    plt.show()

# Plot the optimal policy
plot_optimal_policy(optimal_policy)

# Plot the value function
plot_value_function(V)

# Plot the decision heatmap
plot_decision_heatmap(D)

# Sensitivity Analysis: Effect of total budget
def budget_sensitivity_analysis():
    budget_range = np.linspace(2, 20, 10)  # Test budgets from $2M to $20M
    total_returns = []
    
    for budget in budget_range:
        # Reconfigure the problem with a new budget
        discretization = 100
        step_size = budget / discretization
        
        # Solve the DP problem
        budget_steps = int(budget / step_size) + 1
        V = np.zeros((periods + 1, budget_steps))
        D = np.zeros((periods, budget_steps))
        
        # DP backward recursion
        for t in range(periods - 1, -1, -1):
            for b in range(budget_steps):
                b_amount = b * step_size
                
                best_value = -float('inf')
                best_decision = 0
                
                for i in range(b + 1):
                    investment = i * step_size
                    remaining = b_amount - investment
                    remaining_idx = int(remaining / step_size)
                    
                    current_return = return_functions[t](investment)
                    future_return = V[t + 1][remaining_idx]
                    total_return = current_return + future_return
                    
                    if total_return > best_value:
                        best_value = total_return
                        best_decision = investment
                
                V[t][b] = best_value
                D[t][b] = best_decision
        
        total_returns.append(V[0][budget_steps - 1])
    
    # Plot the results
    plt.figure(figsize=(10, 6))
    plt.plot(budget_range, total_returns, 'o-', linewidth=2)
    plt.title('Sensitivity Analysis: Effect of Total Budget on Maximum Return')
    plt.xlabel('Total Budget ($ Million)')
    plt.ylabel('Maximum Expected Return ($ Million)')
    plt.grid(True)
    plt.tight_layout()
    plt.show()
    
    # Calculate and display ROI for each budget level
    roi = [ret/budget for ret, budget in zip(total_returns, budget_range)]
    
    plt.figure(figsize=(10, 6))
    plt.plot(budget_range, roi, 'o-', linewidth=2, color='green')
    plt.title('Return on Investment (ROI) by Total Budget')
    plt.xlabel('Total Budget ($ Million)')
    plt.ylabel('ROI (Return/Budget)')
    plt.grid(True)
    plt.tight_layout()
    plt.show()

# Perform sensitivity analysis
budget_sensitivity_analysis()

# Problem parameters
total_budget = 10  # Total budget in million dollars
periods = 5  # Number of investment periods (years)
discretization = 100  # How finely we discretize the budget
step_size = total_budget / discretization

# Generate return functions for each period
def create_return_functions():
    # Parameters for return functions (a_t, b_t)
    params = [
        (4.0, 0.2),  # Period 1: Higher initial returns but quicker diminishment
        (3.5, 0.1),  # Period 2: Moderate returns with slower diminishment
        (5.0, 0.3),  # Period 3: Highest initial returns but quick diminishment
        (3.0, 0.05), # Period 4: Lower returns but very slow diminishment
        (4.5, 0.15)  # Period 5: Good returns with moderate diminishment
    ]
    
    # Create return functions for each period
    return_functions = []
    for a, b in params:
        def f(x, a=a, b=b):
            if x == 0:
                return 0
            return a * np.sqrt(x) - b * x
        return_functions.append(f)
    
    return return_functions

def solve_dp():
    # Initialize value table and decision table
    budget_steps = int(total_budget / step_size) + 1
    V = np.zeros((periods + 1, budget_steps))
    D = np.zeros((periods, budget_steps))
    
    # DP backward recursion
    for t in range(periods - 1, -1, -1):
        for b in range(budget_steps):
            budget = b * step_size
            
            # Try different investment amounts for period t
            best_value = -float('inf')
            best_decision = 0
            
            for i in range(b + 1):
                investment = i * step_size
                remaining = budget - investment
                remaining_idx = int(remaining / step_size)
                
                current_return = return_functions[t](investment)
                future_return = V[t + 1][remaining_idx]
                total_return = current_return + future_return
                
                if total_return > best_value:
                    best_value = total_return
                    best_decision = investment
            
            V[t][b] = best_value
            D[t][b] = best_decision
    
    # Extract the optimal policy
    optimal_policy = []
    remaining_budget = total_budget
    
    for t in range(periods):
        b_idx = int(remaining_budget / step_size)
        investment = D[t][b_idx]
        optimal_policy.append(investment)
        remaining_budget -= investment
    
    return V, D, optimal_policy

import pulp
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Problem data
suppliers = ['S1', 'S2']
dcs = ['D1', 'D2']
customers = ['C1', 'C2', 'C3']
demand = {'C1': 100, 'C2': 150, 'C3': 200}
capacity = {'D1': 250, 'D2': 300}
w = 0.5  # Weight for emissions

# Transportation costs
cost_sd = {
    ('S1', 'D1'): 10, ('S1', 'D2'): 15,
    ('S2', 'D1'): 12, ('S2', 'D2'): 8
}
cost_dc = {
    ('D1', 'C1'): 5, ('D1', 'C2'): 7, ('D1', 'C3'): 9,
    ('D2', 'C1'): 6, ('D2', 'C2'): 4, ('D2', 'C3'): 8
}

# Emissions (kg CO2/unit)
emission_sd = {
    ('S1', 'D1'): 0.2, ('S1', 'D2'): 0.3,
    ('S2', 'D1'): 0.25, ('S2', 'D2'): 0.15
}
emission_dc = {
    ('D1', 'C1'): 0.1, ('D1', 'C2'): 0.12, ('D1', 'C3'): 0.14,
    ('D2', 'C1'): 0.11, ('D2', 'C2'): 0.09, ('D2', 'C3'): 0.13
}

# Initialize the model
model = pulp.LpProblem("Supply_Chain_Sustainability", pulp.LpMinimize)

# Decision variables
x = pulp.LpVariable.dicts("x", [(s, d) for s in suppliers for d in dcs], lowBound=0, cat='Continuous')
y = pulp.LpVariable.dicts("y", [(d, c) for d in dcs for c in customers], lowBound=0, cat='Continuous')

# Objective function
model += (
    pulp.lpSum(cost_sd[s, d] * x[s, d] for s in suppliers for d in dcs) +
    pulp.lpSum(cost_dc[d, c] * y[d, c] for d in dcs for c in customers) +
    w * (
        pulp.lpSum(emission_sd[s, d] * x[s, d] for s in suppliers for d in dcs) +
        pulp.lpSum(emission_dc[d, c] * y[d, c] for d in dcs for c in customers)
    )
)

# Constraints
# Demand satisfaction
for c in customers:
    model += pulp.lpSum(y[d, c] for d in dcs) == demand[c]

# Flow balance at DCs
for d in dcs:
    model += pulp.lpSum(x[s, d] for s in suppliers) == pulp.lpSum(y[d, c] for c in customers)

# DC capacity
for d in dcs:
    model += pulp.lpSum(x[s, d] for s in suppliers) <= capacity[d]

# Solve the model
model.solve()

# Extract results
status = pulp.LpStatus[model.status]
total_cost = pulp.value(model.objective)
x_values = {(s, d): x[s, d].varValue for s in suppliers for d in dcs}
y_values = {(d, c): y[d, c].varValue for d in dcs for c in customers}

# Print results
print(f"Status: {status}")
print(f"Total Cost (including emissions penalty): {total_cost:.2f}")
print("\nShipments from Suppliers to DCs:")
for s in suppliers:
    for d in dcs:
        if x_values[s, d] > 0:
            print(f"{s} -> {d}: {x_values[s, d]:.2f} units")
print("\nShipments from DCs to Customers:")
for d in dcs:
    for c in customers:
        if y_values[d, c] > 0:
            print(f"{d} -> {c}: {y_values[d, c]:.2f} units")

# Visualization
# Prepare data for heatmaps
x_matrix = np.array([[x_values[s, d] for d in dcs] for s in suppliers])
y_matrix = np.array([[y_values[d, c] for c in customers] for d in dcs])

# Plot heatmap for supplier to DC shipments
plt.figure(figsize=(10, 4))
sns.heatmap(x_matrix, annot=True, fmt=".1f", xticklabels=dcs, yticklabels=suppliers, cmap="YlGnBu")
plt.title("Shipments from Suppliers to DCs (Units)")
plt.xlabel("Distribution Centers")
plt.ylabel("Suppliers")
plt.savefig("supplier_to_dc.png")
plt.show()

# Plot heatmap for DC to customer shipments
plt.figure(figsize=(10, 4))
sns.heatmap(y_matrix, annot=True, fmt=".1f", xticklabels=customers, yticklabels=dcs, cmap="YlOrRd")
plt.title("Shipments from DCs to Customers (Units)")
plt.xlabel("Customers")
plt.ylabel("Distribution Centers")
plt.savefig("dc_to_customer.png")
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from scipy.optimize import minimize_scalar

# Step 1: Simulated historical data (price vs demand)
np.random.seed(0)
prices = np.linspace(5, 20, 30)
true_a = 100
true_b = 3
demand = true_a - true_b * prices + np.random.normal(0, 5, size=prices.shape)

# Step 2: Fit linear regression: D(p) = a - b * p
X = prices.reshape(-1, 1)
y = demand
model = LinearRegression()
model.fit(X, y)
a = model.intercept_
b = -model.coef_[0]  # reverse sign to match D(p) = a - b*p

# Step 3: Revenue function R(p) = p * D(p) = p * (a - b * p)
def revenue(p):
    return -(p * (a - b * p))  # Negative for minimization

# Step 4: Find optimal price
result = minimize_scalar(revenue, bounds=(5, 20), method='bounded')
optimal_price = result.x
optimal_revenue = -result.fun
optimal_demand = a - b * optimal_price

# Step 5: Plot results
p_range = np.linspace(5, 20, 100)
predicted_demand = a - b * p_range
revenue_curve = p_range * predicted_demand

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

# Demand curve
plt.subplot(1, 2, 1)
plt.scatter(prices, demand, label='Observed Demand')
plt.plot(p_range, predicted_demand, color='green', label='Predicted Demand')
plt.axvline(optimal_price, color='red', linestyle='--', label=f'Optimal Price = ${optimal_price:.2f}')
plt.xlabel('Price')
plt.ylabel('Demand')
plt.title('Demand vs. Price')
plt.legend()

# Revenue curve
plt.subplot(1, 2, 2)
plt.plot(p_range, revenue_curve, color='blue', label='Revenue Curve')
plt.axvline(optimal_price, color='red', linestyle='--', label=f'Optimal Revenue = ${optimal_revenue:.2f}')
plt.xlabel('Price')
plt.ylabel('Revenue')
plt.title('Revenue vs. Price')
plt.legend()

plt.tight_layout()
plt.show()

prices = np.linspace(5, 20, 30)
true_a = 100
true_b = 3
demand = true_a - true_b * prices + np.random.normal(0, 5, size=prices.shape)

X = prices.reshape(-1, 1)
model = LinearRegression()
model.fit(X, y)

def revenue(p):
    return -(p * (a - b * p))  # for minimization

result = minimize_scalar(revenue, bounds=(5, 20), method='bounded')

import numpy as np
import matplotlib.pyplot as plt
from math import factorial

# Parameters
λ = 20  # calls per hour
μ = 5   # service rate per agent (calls per hour)
target_wait = 1 / 60  # target waiting time in hours (1 minute)

def P0(λ, μ, c):
    rho = λ / (c * μ)
    sum1 = sum((λ/μ)**n / factorial(n) for n in range(c))
    sum2 = ((λ/μ)**c / factorial(c)) * (1 / (1 - rho)) if rho < 1 else float('inf')
    return 1 / (sum1 + sum2)

def Lq(λ, μ, c):
    rho = λ / (c * μ)
    if rho >= 1:
        return float('inf')
    p0 = P0(λ, μ, c)
    return (p0 * (λ/μ)**c * rho) / (factorial(c) * (1 - rho)**2)

def Wq(λ, μ, c):
    return Lq(λ, μ, c) / λ

# Try various numbers of staff
staff_range = range(1, 21)
waiting_times = [Wq(λ, μ, c) for c in staff_range]

# Find minimum staff for target waiting time
for c, wq in zip(staff_range, waiting_times):
    if wq < target_wait:
        optimal_staff = c
        break

print(f"✅ Optimal number of staff: {optimal_staff}")

plt.figure(figsize=(10, 6))
plt.plot(staff_range, [w * 60 for w in waiting_times], marker='o')  # convert to minutes
plt.axhline(1, color='red', linestyle='--', label='Target wait time (1 minute)')
plt.axvline(optimal_staff, color='green', linestyle='--', label=f'Optimal staff = {optimal_staff}')
plt.xlabel("Number of Staff")
plt.ylabel("Average Waiting Time in Queue (minutes)")
plt.title("Waiting Time vs. Number of Staff")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define the items
items = [
    {"name": "Water Filter", "cost": 4, "value": 10},
    {"name": "First Aid Kit", "cost": 3, "value": 7},
    {"name": "Tent", "cost": 7, "value": 12},
    {"name": "Map", "cost": 2, "value": 4},
    {"name": "Compass", "cost": 2, "value": 5},
    {"name": "Flashlight", "cost": 3, "value": 6}
]

# Budget
budget = 15
n = len(items)

# Initialize DP table
dp = np.zeros((n + 1, budget + 1), dtype=int)

# Fill DP table
for i in range(1, n + 1):
    for w in range(budget + 1):
        if items[i-1]["cost"] <= w:
            dp[i][w] = max(dp[i-1][w], dp[i-1][w - items[i-1]["cost"]] + items[i-1]["value"])
        else:
            dp[i][w] = dp[i-1][w]

# Traceback to find selected items
w = budget
selected = []
for i in range(n, 0, -1):
    if dp[i][w] != dp[i-1][w]:
        selected.append(items[i-1])
        w -= items[i-1]["cost"]

# Show result
print("Selected items:")
for item in selected:
    print(f"- {item['name']} (Cost: ${item['cost']}, Utility: {item['value']})")

print(f"Total cost: ${sum(i['cost'] for i in selected)}")
print(f"Total utility: {sum(i['value'] for i in selected)}")

plt.figure(figsize=(10, 6))
plt.plot(range(budget + 1), dp[n], marker='o', color='green')
plt.title("Maximum Utility vs. Budget")
plt.xlabel("Budget ($)")
plt.ylabel("Maximum Utility")
plt.grid(True)
plt.show()

import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib.colors import LinearSegmentedColormap

# Set the style for our plots
plt.style.use('fivethirtyeight')
sns.set_palette("viridis")

def create_sample_network():
    """
    Create a sample directed network representing a citation network
    between academic papers or web pages.
    """
    # Create a directed graph
    G = nx.DiGraph()
    
    # Add nodes (representing web pages or academic papers)
    nodes = ["A", "B", "C", "D", "E", "F", "G", "H"]
    G.add_nodes_from(nodes)
    
    # Add edges (representing links between pages)
    edges = [
        ("A", "B"), ("A", "C"), ("A", "D"),
        ("B", "C"), ("B", "E"),
        ("C", "A"), ("C", "F"),
        ("D", "B"), ("D", "G"),
        ("E", "D"), ("E", "H"),
        ("F", "G"), ("F", "H"),
        ("G", "E"),
        ("H", "G")
    ]
    G.add_edges_from(edges)
    
    return G

def calculate_pagerank_matrix(G, alpha=0.85, max_iter=100, tol=1e-6):
    """
    Calculate PageRank using the power iteration method.
    
    Parameters:
    - G: NetworkX DiGraph
    - alpha: damping factor (default 0.85)
    - max_iter: maximum number of iterations (default 100)
    - tol: convergence tolerance (default 1e-6)
    
    Returns:
    - PageRank values dictionary
    - Number of iterations performed
    - Convergence history
    """
    n = len(G)
    
    # Create adjacency matrix
    A = nx.to_numpy_array(G, nodelist=list(G.nodes())).T
    
    # Create out-degree matrix
    out_degrees = np.array([max(G.out_degree(node), 1) for node in G.nodes()])
    D_inv = np.diag(1.0 / out_degrees)
    
    # Create transition probability matrix
    M = A @ D_inv
    
    # Initialize PageRank vector
    pr = np.ones(n) / n
    
    # For tracking convergence
    convergence_history = [np.copy(pr)]
    
    # Power iteration method
    for iteration in range(max_iter):
        prev_pr = np.copy(pr)
        
        # PageRank formula: (1-alpha)/n + alpha * M * pr
        pr = (1 - alpha) / n + alpha * M @ pr
        
        # Normalize to ensure sum equals 1
        pr = pr / np.sum(pr)
        
        # Store for convergence history
        convergence_history.append(np.copy(pr))
        
        # Check for convergence
        if np.linalg.norm(pr - prev_pr, 1) < tol:
            break
    
    # Create a dictionary with node names and their PageRank values
    pagerank = {node: rank for node, rank in zip(G.nodes(), pr)}
    
    return pagerank, iteration + 1, convergence_history

def visualize_network_with_pagerank(G, pagerank):
    """
    Visualize the network with node sizes proportional to PageRank values.
    """
    plt.figure(figsize=(12, 10))
    
    # Create a custom colormap from light to dark blue
    colors = ["#9ecae1", "#6baed6", "#4292c6", "#2171b5", "#084594"]
    custom_cmap = LinearSegmentedColormap.from_list("custom_blues", colors)
    
    # Node positions using spring layout
    pos = nx.spring_layout(G, seed=42)
    
    # Scale PageRank values for visualization
    node_sizes = [pagerank[node] * 10000 for node in G.nodes()]
    node_colors = [pagerank[node] for node in G.nodes()]
    
    # Draw nodes
    nodes = nx.draw_networkx_nodes(
        G, pos, 
        node_size=node_sizes,
        node_color=node_colors,
        cmap=custom_cmap,
        alpha=0.9
    )
    
    # Draw edges with arrows
    nx.draw_networkx_edges(
        G, pos, 
        width=1.5, 
        alpha=0.7, 
        edge_color='gray',
        arrowsize=15,
        connectionstyle='arc3,rad=0.1'
    )
    
    # Draw labels
    nx.draw_networkx_labels(
        G, pos, 
        font_size=12, 
        font_weight='bold',
        font_family='sans-serif'
    )
    
    # Add a colorbar
    plt.colorbar(nodes, label='PageRank Score', shrink=0.8)
    
    plt.title("Network Visualization with PageRank Importance", fontsize=18)
    plt.axis('off')
    plt.tight_layout()
    plt.savefig('network_pagerank.png', dpi=300, bbox_inches='tight')
    plt.show()
    
def visualize_pagerank_distribution(pagerank):
    """
    Create a bar chart showing the PageRank distribution.
    """
    plt.figure(figsize=(12, 6))
    
    # Sort PageRank values
    sorted_pagerank = {k: v for k, v in sorted(pagerank.items(), key=lambda item: item[1], reverse=True)}
    
    # Plot bar chart
    bars = plt.bar(sorted_pagerank.keys(), sorted_pagerank.values(), color=sns.color_palette("viridis", len(pagerank)))
    
    # Add value labels on top of bars
    for bar in bars:
        height = bar.get_height()
        plt.text(
            bar.get_x() + bar.get_width()/2., 
            height * 1.01,
            f'{height:.4f}',
            ha='center', 
            va='bottom', 
            fontsize=10,
            rotation=0
        )
    
    plt.title('PageRank Score Distribution', fontsize=18)
    plt.xlabel('Node', fontsize=14)
    plt.ylabel('PageRank Score', fontsize=14)
    plt.grid(axis='y', linestyle='--', alpha=0.7)
    plt.tight_layout()
    plt.savefig('pagerank_distribution.png', dpi=300, bbox_inches='tight')
    plt.show()

def plot_convergence(convergence_history, node_names):
    """
    Plot the convergence of PageRank values over iterations.
    """
    plt.figure(figsize=(12, 6))
    
    # Convert list of arrays to 2D array for easier indexing
    convergence_data = np.array(convergence_history)
    
    # Plot each node's PageRank convergence
    for i, node in enumerate(node_names):
        plt.plot(convergence_data[:, i], marker='o', markersize=4, label=f'Node {node}')
    
    plt.title('PageRank Convergence Over Iterations', fontsize=18)
    plt.xlabel('Iteration', fontsize=14)
    plt.ylabel('PageRank Value', fontsize=14)
    plt.grid(True, linestyle='--', alpha=0.7)
    plt.legend(loc='center right')
    plt.tight_layout()
    plt.savefig('pagerank_convergence.png', dpi=300, bbox_inches='tight')
    plt.show()

def comparative_analysis(G):
    """
    Compare our PageRank implementation with NetworkX's built-in PageRank.
    """
    # Calculate PageRank using our implementation
    our_pagerank, iterations, _ = calculate_pagerank_matrix(G)
    
    # Calculate PageRank using NetworkX's implementation
    nx_pagerank = nx.pagerank(G, alpha=0.85)
    
    # Create a DataFrame for comparison
    comparison_data = {
        'Node': list(our_pagerank.keys()),
        'Our Implementation': list(our_pagerank.values()),
        'NetworkX Implementation': [nx_pagerank[node] for node in our_pagerank.keys()]
    }
    
    # Plot comparison
    plt.figure(figsize=(12, 6))
    bar_width = 0.35
    x = np.arange(len(comparison_data['Node']))
    
    # Plot bars
    plt.bar(x - bar_width/2, comparison_data['Our Implementation'], bar_width, label='Our Implementation')
    plt.bar(x + bar_width/2, comparison_data['NetworkX Implementation'], bar_width, label='NetworkX Implementation')
    
    # Add details
    plt.title('PageRank Implementation Comparison', fontsize=18)
    plt.xlabel('Node', fontsize=14)
    plt.ylabel('PageRank Score', fontsize=14)
    plt.xticks(x, comparison_data['Node'])
    plt.legend()
    plt.grid(axis='y', linestyle='--', alpha=0.7)
    plt.tight_layout()
    plt.savefig('pagerank_comparison.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    # Calculate absolute difference
    diff = np.abs(np.array(list(our_pagerank.values())) - np.array(list(nx_pagerank.values())))
    max_diff = np.max(diff)
    
    print(f"Maximum absolute difference between implementations: {max_diff:.10f}")
    print(f"Number of iterations until convergence: {iterations}")
    
    return our_pagerank, nx_pagerank

# Main execution
def main():
    # Create our sample network
    G = create_sample_network()
    
    # Calculate PageRank
    pagerank, iterations, convergence_history = calculate_pagerank_matrix(G)
    
    # Print results
    print("PageRank Results:")
    for node, rank in sorted(pagerank.items(), key=lambda x: x[1], reverse=True):
        print(f"Node {node}: {rank:.6f}")
    
    print(f"\nNumber of iterations until convergence: {iterations}")
    
    # Visualize the network with PageRank values
    visualize_network_with_pagerank(G, pagerank)
    
    # Visualize PageRank distribution
    visualize_pagerank_distribution(pagerank)
    
    # Plot convergence
    plot_convergence(convergence_history, list(G.nodes()))
    
    # Compare with NetworkX implementation
    our_pagerank, nx_pagerank = comparative_analysis(G)

if __name__ == "__main__":
    main()```


## Code Explanation

Let's break down the key components of the code:

### Network Creation

We start by creating a sample directed network with 8 nodes (A through H), where each node could represent a web page or an academic paper. 
The edges represent links or citations between these nodes.

```python
def create_sample_network():
    G = nx.DiGraph()
    
    # Add nodes (representing web pages or academic papers)
    nodes = ["A", "B", "C", "D", "E", "F", "G", "H"]
    G.add_nodes_from(nodes)
    
    # Add edges (representing links between pages)
    edges = [
        ("A", "B"), ("A", "C"), ("A", "D"),
        ("B", "C"), ("B", "E"),
        # ... other connections
    ]
    G.add_edges_from(edges)
    
    return G

def calculate_pagerank_matrix(G, alpha=0.85, max_iter=100, tol=1e-6):

# PageRank formula: (1-alpha)/n + alpha * M * pr
pr = (1 - alpha) / n + alpha * M @ pr

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

# Parameters
mu = 100            # average daily demand
sigma = 20          # standard deviation of daily demand
L = 5               # lead time in days
h = 2               # holding cost per unit per day
p = 50              # stockout penalty per unit
service_level = 0.95

# Lead time demand distribution
mu_L = mu * L
sigma_L = sigma * np.sqrt(L)

# Find z-score for target service level
z = stats.norm.ppf(service_level)

# Reorder point R
R = mu_L + z * sigma_L

# Safety stock
safety_stock = z * sigma_L

# Expected shortage (standard normal loss function)
L_z = stats.norm.pdf(z) - z * (1 - service_level)

# Expected units short per cycle
E_short = sigma_L * L_z

# Total expected cost function over order quantity Q
def total_cost(Q):
    cycle_stock = Q / 2
    holding_cost = h * (cycle_stock + safety_stock)
    stockout_cost = p * E_short * (mu / Q)
    return holding_cost + stockout_cost

# Search for optimal order quantity Q
Q_values = np.arange(50, 1000, 10)
costs = [total_cost(Q) for Q in Q_values]
optimal_index = np.argmin(costs)
optimal_Q = Q_values[optimal_index]
min_cost = costs[optimal_index]

print(f"Reorder Point (R): {R:.2f}")
print(f"Safety Stock: {safety_stock:.2f}")
print(f"Optimal Order Quantity (Q): {optimal_Q}")
print(f"Minimum Total Cost: ${min_cost:.2f}")

# Plotting
plt.figure(figsize=(10, 6))
plt.plot(Q_values, costs, label="Total Cost", color='blue')
plt.axvline(optimal_Q, color='red', linestyle='--', label=f"Optimal Q = {optimal_Q}")
plt.title("Total Expected Cost vs. Order Quantity")
plt.xlabel("Order Quantity (Q)")
plt.ylabel("Total Expected Cost ($)")
plt.legend()
plt.grid(True)
plt.show()

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

# Step 1: Define the cost matrix
cost_matrix = np.array([
    [90, 76, 75, 80],   # A
    [35, 85, 55, 65],   # B
    [125, 95, 90, 105], # C
    [45, 110, 95, 115]  # D
])

# Step 2: Solve the assignment problem
row_ind, col_ind = linear_sum_assignment(cost_matrix)
total_cost = cost_matrix[row_ind, col_ind].sum()

# St 3: Output the result
staff = ['A', 'B', 'C', 'D']
jobs = ['Job 1', 'Job 2', 'Job 3', 'Job 4']

print("Optimal Assignment:")
for i in range(len(row_ind)):
    print(f"  Staff {staff[row_ind[i]]} → {jobs[col_ind[i]]} (Cost: {cost_matrix[row_ind[i], col_ind[i]]})")

print(f"\nTotal Minimum Cost: {total_cost}")

# Step 4: Visualize with heatmap
plt.figure(figsize=(8, 6))
ax = sns.heatmap(cost_matrix, annot=True, fmt="d", cmap="YlGnBu", cbar=True,
                 xticklabels=jobs, yticklabels=staff)

# Highlight the assigned cells
for i in range(len(row_ind)):
    ax.add_patch(plt.Rectangle((col_ind[i], row_ind[i]), 1, 1, fill=False, edgecolor='red', lw=3))

plt.title(f"Optimal Assignment (Total Cost = {total_cost})")
plt.xlabel("Jobs")
plt.ylabel("Staff")
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
from scipy.stats import norm
from scipy.spatial import distance_matrix
import random
from itertools import permutations
import networkx as nx
from matplotlib.patches import Patch

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

# Problem parameters
num_customers = 10
num_vehicles = 3
vehicle_capacity = 30
depot_coords = (0, 0)
confidence_level = 0.9  # 90% confidence that routes won't fail due to capacity

# Generate customer locations
customer_x = np.random.uniform(-10, 10, num_customers)
customer_y = np.random.uniform(-10, 10, num_customers)
customer_coords = list(zip(customer_x, customer_y))

# Add depot as node 0
all_coords = [depot_coords] + customer_coords

# Generate stochastic demands (normal distribution)
demand_means = np.random.uniform(5, 15, num_customers)
demand_stds = demand_means * 0.3  # Standard deviation as 30% of mean

# Calculate distance matrix
dist_matrix = np.zeros((num_customers + 1, num_customers + 1))
for i in range(num_customers + 1):
    for j in range(num_customers + 1):
        if i != j:
            dist_matrix[i, j] = np.sqrt((all_coords[i][0] - all_coords[j][0])**2 +
                                        (all_coords[i][1] - all_coords[j][1])**2)

# Function to calculate route length
def calculate_route_length(route, dist_mat):
    length = 0
    prev_node = 0  # Start at depot
    for node in route:
        length += dist_mat[prev_node, node]
        prev_node = node
    length += dist_mat[prev_node, 0]  # Return to depot
    return length

# Function to check if a route satisfies capacity constraints with given confidence
def check_capacity_constraint(route, means, stds, capacity, confidence):
    route_mean_demand = sum(means[i-1] for i in route)
    route_var_demand = sum(stds[i-1]**2 for i in route)
    route_std_demand = np.sqrt(route_var_demand)

    # Calculate the critical value based on confidence level
    z_value = norm.ppf(confidence)

    # The route will satisfy capacity with the given confidence if:
    # mean + z * std <= capacity
    effective_demand = route_mean_demand + z_value * route_std_demand

    return effective_demand <= capacity, effective_demand

# Savings algorithm for SVRP
def savings_algorithm_svrp(dist_mat, means, stds, capacity, confidence, num_nodes):
    # Calculate savings for each pair of customers
    savings = {}
    for i in range(1, num_nodes):
        for j in range(1, num_nodes):
            if i != j:
                # Savings formula: s(i,j) = d(0,i) + d(0,j) - d(i,j)
                savings[(i, j)] = dist_mat[0, i] + dist_mat[0, j] - dist_mat[i, j]

    # Sort savings in descending order
    sorted_savings = sorted(savings.items(), key=lambda x: x[1], reverse=True)

    # Initialize routes: one route per customer
    routes = [[i] for i in range(1, num_nodes)]

    # Try to merge routes based on savings
    for (i, j), saving in sorted_savings:
        # Find routes containing i and j
        route_i = None
        route_j = None
        for r in routes:
            if i in r and i == r[-1]:
                route_i = r
            if j in r and j == r[0]:
                route_j = r

        # If i and j are in different routes and at the right positions
        if route_i and route_j and route_i != route_j:
            # Check if merging would violate capacity constraints
            merged_route = route_i + route_j
            feasible, _ = check_capacity_constraint(merged_route, means, stds, capacity, confidence)

            if feasible:
                # Merge routes
                routes.remove(route_i)
                routes.remove(route_j)
                routes.append(merged_route)

    return routes

# Apply the savings algorithm
routes = savings_algorithm_svrp(dist_matrix, demand_means, demand_stds,
                               vehicle_capacity, confidence_level, num_customers + 1)

# Evaluate the solution
total_distance = 0
route_demands = []
route_risks = []

for route in routes:
    route_length = calculate_route_length(route, dist_matrix)
    total_distance += route_length

    # Calculate route statistics
    feasible, effective_demand = check_capacity_constraint(route, demand_means, demand_stds,
                                                          vehicle_capacity, confidence_level)
    route_mean_demand = sum(demand_means[i-1] for i in route)
    route_std_demand = np.sqrt(sum(demand_stds[i-1]**2 for i in route))

    # Calculate risk of exceeding capacity
    risk = 1 - norm.cdf((vehicle_capacity - route_mean_demand) / route_std_demand)

    route_demands.append((route_mean_demand, route_std_demand, effective_demand))
    route_risks.append(risk)

# Output results
print(f"Number of routes: {len(routes)}")
print(f"Total expected distance: {total_distance:.2f}")

# Create dataframe for route statistics
route_stats = []
for i, route in enumerate(routes):
    mean_demand, std_demand, effective_demand = route_demands[i]
    route_stats.append({
        'Route': i+1,
        'Customers': route,
        'Mean Demand': mean_demand,
        'Std Demand': std_demand,
        'Effective Demand': effective_demand,
        'Risk of Exceeding Capacity': route_risks[i] * 100,  # Convert to percentage
        'Distance': calculate_route_length(route, dist_matrix)
    })

route_stats_df = pd.DataFrame(route_stats)
print("\nRoute Statistics:")
print(route_stats_df[['Route', 'Customers', 'Mean Demand', 'Effective Demand',
                      'Risk of Exceeding Capacity', 'Distance']])

# Visualize the solution
plt.figure(figsize=(12, 10))

# Plot customers and depot
plt.scatter(depot_coords[0], depot_coords[1], c='red', s=200, marker='*', label='Depot')
plt.scatter(customer_x, customer_y, c='blue', s=100, label='Customers')

# Add customer labels
for i, (x, y) in enumerate(customer_coords):
    plt.annotate(f"{i+1}", (x, y), fontsize=12)

# Define colors for routes
colors = ['green', 'purple', 'orange', 'brown', 'pink', 'gray', 'olive', 'cyan']

# Plot routes
for i, route in enumerate(routes):
    color = colors[i % len(colors)]

    # Plot route from depot to first customer
    plt.plot([depot_coords[0], customer_coords[route[0]-1][0]],
             [depot_coords[1], customer_coords[route[0]-1][1]],
             c=color, linewidth=2)

    # Plot route between customers
    for j in range(len(route)-1):
        plt.plot([customer_coords[route[j]-1][0], customer_coords[route[j+1]-1][0]],
                 [customer_coords[route[j]-1][1], customer_coords[route[j+1]-1][1]],
                 c=color, linewidth=2)

    # Plot route from last customer back to depot
    plt.plot([customer_coords[route[-1]-1][0], depot_coords[0]],
             [customer_coords[route[-1]-1][1], depot_coords[1]],
             c=color, linewidth=2,
             label=f"Route {i+1}: {route}")

plt.title('Stochastic Vehicle Routing Problem Solution', fontsize=16)
plt.xlabel('X Coordinate', fontsize=14)
plt.ylabel('Y Coordinate', fontsize=14)
plt.grid(True, linestyle='--', alpha=0.7)
plt.legend(loc='best', fontsize=12)
plt.axis('equal')

# Generate legend for routes
route_patches = []
for i, route in enumerate(routes):
    color = colors[i % len(colors)]
    route_patches.append(Patch(color=color, label=f"Route {i+1}: {route}"))

plt.legend(handles=[
    plt.Line2D([0], [0], marker='*', color='w', markerfacecolor='red', markersize=15, label='Depot'),
    plt.Line2D([0], [0], marker='o', color='w', markerfacecolor='blue', markersize=10, label='Customers')
] + route_patches, loc='best', fontsize=12)

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

# Visualize demand uncertainty
plt.figure(figsize=(14, 8))

# Plot demand distributions for each customer
x = np.linspace(0, 25, 100)
for i in range(num_customers):
    y = norm.pdf(x, demand_means[i], demand_stds[i])
    plt.plot(x, y, label=f"Customer {i+1}")

plt.axvline(x=vehicle_capacity, color='r', linestyle='--', label='Vehicle Capacity')
plt.title('Customer Demand Distributions', fontsize=16)
plt.xlabel('Demand', fontsize=14)
plt.ylabel('Probability Density', fontsize=14)
plt.grid(True, linestyle='--', alpha=0.7)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig('demand_distributions.png', dpi=300)
plt.show()

# Visualize route demand uncertainty
plt.figure(figsize=(12, 8))

# Plot cumulative route demand distributions
x = np.linspace(0, 60, 100)
for i, route in enumerate(routes):
    mean_demand, std_demand, _ = route_demands[i]
    y = norm.pdf(x, mean_demand, std_demand)
    plt.plot(x, y, label=f"Route {i+1}: {route}")

plt.axvline(x=vehicle_capacity, color='r', linestyle='--', label='Vehicle Capacity')
plt.title('Route Demand Distributions', fontsize=16)
plt.xlabel('Total Route Demand', fontsize=14)
plt.ylabel('Probability Density', fontsize=14)
plt.grid(True, linestyle='--', alpha=0.7)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig('route_demand_distributions.png', dpi=300)
plt.show()

# Risk analysis visualization
plt.figure(figsize=(10, 6))
risk_df = pd.DataFrame({
    'Route': [f"Route {i+1}" for i in range(len(routes))],
    'Risk (%)': [r * 100 for r in route_risks]
})
sns.barplot(x='Route', y='Risk (%)', data=risk_df)
plt.axhline(y=(1-confidence_level)*100, color='r', linestyle='--',
            label=f'Acceptable Risk ({(1-confidence_level)*100}%)')
plt.title('Risk of Exceeding Vehicle Capacity by Route', fontsize=16)
plt.xlabel('Route', fontsize=14)
plt.ylabel('Risk (%)', fontsize=14)
plt.grid(True, axis='y', linestyle='--', alpha=0.7)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig('route_risks.png', dpi=300)
plt.show()

# Monte Carlo simulation to validate solution
def simulate_demands(means, stds, num_simulations=1000):
    return np.random.normal(
        np.tile(means, (num_simulations, 1)),
        np.tile(stds, (num_simulations, 1))
    )

# Number of simulations
num_simulations = 10000

# Simulate demands for all customers
simulated_demands = simulate_demands(demand_means, demand_stds, num_simulations)

# Calculate route failures
route_failures = np.zeros(len(routes))
for sim in range(num_simulations):
    for i, route in enumerate(routes):
        # Sum the demands for this route in this simulation
        route_demand = sum(max(0, simulated_demands[sim, j-1]) for j in route)
        if route_demand > vehicle_capacity:
            route_failures[i] += 1

# Calculate empirical failure probabilities
empirical_risks = route_failures / num_simulations * 100

# Compare theoretical and empirical risks
risk_comparison = pd.DataFrame({
    'Route': [f"Route {i+1}" for i in range(len(routes))],
    'Theoretical Risk (%)': [r * 100 for r in route_risks],
    'Empirical Risk (%)': empirical_risks
})

plt.figure(figsize=(12, 6))
risk_comparison_melted = pd.melt(risk_comparison, id_vars=['Route'],
                                var_name='Risk Type', value_name='Risk (%)')
sns.barplot(x='Route', y='Risk (%)', hue='Risk Type', data=risk_comparison_melted)
plt.axhline(y=(1-confidence_level)*100, color='r', linestyle='--',
            label=f'Acceptable Risk ({(1-confidence_level)*100}%)')
plt.title('Theoretical vs. Empirical Risk Comparison', fontsize=16)
plt.xlabel('Route', fontsize=14)
plt.ylabel('Risk (%)', fontsize=14)
plt.grid(True, axis='y', linestyle='--', alpha=0.7)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig('risk_comparison.png', dpi=300)
plt.show()

print("\nRisk Comparison (Theoretical vs. Empirical):")
print(risk_comparison)