import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import linprog
import pandas as pd
from matplotlib.patches import Rectangle
import warnings
warnings.filterwarnings('ignore')

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

# Problem Setup: 10 potential locations for renewable energy installation
n_locations = 10
location_names = [f'Site_{i+1}' for i in range(n_locations)]

# Generate realistic data for each location
# Wind speed data (m/s) - affects wind turbine efficiency
wind_speeds = np.random.uniform(4.5, 12.0, n_locations)

# Solar irradiance data (kWh/m²/day) - affects solar panel efficiency  
solar_irradiance = np.random.uniform(3.5, 6.5, n_locations)

# Wind turbine specifications
turbine_capacity = 2.5  # MW per turbine
turbine_cost = 3.0      # Million $ per turbine
turbine_maintenance = 0.1  # Million $ per year per turbine

# Solar panel specifications  
solar_capacity_per_unit = 1.0  # MW per solar unit
solar_cost = 1.5              # Million $ per solar unit
solar_maintenance = 0.05      # Million $ per year per solar unit

# Grid connection costs (varies by location)
grid_connection_cost = np.random.uniform(0.2, 0.8, n_locations)

# Calculate energy output coefficients
# Wind energy output: P = 0.5 * ρ * A * v³ * Cp (simplified)
# We'll use a simplified efficiency curve
def wind_efficiency(wind_speed):
    """Calculate wind turbine efficiency based on wind speed"""
    if wind_speed < 3:
        return 0
    elif wind_speed < 12:
        return 0.35 * (wind_speed / 12) ** 3
    else:
        return 0.35

# Solar efficiency is more linear with irradiance
def solar_efficiency(irradiance):
    """Calculate solar panel efficiency based on irradiance"""
    return 0.18 * (irradiance / 5.0)

# Calculate annual energy output for each location (GWh/year)
wind_output = np.array([wind_efficiency(ws) * turbine_capacity * 8760 / 1000 
                       for ws in wind_speeds])
solar_output = np.array([solar_efficiency(si) * solar_capacity_per_unit * 8760 / 1000 
                        for si in solar_irradiance])

# Total costs per installation (including 5-year maintenance)
wind_total_cost = turbine_cost + 5 * turbine_maintenance + grid_connection_cost
solar_total_cost = solar_cost + 5 * solar_maintenance + grid_connection_cost

# Budget constraint
total_budget = 50.0  # Million $

print("=== Renewable Energy Optimization Problem ===")
print(f"Budget: ${total_budget} million")
print(f"Number of potential sites: {n_locations}")
print(f"Wind turbine capacity: {turbine_capacity} MW")
print(f"Solar unit capacity: {solar_capacity_per_unit} MW")
print("\n=== Location Data ===")

# Create comprehensive data table
data_table = pd.DataFrame({
    'Location': location_names,
    'Wind_Speed_ms': np.round(wind_speeds, 2),
    'Solar_Irradiance_kWh_m2_day': np.round(solar_irradiance, 2),
    'Wind_Output_GWh_year': np.round(wind_output, 2),
    'Solar_Output_GWh_year': np.round(solar_output, 2),
    'Wind_Cost_Million': np.round(wind_total_cost, 2),
    'Solar_Cost_Million': np.round(solar_total_cost, 2),
    'Grid_Connection_Million': np.round(grid_connection_cost, 2)
})

print(data_table.to_string(index=False))

# Mathematical Formulation
print("\n=== Mathematical Formulation ===")
print("Objective Function:")
print("Maximize: Σ(wind_output[i] * x_wind[i] + solar_output[i] * x_solar[i])")
print("\nSubject to:")
print("Σ(wind_cost[i] * x_wind[i] + solar_cost[i] * x_solar[i]) ≤ Budget")
print("x_wind[i], x_solar[i] ∈ {0, 1} for all i")

# Since scipy.optimize.linprog doesn't handle integer programming directly,
# we'll use a relaxed LP and then apply rounding heuristics
# For exact integer solutions, specialized solvers like PuLP would be better

# Prepare linear programming formulation
# Variables: [x_wind_0, x_wind_1, ..., x_wind_9, x_solar_0, x_solar_1, ..., x_solar_9]
c = np.concatenate([-wind_output, -solar_output])  # Negative for maximization
A_ub = np.concatenate([wind_total_cost, solar_total_cost]).reshape(1, -1)
b_ub = np.array([total_budget])

# Bounds: 0 ≤ x ≤ 1 for relaxed problem
bounds = [(0, 1) for _ in range(2 * n_locations)]

# Solve relaxed linear program
print("\n=== Solving Optimization Problem ===")
result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

if result.success:
    print("Optimization successful!")
    
    # Extract solutions
    x_wind_relaxed = result.x[:n_locations]
    x_solar_relaxed = result.x[n_locations:]
    
    # Apply rounding heuristic for integer solution
    # Sort by efficiency ratio (output/cost) and select greedily
    wind_efficiency_ratio = wind_output / wind_total_cost
    solar_efficiency_ratio = solar_output / solar_total_cost
    
    # Create candidate list with efficiency ratios
    candidates = []
    for i in range(n_locations):
        candidates.append(('wind', i, wind_efficiency_ratio[i], wind_total_cost[i], wind_output[i]))
        candidates.append(('solar', i, solar_efficiency_ratio[i], solar_total_cost[i], solar_output[i]))
    
    # Sort by efficiency ratio (descending)
    candidates.sort(key=lambda x: x[2], reverse=True)
    
    # Greedy selection
    selected_wind = np.zeros(n_locations, dtype=int)
    selected_solar = np.zeros(n_locations, dtype=int)
    remaining_budget = total_budget
    total_output = 0
    
    print("\n=== Greedy Selection Process ===")
    for energy_type, location, ratio, cost, output in candidates:
        if cost <= remaining_budget:
            if energy_type == 'wind':
                selected_wind[location] = 1
            else:
                selected_solar[location] = 1
            remaining_budget -= cost
            total_output += output
            print(f"Selected: {energy_type.upper()} at {location_names[location]} "
                  f"(Efficiency: {ratio:.3f}, Cost: ${cost:.2f}M, Output: {output:.2f} GWh/year)")
    
    print(f"\nRemaining budget: ${remaining_budget:.2f} million")
    print(f"Total annual energy output: {total_output:.2f} GWh/year")
    
    # Calculate total investment
    total_investment = total_budget - remaining_budget
    print(f"Total investment: ${total_investment:.2f} million")
    
else:
    print("Optimization failed!")
    selected_wind = np.zeros(n_locations, dtype=int)
    selected_solar = np.zeros(n_locations, dtype=int)

# Create comprehensive results analysis
print("\n=== Final Solution Analysis ===")
solution_df = pd.DataFrame({
    'Location': location_names,
    'Wind_Selected': selected_wind,
    'Solar_Selected': selected_solar,
    'Wind_Speed': np.round(wind_speeds, 2),
    'Solar_Irradiance': np.round(solar_irradiance, 2),
    'Total_Output_GWh': np.round(
        selected_wind * wind_output + selected_solar * solar_output, 2),
    'Total_Cost_Million': np.round(
        selected_wind * wind_total_cost + selected_solar * solar_total_cost, 2)
})

print(solution_df.to_string(index=False))

# Visualization Section
plt.style.use('seaborn-v0_8')
fig = plt.figure(figsize=(20, 15))

# 1. Location Overview Map-style Plot
ax1 = plt.subplot(2, 3, 1)
scatter_wind = ax1.scatter([i for i in range(n_locations) if selected_wind[i]], 
                         [wind_speeds[i] for i in range(n_locations) if selected_wind[i]],
                         s=[wind_output[i]*50 for i in range(n_locations) if selected_wind[i]], 
                         c='blue', alpha=0.7, label='Wind Turbines', marker='^')

scatter_solar = ax1.scatter([i for i in range(n_locations) if selected_solar[i]], 
                          [solar_irradiance[i] for i in range(n_locations) if selected_solar[i]],
                          s=[solar_output[i]*50 for i in range(n_locations) if selected_solar[i]], 
                          c='orange', alpha=0.7, label='Solar Panels', marker='s')

# Add unselected locations in gray
unselected_idx = [i for i in range(n_locations) 
                 if not selected_wind[i] and not selected_solar[i]]
if unselected_idx:
    ax1.scatter(unselected_idx, [wind_speeds[i] for i in unselected_idx],
               c='gray', alpha=0.3, s=30, label='Unselected')

ax1.set_xlabel('Location Index')
ax1.set_ylabel('Resource Quality')
ax1.set_title('Selected Renewable Energy Sites\n(Bubble size = Annual Output)')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 2. Cost vs Output Analysis
ax2 = plt.subplot(2, 3, 2)
wind_costs = [wind_total_cost[i] if selected_wind[i] else 0 for i in range(n_locations)]
solar_costs = [solar_total_cost[i] if selected_solar[i] else 0 for i in range(n_locations)]
wind_outputs = [wind_output[i] if selected_wind[i] else 0 for i in range(n_locations)]
solar_outputs = [solar_output[i] if selected_solar[i] else 0 for i in range(n_locations)]

x_pos = np.arange(n_locations)
width = 0.35

bars1 = ax2.bar(x_pos - width/2, wind_costs, width, label='Wind Cost', color='lightblue', alpha=0.7)
bars2 = ax2.bar(x_pos + width/2, solar_costs, width, label='Solar Cost', color='lightyellow', alpha=0.7)

ax2_twin = ax2.twinx()
line1 = ax2_twin.plot(x_pos, wind_outputs, 'bo-', label='Wind Output', linewidth=2, markersize=6)
line2 = ax2_twin.plot(x_pos, solar_outputs, 'ro-', label='Solar Output', linewidth=2, markersize=6)

ax2.set_xlabel('Location')
ax2.set_ylabel('Cost (Million $)', color='black')
ax2_twin.set_ylabel('Output (GWh/year)', color='black')
ax2.set_title('Cost vs Output by Location')
ax2.set_xticks(x_pos)
ax2.set_xticklabels([f'S{i+1}' for i in range(n_locations)], rotation=45)

# Combine legends
lines1, labels1 = ax2.get_legend_handles_labels()
lines2, labels2 = ax2_twin.get_legend_handles_labels()
ax2.legend(lines1 + lines2, labels1 + labels2, loc='upper left')
ax2.grid(True, alpha=0.3)

# 3. Resource Quality Distribution
ax3 = plt.subplot(2, 3, 3)
selected_locations = [i for i in range(n_locations) 
                     if selected_wind[i] or selected_solar[i]]
unselected_locations = [i for i in range(n_locations) 
                       if not (selected_wind[i] or selected_solar[i])]

if selected_locations:
    ax3.scatter([wind_speeds[i] for i in selected_locations], 
               [solar_irradiance[i] for i in selected_locations],
               c='green', s=100, alpha=0.7, label='Selected Sites')

if unselected_locations:
    ax3.scatter([wind_speeds[i] for i in unselected_locations], 
               [solar_irradiance[i] for i in unselected_locations],
               c='red', s=100, alpha=0.7, label='Unselected Sites')

ax3.set_xlabel('Wind Speed (m/s)')
ax3.set_ylabel('Solar Irradiance (kWh/m²/day)')
ax3.set_title('Resource Quality: Selected vs Unselected')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 4. Budget Utilization
ax4 = plt.subplot(2, 3, 4)
budget_categories = ['Wind\nInvestment', 'Solar\nInvestment', 'Remaining\nBudget']
budget_values = [
    np.sum(selected_wind * wind_total_cost),
    np.sum(selected_solar * solar_total_cost),
    remaining_budget
]
colors = ['skyblue', 'gold', 'lightgray']

wedges, texts, autotexts = ax4.pie(budget_values, labels=budget_categories, 
                                  autopct='%1.1f%%', startangle=90, colors=colors)
ax4.set_title(f'Budget Utilization\n(Total: ${total_budget}M)')

# 5. Efficiency Comparison
ax5 = plt.subplot(2, 3, 5)
efficiency_data = {
    'Technology': ['Wind Turbines', 'Solar Panels'],
    'Avg_Efficiency_Ratio': [
        np.mean([wind_efficiency_ratio[i] for i in range(n_locations) if selected_wind[i]]) 
        if np.sum(selected_wind) > 0 else 0,
        np.mean([solar_efficiency_ratio[i] for i in range(n_locations) if selected_solar[i]]) 
        if np.sum(selected_solar) > 0 else 0
    ],
    'Count': [np.sum(selected_wind), np.sum(selected_solar)]
}

bars = ax5.bar(efficiency_data['Technology'], efficiency_data['Avg_Efficiency_Ratio'], 
               color=['lightblue', 'lightyellow'], alpha=0.8)

# Add count labels on bars
for bar, count in zip(bars, efficiency_data['Count']):
    height = bar.get_height()
    ax5.text(bar.get_x() + bar.get_width()/2., height + 0.01,
             f'Count: {count}', ha='center', va='bottom')

ax5.set_ylabel('Average Efficiency Ratio (GWh/$M)')
ax5.set_title('Technology Efficiency Comparison')
ax5.grid(True, alpha=0.3)

# 6. Timeline Projection
ax6 = plt.subplot(2, 3, 6)
years = np.arange(1, 11)  # 10-year projection
annual_output = total_output
cumulative_output = annual_output * years
cumulative_investment = total_investment + np.arange(1, 11) * 0.1 * total_investment  # 10% annual maintenance

ax6.plot(years, cumulative_output, 'g-', linewidth=3, label='Cumulative Energy Output (GWh)')
ax6_twin = ax6.twinx()
ax6_twin.plot(years, cumulative_investment, 'r--', linewidth=2, label='Cumulative Investment ($M)')

ax6.set_xlabel('Years')
ax6.set_ylabel('Cumulative Energy (GWh)', color='green')
ax6_twin.set_ylabel('Cumulative Cost ($M)', color='red')
ax6.set_title('10-Year Energy & Investment Projection')
ax6.grid(True, alpha=0.3)

# Combine legends
lines1, labels1 = ax6.get_legend_handles_labels()
lines2, labels2 = ax6_twin.get_legend_handles_labels()
ax6.legend(lines1 + lines2, labels1 + labels2, loc='upper left')

plt.tight_layout()
plt.show()

# Mathematical Summary
print("\n" + "="*60)
print("MATHEMATICAL FORMULATION SUMMARY")
print("="*60)

print("\nObjective Function (Maximization):")
print("$$\\max \\sum_{i=1}^{n} \\left( E_{wind}^i \\cdot x_{wind}^i + E_{solar}^i \\cdot x_{solar}^i \\right)$$")

print("\nWhere:")
print("- $E_{wind}^i$: Annual wind energy output at location $i$")
print("- $E_{solar}^i$: Annual solar energy output at location $i$")  
print("- $x_{wind}^i$: Binary decision variable for wind installation at location $i$")
print("- $x_{solar}^i$: Binary decision variable for solar installation at location $i$")

print("\nConstraints:")
print("Budget Constraint:")
print("$$\\sum_{i=1}^{n} \\left( C_{wind}^i \\cdot x_{wind}^i + C_{solar}^i \\cdot x_{solar}^i \\right) \\leq B$$")

print("\nBinary Constraints:")
print("$$x_{wind}^i, x_{solar}^i \\in \\{0,1\\} \\quad \\forall i \\in \\{1,2,...,n\\}$$")

print(f"\nProblem Instance:")
print(f"- $n = {n_locations}$ (number of locations)")
print(f"- $B = {total_budget}$ million USD (budget)")
print(f"- Total variables: {2*n_locations} (binary)")

print(f"\nSolution Quality:")
print(f"- Total annual output: {total_output:.2f} GWh/year")
print(f"- Budget utilization: {(total_investment/total_budget)*100:.1f}%")
print(f"- Wind installations: {np.sum(selected_wind)}")
print(f"- Solar installations: {np.sum(selected_solar)}")
print(f"- Average ROI: {total_output/total_investment:.2f} GWh per million USD")

print("\n" + "="*60)
print("OPTIMIZATION COMPLETE")
print("="*60)

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

# Given constants
P = 1000          # Load in N
L = 2.0           # Length in meters
rho = 7850        # Density of steel (kg/m^3)
g = 9.81          # Gravity (m/s^2)
sigma_yield = 250e6  # Yield stress in Pa
b = 0.05          # Fixed width in meters

# Objective function: weight of the beam
def weight(x):
    h = x[0]
    volume = b * h * L
    return rho * g * volume  # Weight = mass * gravity

# Constraint: maximum bending stress must be less than yield
def stress_constraint(x):
    h = x[0]
    sigma_max = 6 * P * L / (b * h**2)
    return sigma_yield - sigma_max

# Bounds and constraints
bounds = [(0.01, 0.5)]  # height h should be between 1cm and 50cm
constraints = [{'type': 'ineq', 'fun': stress_constraint}]

# Initial guess
x0 = [0.1]

# Run optimization
result = minimize(weight, x0, method='SLSQP', bounds=bounds, constraints=constraints)

# Extract result
h_opt = result.x[0]
w_opt = weight([h_opt])

print(f"Optimal height h = {h_opt:.4f} m")
print(f"Minimum weight = {w_opt:.2f} N")

def weight(x):
    h = x[0]
    volume = b * h * L
    return rho * g * volume

def stress_constraint(x):
    h = x[0]
    sigma_max = 6 * P * L / (b * h**2)
    return sigma_yield - sigma_max

h_vals = np.linspace(0.01, 0.5, 100)
weights = rho * g * b * h_vals * L
stresses = 6 * P * L / (b * h_vals**2)

fig, ax1 = plt.subplots()

ax1.set_xlabel('Height h (m)')
ax1.set_ylabel('Weight (N)', color='tab:blue')
ax1.plot(h_vals, weights, label='Weight', color='tab:blue')
ax1.tick_params(axis='y', labelcolor='tab:blue')

ax2 = ax1.twinx()
ax2.set_ylabel('Max Bending Stress (Pa)', color='tab:red')
ax2.plot(h_vals, stresses, label='Stress', color='tab:red')
ax2.axhline(sigma_yield, color='gray', linestyle='--', label='Yield Stress')
ax2.tick_params(axis='y', labelcolor='tab:red')

plt.title("Beam Weight and Stress vs Height")
fig.tight_layout()
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import random
from collections import defaultdict
import seaborn as sns
from scipy.optimize import differential_evolution
import warnings
warnings.filterwarnings('ignore')

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

class HPProteinFolder:
    """
    HP Model Protein Folding Optimizer
    
    This class implements a simplified protein folding model where:
    - H represents hydrophobic amino acids
    - P represents polar amino acids
    - Goal: Maximize hydrophobic contacts while avoiding collisions
    """
    
    def __init__(self, sequence):
        """
        Initialize the protein folder with an HP sequence
        
        Args:
            sequence (str): String of 'H' and 'P' characters representing the protein
        """
        self.sequence = sequence
        self.length = len(sequence)
        self.moves = [(0, 1), (1, 0), (0, -1), (-1, 0)]  # Up, Right, Down, Left
        self.move_names = ['U', 'R', 'D', 'L']
        
    def decode_folding(self, angles):
        """
        Convert angle representation to 2D coordinates
        
        Args:
            angles (list): List of angles (0-3) representing folding directions
            
        Returns:
            tuple: (coordinates, valid_folding_flag)
        """
        if len(angles) != self.length - 1:
            return None, False
            
        coordinates = [(0, 0)]  # Start at origin
        current_pos = (0, 0)
        occupied = {(0, 0)}
        
        for i, angle in enumerate(angles):
            # Convert continuous angle to discrete move
            move_idx = int(angle) % 4
            dx, dy = self.moves[move_idx]
            new_pos = (current_pos[0] + dx, current_pos[1] + dy)
            
            # Check for collision (self-intersection)
            if new_pos in occupied:
                return None, False
                
            coordinates.append(new_pos)
            occupied.add(new_pos)
            current_pos = new_pos
            
        return coordinates, True
    
    def calculate_energy(self, coordinates):
        """
        Calculate the energy of a folding configuration
        
        Energy = -1 * (number of H-H contacts that are adjacent in space but not in sequence)
        
        Args:
            coordinates (list): List of (x, y) coordinates for each amino acid
            
        Returns:
            float: Energy value (lower is better)
        """
        if coordinates is None:
            return float('inf')  # Invalid folding has infinite energy
            
        energy = 0
        
        # Check all pairs of amino acids
        for i in range(self.length):
            for j in range(i + 2, self.length):  # Skip adjacent in sequence
                # Only count H-H interactions
                if self.sequence[i] == 'H' and self.sequence[j] == 'H':
                    # Check if they are adjacent in 2D space
                    dist = abs(coordinates[i][0] - coordinates[j][0]) + abs(coordinates[i][1] - coordinates[j][1])
                    if dist == 1:  # Manhattan distance = 1 means adjacent
                        energy -= 1  # Favorable interaction
                        
        return energy
    
    def objective_function(self, angles):
        """
        Objective function for optimization
        
        Args:
            angles (array): Array of folding angles
            
        Returns:
            float: Energy to minimize
        """
        coordinates, valid = self.decode_folding(angles)
        if not valid:
            return 1000  # Heavy penalty for invalid folding
            
        return self.calculate_energy(coordinates)
    
    def optimize_folding(self, method='differential_evolution', max_iterations=1000):
        """
        Optimize protein folding using specified optimization method
        
        Args:
            method (str): Optimization method to use
            max_iterations (int): Maximum number of iterations
            
        Returns:
            tuple: (best_angles, best_energy, optimization_history)
        """
        bounds = [(0, 3.99) for _ in range(self.length - 1)]
        history = []
        
        def callback(xk, convergence=None):
            energy = self.objective_function(xk)
            history.append(energy)
            return False
        
        if method == 'differential_evolution':
            result = differential_evolution(
                self.objective_function,
                bounds,
                maxiter=max_iterations // 10,
                popsize=15,
                callback=callback,
                seed=42
            )
            best_angles = result.x
            best_energy = result.fun
            
        return best_angles, best_energy, history
    
    def monte_carlo_optimization(self, max_iterations=10000, temperature=1.0, cooling_rate=0.995):
        """
        Monte Carlo optimization with simulated annealing
        
        Args:
            max_iterations (int): Maximum number of iterations
            temperature (float): Initial temperature
            cooling_rate (float): Temperature cooling rate
            
        Returns:
            tuple: (best_angles, best_energy, history)
        """
        current_angles = [random.uniform(0, 3.99) for _ in range(self.length - 1)]
        current_energy = self.objective_function(current_angles)
        
        best_angles = current_angles.copy()
        best_energy = current_energy
        history = [current_energy]
        
        temp = temperature
        
        for iteration in range(max_iterations):
            # Generate neighbor solution
            new_angles = current_angles.copy()
            idx = random.randint(0, len(new_angles) - 1)
            new_angles[idx] = random.uniform(0, 3.99)
            
            new_energy = self.objective_function(new_angles)
            
            # Accept or reject based on Metropolis criterion
            if new_energy < current_energy or random.random() < np.exp(-(new_energy - current_energy) / temp):
                current_angles = new_angles
                current_energy = new_energy
                
                if current_energy < best_energy:
                    best_angles = current_angles.copy()
                    best_energy = current_energy
            
            history.append(best_energy)
            temp *= cooling_rate
            
            if iteration % 1000 == 0:
                print(f"Iteration {iteration}: Best Energy = {best_energy:.3f}, Temp = {temp:.4f}")
        
        return best_angles, best_energy, history

# Example protein sequences for testing
sequences = {
    'simple': 'HPHPPHHPHH',  # 10 amino acids
    'medium': 'HPHPPHHPHHPPHPPHHPHH',  # 20 amino acids  
    'complex': 'HHPPHPHPHPPHPPHPPHPHPHHPPHPPH'  # 27 amino acids
}

print("=== Protein Folding Optimization Results ===\n")

# Let's work with the medium complexity sequence
sequence = sequences['medium']
print(f"Protein sequence: {sequence}")
print(f"Length: {len(sequence)} amino acids")
print(f"Hydrophobic residues: {sequence.count('H')}")
print(f"Polar residues: {sequence.count('P')}\n")

# Initialize the protein folder
folder = HPProteinFolder(sequence)

# Optimize using differential evolution
print("Running Differential Evolution optimization...")
de_angles, de_energy, de_history = folder.optimize_folding(method='differential_evolution')
de_coords, de_valid = folder.decode_folding(de_angles)

print(f"DE - Best energy: {de_energy:.3f}")
print(f"DE - Valid folding: {de_valid}\n")

# Optimize using Monte Carlo with simulated annealing
print("Running Monte Carlo with Simulated Annealing...")
mc_angles, mc_energy, mc_history = folder.monte_carlo_optimization(max_iterations=5000)
mc_coords, mc_valid = folder.decode_folding(mc_angles)

print(f"MC - Best energy: {mc_energy:.3f}")
print(f"MC - Valid folding: {mc_valid}\n")

# Analysis and visualization
def analyze_folding(coordinates, sequence, title):
    """
    Analyze and display folding statistics
    """
    if coordinates is None:
        print(f"{title}: Invalid folding")
        return
        
    # Calculate statistics
    h_positions = [i for i, aa in enumerate(sequence) if aa == 'H']
    p_positions = [i for i, aa in enumerate(sequence) if aa == 'P']
    
    # Count H-H contacts
    hh_contacts = 0
    for i in range(len(sequence)):
        for j in range(i + 2, len(sequence)):
            if sequence[i] == 'H' and sequence[j] == 'H':
                dist = abs(coordinates[i][0] - coordinates[j][0]) + abs(coordinates[i][1] - coordinates[j][1])
                if dist == 1:
                    hh_contacts += 1
    
    # Calculate compactness (radius of gyration)
    center_x = np.mean([coord[0] for coord in coordinates])
    center_y = np.mean([coord[1] for coord in coordinates])
    rg = np.sqrt(np.mean([(coord[0] - center_x)**2 + (coord[1] - center_y)**2 for coord in coordinates]))
    
    print(f"{title} Analysis:")
    print(f"  H-H contacts: {hh_contacts}")
    print(f"  Energy: {-hh_contacts}")
    print(f"  Radius of gyration: {rg:.3f}")
    print(f"  Span X: {max(coord[0] for coord in coordinates) - min(coord[0] for coord in coordinates)}")
    print(f"  Span Y: {max(coord[1] for coord in coordinates) - min(coord[1] for coord in coordinates)}")
    print()

# Analyze both solutions
analyze_folding(de_coords, sequence, "Differential Evolution")
analyze_folding(mc_coords, sequence, "Monte Carlo")

# Create comprehensive visualization
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Protein Folding Optimization Results', fontsize=16, fontweight='bold')

# Plot 1: DE Folding Structure
if de_coords and de_valid:
    ax = axes[0, 0]
    x_coords = [coord[0] for coord in de_coords]
    y_coords = [coord[1] for coord in de_coords]
    
    # Plot backbone
    ax.plot(x_coords, y_coords, 'k-', linewidth=2, alpha=0.7, label='Backbone')
    
    # Plot amino acids with different colors
    for i, (x, y) in enumerate(de_coords):
        color = 'red' if sequence[i] == 'H' else 'blue'
        marker = 'o' if sequence[i] == 'H' else 's'
        ax.scatter(x, y, c=color, s=150, marker=marker, 
                  edgecolors='black', linewidth=1, zorder=5)
        ax.annotate(f'{i}', (x, y), xytext=(5, 5), textcoords='offset points', 
                   fontsize=8, ha='left')
    
    ax.set_title(f'Differential Evolution\nEnergy: {de_energy:.1f}', fontweight='bold')
    ax.set_xlabel('X coordinate')
    ax.set_ylabel('Y coordinate')
    ax.grid(True, alpha=0.3)
    ax.legend(['Backbone', 'Hydrophobic (H)', 'Polar (P)'])
    ax.set_aspect('equal')

# Plot 2: MC Folding Structure  
if mc_coords and mc_valid:
    ax = axes[0, 1]
    x_coords = [coord[0] for coord in mc_coords]
    y_coords = [coord[1] for coord in mc_coords]
    
    # Plot backbone
    ax.plot(x_coords, y_coords, 'k-', linewidth=2, alpha=0.7)
    
    # Plot amino acids
    for i, (x, y) in enumerate(mc_coords):
        color = 'red' if sequence[i] == 'H' else 'blue'
        marker = 'o' if sequence[i] == 'H' else 's'
        ax.scatter(x, y, c=color, s=150, marker=marker, 
                  edgecolors='black', linewidth=1, zorder=5)
        ax.annotate(f'{i}', (x, y), xytext=(5, 5), textcoords='offset points', 
                   fontsize=8, ha='left')
    
    ax.set_title(f'Monte Carlo\nEnergy: {mc_energy:.1f}', fontweight='bold')
    ax.set_xlabel('X coordinate')
    ax.set_ylabel('Y coordinate')
    ax.grid(True, alpha=0.3)
    ax.set_aspect('equal')

# Plot 3: Convergence Comparison
ax = axes[0, 2]
if len(de_history) > 0:
    ax.plot(range(len(de_history)), de_history, 'g-', linewidth=2, label='Differential Evolution')
if len(mc_history) > 0:
    # Sample MC history for better visualization
    sample_indices = np.linspace(0, len(mc_history)-1, min(200, len(mc_history)), dtype=int)
    sampled_history = [mc_history[i] for i in sample_indices]
    ax.plot(sample_indices, sampled_history, 'r-', linewidth=2, label='Monte Carlo')

ax.set_title('Optimization Convergence', fontweight='bold')
ax.set_xlabel('Iteration')
ax.set_ylabel('Best Energy Found')
ax.legend()
ax.grid(True, alpha=0.3)

# Plot 4: Energy Landscape Analysis
ax = axes[1, 0]
# Generate random foldings to show energy distribution
random_energies = []
for _ in range(1000):
    random_angles = [random.uniform(0, 3.99) for _ in range(len(sequence) - 1)]
    energy = folder.objective_function(random_angles)
    if energy < 100:  # Only valid foldings
        random_energies.append(energy)

if random_energies:
    ax.hist(random_energies, bins=30, alpha=0.7, color='skyblue', edgecolor='black')
    ax.axvline(de_energy, color='green', linestyle='--', linewidth=2, label=f'DE: {de_energy:.1f}')
    ax.axvline(mc_energy, color='red', linestyle='--', linewidth=2, label=f'MC: {mc_energy:.1f}')
    ax.set_title('Energy Distribution\n(Random vs Optimized)', fontweight='bold')
    ax.set_xlabel('Energy')
    ax.set_ylabel('Frequency')
    ax.legend()
    ax.grid(True, alpha=0.3)

# Plot 5: Sequence Analysis
ax = axes[1, 1]
positions = list(range(len(sequence)))
colors = ['red' if aa == 'H' else 'blue' for aa in sequence]
bars = ax.bar(positions, [1]*len(sequence), color=colors, alpha=0.7, edgecolor='black')

ax.set_title('Protein Sequence', fontweight='bold')
ax.set_xlabel('Position')
ax.set_ylabel('Amino Acid Type')
ax.set_xticks(positions[::2])
ax.set_xticklabels(positions[::2])
ax.set_yticks([0.5, 1])
ax.set_yticklabels(['', 'H/P'])

# Add legend
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor='red', alpha=0.7, label='Hydrophobic (H)'),
                  Patch(facecolor='blue', alpha=0.7, label='Polar (P)')]
ax.legend(handles=legend_elements)

# Plot 6: Contact Map
ax = axes[1, 2]
contact_matrix = np.zeros((len(sequence), len(sequence)))

# Fill contact matrix for the best solution (MC in this case)
best_coords = mc_coords if mc_valid else de_coords
if best_coords:
    for i in range(len(sequence)):
        for j in range(len(sequence)):
            if i != j:
                dist = abs(best_coords[i][0] - best_coords[j][0]) + abs(best_coords[i][1] - best_coords[j][1])
                if dist == 1:  # Adjacent in space
                    contact_matrix[i, j] = 1

im = ax.imshow(contact_matrix, cmap='RdYlBu_r', aspect='equal')
ax.set_title('Contact Map\n(Best Solution)', fontweight='bold')
ax.set_xlabel('Residue Index')
ax.set_ylabel('Residue Index')

# Add colorbar
cbar = plt.colorbar(im, ax=ax, shrink=0.8)
cbar.set_label('Contact')

plt.tight_layout()
plt.show()

# Performance comparison
print("=== Performance Comparison ===")
print(f"Differential Evolution:")
print(f"  Final Energy: {de_energy:.3f}")
print(f"  Convergence Steps: {len(de_history)}")

print(f"Monte Carlo + Simulated Annealing:")
print(f"  Final Energy: {mc_energy:.3f}")
print(f"  Total Steps: {len(mc_history)}")

# Determine the winner
if de_valid and mc_valid:
    if de_energy < mc_energy:
        print(f"\n🏆 Winner: Differential Evolution (Energy: {de_energy:.3f})")
    elif mc_energy < de_energy:
        print(f"\n🏆 Winner: Monte Carlo (Energy: {mc_energy:.3f})")
    else:
        print(f"\n🤝 Tie! Both methods achieved energy: {de_energy:.3f}")
else:
    print(f"\n⚠️  Some optimizations produced invalid foldings")

print("\n=== Mathematical Analysis ===")
print("Energy Function:")
print("E = -∑(i,j) δ_HiHj * contact(i,j)")
print("where δ_HiHj = 1 if both residues i,j are hydrophobic and adjacent")
print("Lower energy indicates better folding (more H-H contacts)")

def calculate_energy(self, coordinates):
    energy = 0
    for i in range(self.length):
        for j in range(i + 2, self.length):  # Skip adjacent in sequence
            if self.sequence[i] == 'H' and self.sequence[j] == 'H':
                dist = abs(coordinates[i][0] - coordinates[j][0]) + abs(coordinates[i][1] - coordinates[j][1])
                if dist == 1:  # Manhattan distance = 1 means adjacent
                    energy -= 1  # Favorable interaction
    return energy

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

# Constants
CYCLE_TIME = 60  # Total signal cycle time (seconds)
lambda_NS = 0.6  # Arrival rate per second from North-South
lambda_EW = 0.4  # Arrival rate per second from East-West

# Cost function to minimize: total waiting time
def total_waiting_time(g):
    if g <= 0 or g >= CYCLE_TIME:
        return np.inf  # avoid invalid green times
    red_NS = CYCLE_TIME - g
    red_EW = g
    wait_NS = (lambda_NS * red_NS ** 2) / (2 * g)
    wait_EW = (lambda_EW * red_EW ** 2) / (2 * (CYCLE_TIME - g))
    return wait_NS + wait_EW

# Find the optimal green time for NS direction
result = minimize_scalar(total_waiting_time, bounds=(1, CYCLE_TIME - 1), method='bounded')

optimal_g = result.x
min_wait = result.fun

print(f"✅ Optimal green time for NS: {optimal_g:.2f} seconds")
print(f"⏱️ Minimum total waiting time: {min_wait:.2f} vehicle-seconds")

green_times = np.linspace(1, CYCLE_TIME - 1, 100)
waiting_times = [total_waiting_time(g) for g in green_times]

plt.figure(figsize=(10, 6))
plt.plot(green_times, waiting_times, label='Total Waiting Time', color='blue')
plt.axvline(optimal_g, color='red', linestyle='--', label=f'Optimal g = {optimal_g:.2f}s')
plt.title('🚦 Traffic Signal Optimization: Total Waiting Time vs Green Time')
plt.xlabel('Green Time for NS (seconds)')
plt.ylabel('Total Waiting Time (vehicle-seconds)')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import linprog
import warnings
warnings.filterwarnings('ignore')

# Set up the problem parameters
class BasketballOptimizer:
    def __init__(self):
        # Player data: [Name, Offensive Rating, Defensive Rating, Fatigue Factor, Max Minutes]
        self.players_data = [
            ['LeBron James', 118, 110, 0.8, 42],
            ['Stephen Curry', 125, 105, 1.2, 40],
            ['Kawhi Leonard', 115, 115, 1.0, 35],
            ['Giannis Antetokounmpo', 120, 112, 0.9, 38],
            ['Kevin Durant', 122, 108, 1.1, 36],
            ['Luka Doncic', 116, 104, 1.3, 37],
            ['Jayson Tatum', 114, 109, 1.0, 39],
            ['Joel Embiid', 119, 113, 1.4, 32],
            ['Nikola Jokic', 121, 106, 0.7, 36],
            ['Jimmy Butler', 113, 111, 0.9, 35]
        ]
        
        # Create DataFrame
        self.df = pd.DataFrame(self.players_data, 
                              columns=['Player', 'Off_Rating', 'Def_Rating', 'Fatigue_Factor', 'Max_Minutes'])
        
        # Game parameters
        self.game_duration = 48  # minutes
        self.time_periods = 12   # 4-minute periods
        self.players_on_court = 5
        
    def calculate_performance_matrix(self):
        """Calculate performance for each player in each time period"""
        n_players = len(self.df)
        performance_matrix = np.zeros((n_players, self.time_periods))
        
        for i in range(n_players):
            base_performance = (self.df.iloc[i]['Off_Rating'] + self.df.iloc[i]['Def_Rating']) / 2
            fatigue_factor = self.df.iloc[i]['Fatigue_Factor']
            
            for t in range(self.time_periods):
                # Performance decreases with fatigue over time
                fatigue_penalty = fatigue_factor * (t + 1) * 0.5
                performance_matrix[i, t] = max(base_performance - fatigue_penalty, 50)
                
        return performance_matrix
    
    def solve_optimization(self):
        """Solve the player rotation optimization problem"""
        n_players = len(self.df)
        n_periods = self.time_periods
        
        # Get performance matrix
        performance_matrix = self.calculate_performance_matrix()
        
        # Create decision variables: x[i,t] for player i in period t
        # We'll use a simple greedy approach for demonstration
        rotation_schedule = np.zeros((n_players, n_periods))
        player_minutes = np.zeros(n_players)
        
        # For each time period, select top 5 available players
        for t in range(n_periods):
            # Calculate current performance considering fatigue and availability
            current_scores = []
            for i in range(n_players):
                if player_minutes[i] < self.df.iloc[i]['Max_Minutes']:
                    current_scores.append((performance_matrix[i, t], i))
                else:
                    current_scores.append((0, i))
            
            # Sort by performance and select top 5
            current_scores.sort(reverse=True)
            selected_players = [idx for _, idx in current_scores[:self.players_on_court]]
            
            # Update rotation schedule and minutes
            for player_idx in selected_players:
                rotation_schedule[player_idx, t] = 1
                player_minutes[player_idx] += 4  # 4 minutes per period
        
        return rotation_schedule, performance_matrix
    
    def analyze_results(self, rotation_schedule, performance_matrix):
        """Analyze optimization results"""
        results = {}
        
        # Calculate total minutes per player
        total_minutes = np.sum(rotation_schedule, axis=1) * 4
        
        # Calculate total performance
        total_performance = np.sum(rotation_schedule * performance_matrix)
        
        # Create results dataframe
        results_df = self.df.copy()
        results_df['Minutes_Played'] = total_minutes
        results_df['Utilization'] = (total_minutes / results_df['Max_Minutes'] * 100).round(1)
        
        # Calculate performance contribution
        performance_contribution = np.sum(rotation_schedule * performance_matrix, axis=1)
        results_df['Performance_Contribution'] = performance_contribution.round(1)
        
        results['summary'] = results_df
        results['total_performance'] = total_performance
        results['rotation_matrix'] = rotation_schedule
        
        return results
    
    def create_visualizations(self, results):
        """Create comprehensive visualizations"""
        fig, axes = plt.subplots(2, 2, figsize=(16, 12))
        
        # 1. Player Minutes Distribution
        ax1 = axes[0, 0]
        players = results['summary']['Player']
        minutes = results['summary']['Minutes_Played']
        max_minutes = results['summary']['Max_Minutes']
        
        x_pos = np.arange(len(players))
        bars1 = ax1.bar(x_pos, minutes, alpha=0.7, label='Minutes Played', color='skyblue')
        bars2 = ax1.bar(x_pos, max_minutes, alpha=0.3, label='Max Minutes', color='red')
        
        ax1.set_xlabel('Players')
        ax1.set_ylabel('Minutes')
        ax1.set_title('Player Minutes: Played vs Maximum Available')
        ax1.set_xticks(x_pos)
        ax1.set_xticklabels([name.split()[-1] for name in players], rotation=45)
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        # 2. Performance Contribution
        ax2 = axes[0, 1]
        performance = results['summary']['Performance_Contribution']
        bars = ax2.bar(range(len(players)), performance, color='lightgreen', alpha=0.7)
        
        ax2.set_xlabel('Players')
        ax2.set_ylabel('Performance Contribution')
        ax2.set_title('Individual Player Performance Contribution')
        ax2.set_xticks(range(len(players)))
        ax2.set_xticklabels([name.split()[-1] for name in players], rotation=45)
        ax2.grid(True, alpha=0.3)
        
        # Add value labels on bars
        for bar, value in zip(bars, performance):
            ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1,
                    f'{value:.0f}', ha='center', va='bottom')
        
        # 3. Rotation Schedule Heatmap
        ax3 = axes[1, 0]
        rotation_matrix = results['rotation_matrix']
        
        # Create heatmap
        im = ax3.imshow(rotation_matrix, cmap='YlOrRd', aspect='auto')
        ax3.set_xlabel('Time Periods (4-min each)')
        ax3.set_ylabel('Players')
        ax3.set_title('Player Rotation Schedule')
        ax3.set_xticks(range(self.time_periods))
        ax3.set_xticklabels([f'P{i+1}' for i in range(self.time_periods)])
        ax3.set_yticks(range(len(players)))
        ax3.set_yticklabels([name.split()[-1] for name in players])
        
        # Add colorbar
        plt.colorbar(im, ax=ax3, label='Playing (1) / Resting (0)')
        
        # 4. Utilization Rate
        ax4 = axes[1, 1]
        utilization = results['summary']['Utilization']
        colors = ['red' if x > 90 else 'orange' if x > 75 else 'green' for x in utilization]
        
        bars = ax4.bar(range(len(players)), utilization, color=colors, alpha=0.7)
        ax4.axhline(y=100, color='red', linestyle='--', alpha=0.5, label='Max Capacity')
        ax4.axhline(y=75, color='orange', linestyle='--', alpha=0.5, label='High Utilization')
        
        ax4.set_xlabel('Players')
        ax4.set_ylabel('Utilization Rate (%)')
        ax4.set_title('Player Utilization Rate')
        ax4.set_xticks(range(len(players)))
        ax4.set_xticklabels([name.split()[-1] for name in players], rotation=45)
        ax4.legend()
        ax4.grid(True, alpha=0.3)
        
        # Add percentage labels
        for bar, value in zip(bars, utilization):
            ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1,
                    f'{value:.1f}%', ha='center', va='bottom')
        
        plt.tight_layout()
        plt.show()
        
        return fig
    
    def generate_insights(self, results):
        """Generate strategic insights from the optimization"""
        summary = results['summary']
        
        print("=== BASKETBALL TEAM ROTATION OPTIMIZATION RESULTS ===\n")
        
        print("📊 PERFORMANCE SUMMARY:")
        print(f"Total Team Performance Score: {results['total_performance']:.1f}")
        print(f"Average Player Utilization: {summary['Utilization'].mean():.1f}%")
        print(f"Players at >90% Utilization: {sum(summary['Utilization'] > 90)}")
        
        print("\n🏀 TOP PERFORMERS:")
        top_performers = summary.nlargest(3, 'Performance_Contribution')
        for idx, (_, player) in enumerate(top_performers.iterrows()):
            print(f"{idx+1}. {player['Player']}: {player['Performance_Contribution']:.1f} points "
                  f"({player['Minutes_Played']} min, {player['Utilization']:.1f}% utilization)")
        
        print("\n⚠️  STRATEGIC RECOMMENDATIONS:")
        
        # High utilization players
        high_util = summary[summary['Utilization'] > 85]
        if len(high_util) > 0:
            print("• High Utilization Alert:")
            for _, player in high_util.iterrows():
                print(f"  - {player['Player']}: {player['Utilization']:.1f}% utilization - consider rest periods")
        
        # Underutilized players
        low_util = summary[summary['Utilization'] < 60]
        if len(low_util) > 0:
            print("• Underutilized Players:")
            for _, player in low_util.iterrows():
                print(f"  - {player['Player']}: {player['Utilization']:.1f}% utilization - opportunity for more minutes")
        
        # Performance efficiency
        summary['Efficiency'] = summary['Performance_Contribution'] / summary['Minutes_Played']
        top_efficiency = summary.nlargest(3, 'Efficiency')
        print("\n💡 MOST EFFICIENT PLAYERS (Performance per Minute):")
        for idx, (_, player) in enumerate(top_efficiency.iterrows()):
            print(f"{idx+1}. {player['Player']}: {player['Efficiency']:.2f} points/minute")
        
        return summary

# Execute the optimization
def main():
    # Initialize optimizer
    optimizer = BasketballOptimizer()
    
    print("🏀 Basketball Team Rotation Optimization")
    print("=" * 50)
    
    # Display initial team data
    print("\n📋 TEAM ROSTER:")
    print(optimizer.df.to_string(index=False))
    
    # Solve optimization
    print("\n🔄 Solving optimization problem...")
    rotation_schedule, performance_matrix = optimizer.solve_optimization()
    
    # Analyze results
    results = optimizer.analyze_results(rotation_schedule, performance_matrix)
    
    # Generate insights
    insights = optimizer.generate_insights(results)
    
    # Create visualizations
    print("\n📈 Generating visualizations...")
    fig = optimizer.create_visualizations(results)
    
    # Additional analysis: Performance over time
    plt.figure(figsize=(14, 8))
    
    # Plot team performance over time periods
    plt.subplot(2, 1, 1)
    period_performance = np.sum(rotation_schedule * performance_matrix, axis=0)
    periods = [f'Q{(i//3)+1}-{(i%3)*4+4}min' for i in range(optimizer.time_periods)]
    
    plt.plot(periods, period_performance, marker='o', linewidth=2, markersize=8, color='blue')
    plt.title('Team Performance Throughout the Game', fontsize=14, fontweight='bold')
    plt.xlabel('Game Period')
    plt.ylabel('Team Performance Score')
    plt.grid(True, alpha=0.3)
    plt.xticks(rotation=45)
    
    # Plot active players per period
    plt.subplot(2, 1, 2)
    active_players = np.sum(rotation_schedule, axis=0)
    plt.bar(periods, active_players, alpha=0.7, color='lightcoral')
    plt.title('Number of Active Players per Period', fontsize=14, fontweight='bold')
    plt.xlabel('Game Period')
    plt.ylabel('Active Players')
    plt.ylim(0, 6)
    plt.grid(True, alpha=0.3)
    plt.xticks(rotation=45)
    
    # Add horizontal line at 5 players
    plt.axhline(y=5, color='red', linestyle='--', alpha=0.7, label='Required Players')
    plt.legend()
    
    plt.tight_layout()
    plt.show()
    
    return optimizer, results

# Run the main function
if __name__ == "__main__":
    optimizer, results = main()

fatigue_penalty = fatigue_factor * (t + 1) * 0.5
performance_matrix[i, t] = max(base_performance - fatigue_penalty, 50)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from datetime import datetime, timedelta

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

# Generate synthetic time series data
dates = [datetime(2025, 6, 1) + timedelta(hours=i) for i in range(168)]  # One week of hourly data
cpu_usage = np.random.normal(loc=50, scale=10, size=len(dates))  # Normal CPU usage ~ N(50, 10)
cpu_usage[50:52] = [90, 95]  # Introduce spike anomalies
cpu_usage[120:122] = [10, 5]  # Introduce drop anomalies

# Create DataFrame
df = pd.DataFrame({'timestamp': dates, 'cpu_usage': cpu_usage})

# Define anomaly detection function
def detect_anomalies(df, window=24, z_threshold=2.5):
    # Calculate rolling mean and standard deviation
    df['rolling_mean'] = df['cpu_usage'].rolling(window=window, center=True).mean()
    df['rolling_std'] = df['cpu_usage'].rolling(window=window, center=True).std()
    
    # Calculate z-score
    df['z_score'] = (df['cpu_usage'] - df['rolling_mean']) / df['rolling_std']
    
    # Identify anomalies
    df['anomaly'] = df['z_score'].abs() > z_threshold
    
    return df

# Apply anomaly detection
df = detect_anomalies(df)

# Plot the results
plt.figure(figsize=(12, 6))
plt.plot(df['timestamp'], df['cpu_usage'], label='CPU Usage', color='blue')
plt.plot(df['timestamp'], df['rolling_mean'], label='Rolling Mean', color='orange', linestyle='--')
plt.fill_between(df['timestamp'], 
                 df['rolling_mean'] - 2.5 * df['rolling_std'], 
                 df['rolling_mean'] + 2.5 * df['rolling_std'], 
                 color='gray', alpha=0.2, label='Normal Range (±2.5σ)')
plt.scatter(df[df['anomaly']]['timestamp'], df[df['anomaly']]['cpu_usage'], 
            color='red', label='Anomalies', s=100, marker='o')
plt.title('Server CPU Usage with Anomaly Detection')
plt.xlabel('Timestamp')
plt.ylabel('CPU Usage (%)')
plt.legend()
plt.grid(True)
plt.xticks(rotation=45)
plt.tight_layout()

# Save the plot
plt.savefig('cpu_usage_anomaly_detection.png')