# Import necessary libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import load_wine
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from mpl_toolkits.mplot3d import Axes3D

# Set style for better visualizations
plt.style.use('seaborn-v0_8')
sns.set_palette("husl")

# Load the Wine dataset
wine_data = load_wine()
X = wine_data.data
y = wine_data.target
feature_names = wine_data.feature_names
target_names = wine_data.target_names

print("Dataset Information:")
print(f"Number of samples: {X.shape[0]}")
print(f"Number of features: {X.shape[1]}")
print(f"Target classes: {target_names}")
print(f"Class distribution: {np.bincount(y)}")

# Convert to DataFrame for easier handling
df = pd.DataFrame(X, columns=feature_names)
df['target'] = y
df['target_name'] = [target_names[i] for i in y]

print("\nFirst 5 rows of the dataset:")
print(df.head())

# Display basic statistics
print("\nDataset Statistics:")
print(df.describe())

# Step 1: Data Standardization
print("\n" + "="*50)
print("STEP 1: DATA STANDARDIZATION")
print("="*50)

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

print("Original data statistics:")
print(f"Mean: {np.mean(X, axis=0)[:3]}... (showing first 3 features)")
print(f"Std:  {np.std(X, axis=0)[:3]}... (showing first 3 features)")

print("\nStandardized data statistics:")
print(f"Mean: {np.mean(X_scaled, axis=0)[:3]}... (should be ~0)")
print(f"Std:  {np.std(X_scaled, axis=0)[:3]}... (should be ~1)")

# Step 2: Apply PCA
print("\n" + "="*50)
print("STEP 2: PRINCIPAL COMPONENT ANALYSIS")
print("="*50)

# Perform PCA with all components first to analyze explained variance
pca_full = PCA()
X_pca_full = pca_full.fit_transform(X_scaled)

# Calculate cumulative explained variance
explained_variance_ratio = pca_full.explained_variance_ratio_
cumulative_variance = np.cumsum(explained_variance_ratio)

print("Explained Variance by each Principal Component:")
for i, (var, cum_var) in enumerate(zip(explained_variance_ratio, cumulative_variance)):
    if i < 10:  # Show first 10 components
        print(f"PC{i+1}: {var:.4f} ({var*100:.2f}%) - Cumulative: {cum_var:.4f} ({cum_var*100:.2f}%)")

# Determine optimal number of components (e.g., 95% variance)
n_components_95 = np.argmax(cumulative_variance >= 0.95) + 1
n_components_90 = np.argmax(cumulative_variance >= 0.90) + 1

print(f"\nComponents needed for 90% variance: {n_components_90}")
print(f"Components needed for 95% variance: {n_components_95}")

# Apply PCA with reduced dimensions (let's use 3 components for visualization)
pca_3d = PCA(n_components=3)
X_pca_3d = pca_3d.fit_transform(X_scaled)

print(f"\nExplained variance with 3 components: {sum(pca_3d.explained_variance_ratio_):.4f} ({sum(pca_3d.explained_variance_ratio_)*100:.2f}%)")

# Step 3: Analyze Principal Components
print("\n" + "="*50)
print("STEP 3: PRINCIPAL COMPONENTS ANALYSIS")
print("="*50)

# Create DataFrame for PC loadings
pc_loadings = pd.DataFrame(
    pca_3d.components_.T,
    columns=['PC1', 'PC2', 'PC3'],
    index=feature_names
)

print("Principal Component Loadings (Top 5 features for each PC):")
for pc in ['PC1', 'PC2', 'PC3']:
    print(f"\n{pc} - Top contributors:")
    top_features = pc_loadings[pc].abs().sort_values(ascending=False).head(5)
    for feature, loading in top_features.items():
        direction = "+" if pc_loadings.loc[feature, pc] > 0 else "-"
        print(f"  {direction} {feature}: {abs(loading):.3f}")

# Step 4: Visualization
print("\n" + "="*50)
print("STEP 4: CREATING VISUALIZATIONS")
print("="*50)

# Create comprehensive visualization
fig = plt.figure(figsize=(20, 15))

# 1. Explained Variance Plot
ax1 = plt.subplot(2, 3, 1)
plt.bar(range(1, 14), explained_variance_ratio, alpha=0.7, label='Individual')
plt.plot(range(1, 14), cumulative_variance, 'ro-', label='Cumulative')
plt.xlabel('Principal Component')
plt.ylabel('Explained Variance Ratio')
plt.title('Explained Variance by Principal Components')
plt.legend()
plt.grid(True, alpha=0.3)

# Add annotation for 95% threshold
plt.axhline(y=0.95, color='red', linestyle='--', alpha=0.7)
plt.text(8, 0.96, '95% threshold', fontsize=10)

# 2. 2D PCA Scatter Plot (PC1 vs PC2)
ax2 = plt.subplot(2, 3, 2)
colors = ['red', 'blue', 'green']
for i, (target_class, color) in enumerate(zip(target_names, colors)):
    mask = y == i
    plt.scatter(X_pca_3d[mask, 0], X_pca_3d[mask, 1], 
               c=color, label=target_class, alpha=0.7, s=50)

plt.xlabel(f'PC1 ({pca_3d.explained_variance_ratio_[0]:.2%} variance)')
plt.ylabel(f'PC2 ({pca_3d.explained_variance_ratio_[1]:.2%} variance)')
plt.title('PCA: First Two Principal Components')
plt.legend()
plt.grid(True, alpha=0.3)

# 3. 3D PCA Scatter Plot
ax3 = plt.subplot(2, 3, 3, projection='3d')
for i, (target_class, color) in enumerate(zip(target_names, colors)):
    mask = y == i
    ax3.scatter(X_pca_3d[mask, 0], X_pca_3d[mask, 1], X_pca_3d[mask, 2],
               c=color, label=target_class, alpha=0.7, s=30)

ax3.set_xlabel(f'PC1 ({pca_3d.explained_variance_ratio_[0]:.2%})')
ax3.set_ylabel(f'PC2 ({pca_3d.explained_variance_ratio_[1]:.2%})')
ax3.set_zlabel(f'PC3 ({pca_3d.explained_variance_ratio_[2]:.2%})')
ax3.set_title('3D PCA Visualization')
ax3.legend()

# 4. Feature Contributions Heatmap
ax4 = plt.subplot(2, 3, 4)
sns.heatmap(pc_loadings.T, annot=False, cmap='RdBu_r', center=0,
            xticklabels=[name.replace('_', ' ').title() for name in feature_names],
            yticklabels=['PC1', 'PC2', 'PC3'])
plt.title('Principal Component Loadings Heatmap')
plt.xlabel('Features')
plt.ylabel('Principal Components')
plt.xticks(rotation=45, ha='right')

# 5. PC1 vs PC3
ax5 = plt.subplot(2, 3, 5)
for i, (target_class, color) in enumerate(zip(target_names, colors)):
    mask = y == i
    plt.scatter(X_pca_3d[mask, 0], X_pca_3d[mask, 2], 
               c=color, label=target_class, alpha=0.7, s=50)

plt.xlabel(f'PC1 ({pca_3d.explained_variance_ratio_[0]:.2%} variance)')
plt.ylabel(f'PC3 ({pca_3d.explained_variance_ratio_[2]:.2%} variance)')
plt.title('PCA: PC1 vs PC3')
plt.legend()
plt.grid(True, alpha=0.3)

# 6. PC2 vs PC3
ax6 = plt.subplot(2, 3, 6)
for i, (target_class, color) in enumerate(zip(target_names, colors)):
    mask = y == i
    plt.scatter(X_pca_3d[mask, 1], X_pca_3d[mask, 2], 
               c=color, label=target_class, alpha=0.7, s=50)

plt.xlabel(f'PC2 ({pca_3d.explained_variance_ratio_[1]:.2%} variance)')
plt.ylabel(f'PC3 ({pca_3d.explained_variance_ratio_[2]:.2%} variance)')
plt.title('PCA: PC2 vs PC3')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Additional Analysis: Biplot
print("\nCreating Biplot for detailed feature analysis...")
fig, ax = plt.subplots(figsize=(12, 10))

# Plot data points
for i, (target_class, color) in enumerate(zip(target_names, colors)):
    mask = y == i
    ax.scatter(X_pca_3d[mask, 0], X_pca_3d[mask, 1], 
              c=color, label=target_class, alpha=0.6, s=50)

# Plot feature vectors
feature_vectors = pca_3d.components_[:2].T * 3  # Scale for visibility
for i, (feature, vector) in enumerate(zip(feature_names, feature_vectors)):
    ax.arrow(0, 0, vector[0], vector[1], head_width=0.1, head_length=0.1, 
             fc='black', ec='black', alpha=0.7)
    ax.text(vector[0]*1.1, vector[1]*1.1, feature.replace('_', ' ').title(), 
            fontsize=8, ha='center', va='center')

ax.set_xlabel(f'PC1 ({pca_3d.explained_variance_ratio_[0]:.2%} variance)')
ax.set_ylabel(f'PC2 ({pca_3d.explained_variance_ratio_[1]:.2%} variance)')
ax.set_title('PCA Biplot: Data Points and Feature Contributions')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)

plt.tight_layout()
plt.show()

# Step 5: Reconstruction Analysis
print("\n" + "="*50)
print("STEP 5: RECONSTRUCTION ANALYSIS")
print("="*50)

# Reconstruct data using different numbers of components
n_components_list = [1, 2, 3, 5, 7, 10, 13]
reconstruction_errors = []

for n_comp in n_components_list:
    pca_temp = PCA(n_components=n_comp)
    X_transformed = pca_temp.fit_transform(X_scaled)
    X_reconstructed = pca_temp.inverse_transform(X_transformed)
    
    # Calculate reconstruction error (MSE)
    mse = np.mean((X_scaled - X_reconstructed) ** 2)
    reconstruction_errors.append(mse)
    
    if n_comp <= 5:
        print(f"Components: {n_comp:2d}, Explained Variance: {sum(pca_temp.explained_variance_ratio_):.4f}, Reconstruction Error: {mse:.6f}")

# Plot reconstruction error
plt.figure(figsize=(10, 6))
plt.plot(n_components_list, reconstruction_errors, 'bo-', linewidth=2, markersize=8)
plt.xlabel('Number of Principal Components')
plt.ylabel('Reconstruction Error (MSE)')
plt.title('Reconstruction Error vs Number of Components')
plt.grid(True, alpha=0.3)
plt.yscale('log')

# Add annotations
for i, (n_comp, error) in enumerate(zip(n_components_list, reconstruction_errors)):
    if n_comp in [2, 3, 5]:
        plt.annotate(f'{error:.4f}', (n_comp, error), 
                    textcoords="offset points", xytext=(0,10), ha='center')

plt.tight_layout()
plt.show()

# Final Summary
print("\n" + "="*50)
print("ANALYSIS SUMMARY")
print("="*50)
print(f"• Original dataset: {X.shape[0]} samples × {X.shape[1]} features")
print(f"• First 3 PCs explain {sum(pca_3d.explained_variance_ratio_)*100:.1f}% of total variance")
print(f"• PC1 captures {pca_3d.explained_variance_ratio_[0]*100:.1f}% (likely represents overall wine quality/intensity)")
print(f"• PC2 captures {pca_3d.explained_variance_ratio_[1]*100:.1f}% (likely represents color-related properties)")
print(f"• PC3 captures {pca_3d.explained_variance_ratio_[2]*100:.1f}% (likely represents acidity/pH balance)")
print(f"• The three wine classes show good separation in the PC space")
print(f"• Dimensionality reduction from 13D to 3D maintains {sum(pca_3d.explained_variance_ratio_)*100:.1f}% of information")

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler

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

# Generate synthetic dataset
n_samples = 100
house_size = np.random.normal(2000, 500, n_samples)  # Size in sq ft
bedrooms = house_size / 500 + np.random.normal(0, 0.5, n_samples)  # Correlated with size
age = np.random.normal(20, 10, n_samples)  # Age in years
X = np.column_stack((house_size, bedrooms, age))

# True coefficients: price = 100 * size + 50 * bedrooms - 20 * age + noise
true_coefficients = [100, 50, -20]
y = X @ true_coefficients + np.random.normal(0, 10000, n_samples)

# Scale features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Fit Linear Regression
lr = LinearRegression()
lr.fit(X_scaled, y)

# Fit Ridge Regression with different alpha values
alphas = [0.1, 1.0, 10.0]
ridge_models = {alpha: Ridge(alpha=alpha).fit(X_scaled, y) for alpha in alphas}

# Predictions
y_pred_lr = lr.predict(X_scaled)
y_pred_ridge = {alpha: model.predict(X_scaled) for alpha, model in ridge_models.items()}

# Calculate MSE
mse_lr = mean_squared_error(y, y_pred_lr)
mse_ridge = {alpha: mean_squared_error(y, y_pred_ridge[alpha]) for alpha in alphas}

# Plotting
plt.figure(figsize=(12, 6))

# Plot 1: Actual vs Predicted for Linear Regression
plt.subplot(1, 2, 1)
plt.scatter(y, y_pred_lr, color='blue', alpha=0.5, label='Predicted')
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'r--', lw=2, label='Ideal')
plt.xlabel('Actual Price ($)')
plt.ylabel('Predicted Price ($)')
plt.title(f'Linear Regression\nMSE: {mse_lr:.2f}')
plt.legend()

# Plot 2: Actual vs Predicted for Ridge Regression (alpha=1.0)
plt.subplot(1, 2, 2)
plt.scatter(y, y_pred_ridge[1.0], color='green', alpha=0.5, label='Predicted')
plt.plot([y.min(), y.max()], [y.min(), y.max()], 'r--', lw=2, label='Ideal')
plt.xlabel('Actual Price ($)')
plt.ylabel('Predicted Price ($)')
plt.title(f'Ridge Regression (α=1.0)\nMSE: {mse_ridge[1.0]:.2f}')
plt.legend()

plt.tight_layout()
plt.savefig('ridge_regression_plot.png')
plt.show()

# Print coefficients
print("Linear Regression Coefficients:", lr.coef_)
print("Ridge Regression Coefficients:")
for alpha, model in ridge_models.items():
    print(f"Alpha = {alpha}: {model.coef_}")

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from itertools import permutations
import pandas as pd
from scipy.optimize import linear_sum_assignment
import networkx as nx
import math

class QAPSolver:
    def __init__(self, flow_matrix, distance_matrix, facility_names, location_names):
        """
        Initialize QAP Solver
        
        Parameters:
        flow_matrix: Matrix representing flow between facilities
        distance_matrix: Matrix representing distances between locations  
        facility_names: Names of facilities
        location_names: Names of locations
        """
        self.flow_matrix = np.array(flow_matrix)
        self.distance_matrix = np.array(distance_matrix)
        self.facility_names = facility_names
        self.location_names = location_names
        self.n = len(facility_names)
        self.best_assignment = None
        self.best_cost = float('inf')
        self.all_solutions = []
        
    def calculate_cost(self, assignment):
        """
        Calculate total cost for a given assignment
        
        Parameters:
        assignment: List where assignment[i] is the location assigned to facility i
        
        Returns:
        Total cost of the assignment
        """
        total_cost = 0
        for i in range(self.n):
            for j in range(self.n):
                # Flow between facilities i and j
                flow = self.flow_matrix[i][j]
                # Distance between assigned locations
                distance = self.distance_matrix[assignment[i]][assignment[j]]
                total_cost += flow * distance
        return total_cost
    
    def solve_exhaustive(self):
        """
        Solve QAP using exhaustive search (suitable for small instances)
        """
        print("Solving QAP using exhaustive search...")
        
        # Generate all possible permutations
        for perm in permutations(range(self.n)):
            cost = self.calculate_cost(perm)
            self.all_solutions.append((list(perm), cost))
            
            if cost < self.best_cost:
                self.best_cost = cost
                self.best_assignment = list(perm)
        
        # Sort solutions by cost
        self.all_solutions.sort(key=lambda x: x[1])
        
        print(f"Optimal cost: {self.best_cost}")
        print(f"Optimal assignment: {self.best_assignment}")
        return self.best_assignment, self.best_cost
    
    def solve_greedy_heuristic(self):
        """
        Solve QAP using a greedy heuristic approach
        """
        print("Solving QAP using greedy heuristic...")
        
        # Calculate assignment priorities based on flow*distance products
        assignment_costs = np.zeros((self.n, self.n))
        
        for i in range(self.n):
            for j in range(self.n):
                # Calculate cost if facility i is assigned to location j
                cost = 0
                for k in range(self.n):
                    for l in range(self.n):
                        if k != i:  # Other facilities
                            # Estimate cost assuming worst case for other assignments  
                            flow = self.flow_matrix[i][k] + self.flow_matrix[k][i]
                            distance = self.distance_matrix[j][l]
                            cost += flow * distance
                assignment_costs[i][j] = cost
        
        # Use Hungarian algorithm for approximation
        row_ind, col_ind = linear_sum_assignment(assignment_costs)
        heuristic_assignment = col_ind.tolist()
        heuristic_cost = self.calculate_cost(heuristic_assignment)
        
        print(f"Heuristic cost: {heuristic_cost}")
        print(f"Heuristic assignment: {heuristic_assignment}")
        
        return heuristic_assignment, heuristic_cost
    
    def display_results(self):
        """Display the results in a formatted way"""
        print("\n" + "="*60)
        print("QUADRATIC ASSIGNMENT PROBLEM RESULTS")
        print("="*60)
        
        print(f"\nProblem Size: {self.n} facilities, {self.n} locations")
        print(f"Total possible assignments: {math.factorial(self.n)}")
        
        print(f"\nOptimal Assignment:")
        for i, loc in enumerate(self.best_assignment):
            print(f"  {self.facility_names[i]} → {self.location_names[loc]}")
        
        print(f"\nOptimal Cost: {self.best_cost}")
        
        # Show top 5 best and worst solutions
        print(f"\nTop 5 Best Solutions:")
        for i, (assignment, cost) in enumerate(self.all_solutions[:5]):
            assignment_str = " → ".join([f"{self.facility_names[j]}:{self.location_names[assignment[j]]}" 
                                       for j in range(self.n)])
            print(f"  {i+1}. Cost: {cost:6.1f} | {assignment_str}")
            
        print(f"\nTop 5 Worst Solutions:")
        for i, (assignment, cost) in enumerate(self.all_solutions[-5:]):
            assignment_str = " → ".join([f"{self.facility_names[j]}:{self.location_names[assignment[j]]}" 
                                       for j in range(self.n)])
            print(f"  {i+1}. Cost: {cost:6.1f} | {assignment_str}")

    def visualize_results(self):
        """Create comprehensive visualizations of the QAP solution"""
        
        # Set up the plotting style
        plt.style.use('default')
        sns.set_palette("husl")
        
        # Create figure with subplots
        fig = plt.figure(figsize=(20, 15))
        
        # 1. Flow Matrix Heatmap
        ax1 = plt.subplot(2, 4, 1)
        sns.heatmap(self.flow_matrix, annot=True, fmt='g', cmap='YlOrRd',
                   xticklabels=self.facility_names, yticklabels=self.facility_names,
                   cbar_kws={'label': 'Flow Units'})
        plt.title('Flow Matrix Between Facilities', fontsize=12, fontweight='bold')
        plt.xlabel('To Facility')
        plt.ylabel('From Facility')
        
        # 2. Distance Matrix Heatmap  
        ax2 = plt.subplot(2, 4, 2)
        sns.heatmap(self.distance_matrix, annot=True, fmt='g', cmap='Blues',
                   xticklabels=self.location_names, yticklabels=self.location_names,
                   cbar_kws={'label': 'Distance Units'})
        plt.title('Distance Matrix Between Locations', fontsize=12, fontweight='bold')
        plt.xlabel('To Location')
        plt.ylabel('From Location')
        
        # 3. Cost Distribution
        ax3 = plt.subplot(2, 4, 3)
        costs = [sol[1] for sol in self.all_solutions]
        plt.hist(costs, bins=min(20, len(costs)//2), alpha=0.7, color='skyblue', edgecolor='black')
        plt.axvline(self.best_cost, color='red', linestyle='--', linewidth=2, label=f'Optimal: {self.best_cost}')
        plt.xlabel('Total Cost')
        plt.ylabel('Number of Solutions')
        plt.title('Distribution of Solution Costs', fontsize=12, fontweight='bold')
        plt.legend()
        plt.grid(True, alpha=0.3)
        
        # 4. Assignment Network Graph
        ax4 = plt.subplot(2, 4, 4)
        G = nx.Graph()
        
        # Add nodes for facilities and locations
        facility_nodes = [f"F_{i}_{self.facility_names[i][:3]}" for i in range(self.n)]
        location_nodes = [f"L_{i}_{self.location_names[i][:3]}" for i in range(self.n)]
        
        G.add_nodes_from(facility_nodes)
        G.add_nodes_from(location_nodes)
        
        # Add edges for optimal assignment
        for i, loc in enumerate(self.best_assignment):
            G.add_edge(facility_nodes[i], location_nodes[loc])
        
        # Position nodes
        pos = {}
        for i, node in enumerate(facility_nodes):
            pos[node] = (0, i)
        for i, node in enumerate(location_nodes):
            pos[node] = (2, i)
            
        nx.draw(G, pos, ax=ax4, with_labels=True, node_color=['lightcoral']*self.n + ['lightblue']*self.n,
                node_size=1500, font_size=8, font_weight='bold')
        plt.title('Optimal Assignment Network', fontsize=12, fontweight='bold')
        
        # 5. Cost Contribution Matrix
        ax5 = plt.subplot(2, 4, 5)
        cost_matrix = np.zeros((self.n, self.n))
        for i in range(self.n):
            for j in range(self.n):
                flow = self.flow_matrix[i][j]
                distance = self.distance_matrix[self.best_assignment[i]][self.best_assignment[j]]
                cost_matrix[i][j] = flow * distance
                
        sns.heatmap(cost_matrix, annot=True, fmt='.1f', cmap='Reds',
                   xticklabels=self.facility_names, yticklabels=self.facility_names,
                   cbar_kws={'label': 'Cost Contribution'})
        plt.title('Cost Contribution Matrix\n(Optimal Assignment)', fontsize=12, fontweight='bold')
        plt.xlabel('To Facility')
        plt.ylabel('From Facility')
        
        # 6. Solution Quality Comparison
        ax6 = plt.subplot(2, 4, 6)
        x = range(len(self.all_solutions))
        y = [sol[1] for sol in self.all_solutions]
        plt.plot(x, y, 'b-', alpha=0.7, linewidth=1, label='All Solutions')
        plt.scatter(0, self.best_cost, color='red', s=100, zorder=5, label=f'Optimal: {self.best_cost}')
        plt.xlabel('Solution Rank')
        plt.ylabel('Total Cost')
        plt.title('Solution Quality Ranking', fontsize=12, fontweight='bold')
        plt.legend()
        plt.grid(True, alpha=0.3)
        
        # 7. Facility-Location Assignment Bar Chart
        ax7 = plt.subplot(2, 4, 7)
        y_pos = np.arange(self.n)
        colors = sns.color_palette("husl", self.n)
        
        bars = plt.barh(y_pos, [self.best_assignment[i] for i in range(self.n)], color=colors)
        plt.yticks(y_pos, self.facility_names)
        plt.xlabel('Assigned Location Index')
        plt.title('Optimal Facility Assignments', fontsize=12, fontweight='bold')
        plt.grid(True, alpha=0.3, axis='x')
        
        # Add location names as text
        for i, bar in enumerate(bars):
            width = bar.get_width()
            location_name = self.location_names[int(width)]
            plt.text(width + 0.05, bar.get_y() + bar.get_height()/2, 
                    location_name, ha='left', va='center', fontweight='bold')
        
        # 8. Performance Metrics
        ax8 = plt.subplot(2, 4, 8)
        ax8.axis('off')
        
        # Calculate some performance metrics
        total_flow = np.sum(self.flow_matrix)
        total_distance = np.sum(self.distance_matrix)
        avg_cost = np.mean(costs)
        std_cost = np.std(costs)
        worst_cost = max(costs)
        gap = ((worst_cost - self.best_cost) / self.best_cost) * 100
        
        metrics_text = f"""
        PERFORMANCE METRICS
        ─────────────────────
        Total Flow: {total_flow:.1f}
        Total Distance: {total_distance:.1f}
        
        Optimal Cost: {self.best_cost:.1f}
        Average Cost: {avg_cost:.1f}
        Worst Cost: {worst_cost:.1f}
        Standard Deviation: {std_cost:.1f}
        
        Optimality Gap: {gap:.1f}%
        Solutions Evaluated: {len(self.all_solutions)}
        """
        
        plt.text(0.1, 0.9, metrics_text, transform=ax8.transAxes, fontsize=11,
                verticalalignment='top', fontfamily='monospace',
                bbox=dict(boxstyle="round,pad=0.3", facecolor="lightgray", alpha=0.8))
        
        plt.tight_layout()
        plt.suptitle('Quadratic Assignment Problem: Complete Analysis', 
                    fontsize=16, fontweight='bold', y=0.98)
        plt.show()
        
        # Create a detailed assignment comparison table
        self.create_assignment_table()
    
    def create_assignment_table(self):
        """Create a detailed table comparing different assignments"""
        print("\n" + "="*80)
        print("DETAILED ASSIGNMENT ANALYSIS")
        print("="*80)
        
        # Create DataFrame for analysis
        analysis_data = []
        
        for rank, (assignment, cost) in enumerate(self.all_solutions[:10]):  # Top 10 solutions
            row = {
                'Rank': rank + 1,
                'Total_Cost': cost,
                'Cost_Gap_%': ((cost - self.best_cost) / self.best_cost) * 100
            }
            
            # Add individual assignments
            for i, loc in enumerate(assignment):
                row[f'{self.facility_names[i]}'] = self.location_names[loc]
            
            analysis_data.append(row)
            
        df = pd.DataFrame(analysis_data)
        print(df.to_string(index=False))
        
        # Flow-Distance Analysis
        print(f"\n\nFLOW-DISTANCE INTERACTION ANALYSIS (Optimal Solution)")
        print("-" * 60)
        
        for i in range(self.n):
            for j in range(self.n):
                if i != j:  # Skip diagonal
                    flow = self.flow_matrix[i][j]
                    if flow > 0:  # Only show non-zero flows
                        loc_i = self.best_assignment[i]
                        loc_j = self.best_assignment[j]
                        distance = self.distance_matrix[loc_i][loc_j]
                        cost_contribution = flow * distance
                        
                        print(f"{self.facility_names[i]:>12} → {self.facility_names[j]:<12} | "
                              f"Flow: {flow:>4.1f} × Distance: {distance:>4.1f} = Cost: {cost_contribution:>6.1f}")

# Define the hospital department layout problem
def create_hospital_qap():
    """
    Create a hospital department layout QAP instance
    """
    
    # Flow matrix (patients/staff movement per day between departments)
    # Rows/Cols: ER, Surgery, Laboratory, Pharmacy
    flow_matrix = [
        [0,   25,  30,  40],  # ER to others
        [25,  0,   35,  20],  # Surgery to others  
        [30,  35,  0,   15],  # Laboratory to others
        [40,  20,  15,  0]    # Pharmacy to others
    ]
    
    # Distance matrix between locations (in meters)
    # Rows/Cols: Location A, B, C, D
    distance_matrix = [
        [0,   50,  80,  120], # Location A to others
        [50,  0,   60,  90],  # Location B to others
        [80,  60,  0,   40],  # Location C to others  
        [120, 90,  40,  0]    # Location D to others
    ]
    
    facility_names = ['Emergency Room', 'Surgery Dept', 'Laboratory', 'Pharmacy']
    location_names = ['Location A', 'Location B', 'Location C', 'Location D'] 
    
    return QAPSolver(flow_matrix, distance_matrix, facility_names, location_names)

# Main execution
if __name__ == "__main__":
    print("QUADRATIC ASSIGNMENT PROBLEM SOLVER")
    print("Hospital Department Layout Optimization")
    print("="*50)
    
    # Create and solve the hospital QAP instance
    qap = create_hospital_qap()
    
    # Solve using exhaustive search
    optimal_assignment, optimal_cost = qap.solve_exhaustive()
    
    # Solve using heuristic for comparison  
    heuristic_assignment, heuristic_cost = qap.solve_greedy_heuristic()
    
    # Display results
    qap.display_results()
    
    # Compare heuristic vs optimal
    if optimal_cost != heuristic_cost:
        gap = ((heuristic_cost - optimal_cost) / optimal_cost) * 100
        print(f"\nHeuristic vs Optimal Comparison:")
        print(f"Optimal Cost: {optimal_cost}")
        print(f"Heuristic Cost: {heuristic_cost}")
        print(f"Gap: {gap:.2f}%")
    else:
        print(f"\nHeuristic found the optimal solution!")
    
    # Create visualizations
    qap.visualize_results()

def calculate_cost(self, assignment):
    total_cost = 0
    for i in range(self.n):
        for j in range(self.n):
            flow = self.flow_matrix[i][j]
            distance = self.distance_matrix[assignment[i]][assignment[j]]
            total_cost += flow * distance
    return total_cost

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error
from sklearn.preprocessing import StandardScaler
from scipy import optimize
import warnings
warnings.filterwarnings('ignore')

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

class DemandForecaster:
    """Machine Learning based demand forecaster with uncertainty quantification"""
    
    def __init__(self):
        self.model = GradientBoostingRegressor(
            n_estimators=100,
            learning_rate=0.1,
            max_depth=4,
            random_state=42
        )
        self.scaler = StandardScaler()
        self.is_fitted = False
        
    def create_features(self, data, lookback=7):
        """Create time series features for ML model"""
        features = []
        targets = []
        
        for i in range(lookback, len(data)):
            # Lagged demand features
            lag_features = data[i-lookback:i].tolist()
            
            # Additional features
            day_of_week = i % 7
            week_of_year = (i // 7) % 52
            trend = i / len(data)  # Linear trend
            
            # Seasonal features
            seasonal = np.sin(2 * np.pi * i / 7)  # Weekly seasonality
            
            feature_vector = lag_features + [day_of_week, week_of_year, trend, seasonal]
            features.append(feature_vector)
            targets.append(data[i])
            
        return np.array(features), np.array(targets)
    
    def fit(self, demand_data):
        """Train the demand forecasting model"""
        X, y = self.create_features(demand_data)
        X_scaled = self.scaler.fit_transform(X)
        self.model.fit(X_scaled, y)
        self.is_fitted = True
        return self
    
    def predict_with_uncertainty(self, demand_data, horizon=30, n_simulations=100):
        """Predict demand with uncertainty bounds using quantile regression approach"""
        if not self.is_fitted:
            raise ValueError("Model must be fitted before prediction")
            
        X, _ = self.create_features(demand_data)
        X_scaled = self.scaler.transform(X[-1:])  # Use last observation
        
        predictions = []
        uncertainties = []
        
        # Generate probabilistic forecasts
        for h in range(horizon):
            # Base prediction
            base_pred = self.model.predict(X_scaled)[0]
            
            # Add uncertainty based on model's training error patterns
            # In practice, you'd use more sophisticated uncertainty quantification
            residual_std = np.std(demand_data) * (1 + 0.1 * h)  # Increasing uncertainty over time
            pred_samples = np.random.normal(base_pred, residual_std, n_simulations)
            pred_samples = np.maximum(pred_samples, 0)  # Demand can't be negative
            
            predictions.append({
                'mean': base_pred,
                'median': np.median(pred_samples),
                'lower_80': np.percentile(pred_samples, 10),
                'upper_80': np.percentile(pred_samples, 90),
                'lower_95': np.percentile(pred_samples, 2.5),
                'upper_95': np.percentile(pred_samples, 97.5),
                'samples': pred_samples
            })
            
            # Update features for next prediction (simplified)
            # In practice, you'd properly roll the window
            X_scaled = np.roll(X_scaled, -1)
            X_scaled[0, -4:] = [h % 7, (h // 7) % 52, h / horizon, np.sin(2 * np.pi * h / 7)]
        
        return predictions

class InventoryOptimizer:
    """Operations Research based inventory optimizer"""
    
    def __init__(self, holding_cost=1.0, stockout_cost=10.0, order_cost=0.5):
        self.holding_cost = holding_cost
        self.stockout_cost = stockout_cost
        self.order_cost = order_cost
    
    def newsvendor_optimal_quantity(self, demand_samples):
        """Calculate optimal order quantity using newsvendor model"""
        # Critical ratio for newsvendor problem
        critical_ratio = self.stockout_cost / (self.stockout_cost + self.holding_cost)
        optimal_quantity = np.percentile(demand_samples, critical_ratio * 100)
        return max(0, optimal_quantity)
    
    def calculate_expected_cost(self, order_quantity, demand_samples):
        """Calculate expected total cost for given order quantity"""
        inventory_levels = order_quantity - demand_samples
        
        # Holding costs (positive inventory)
        holding_costs = np.maximum(inventory_levels, 0) * self.holding_cost
        
        # Stockout costs (negative inventory)
        stockout_costs = np.maximum(-inventory_levels, 0) * self.stockout_cost
        
        # Order costs
        order_costs = order_quantity * self.order_cost
        
        total_cost = np.mean(holding_costs + stockout_costs) + order_costs
        return total_cost
    
    def optimize_inventory(self, demand_predictions, initial_inventory=0):
        """Optimize inventory decisions over planning horizon"""
        horizon = len(demand_predictions)
        optimal_orders = []
        expected_costs = []
        inventory_levels = [initial_inventory]
        
        current_inventory = initial_inventory
        
        for t in range(horizon):
            demand_samples = demand_predictions[t]['samples']
            
            # Calculate net demand (demand - current inventory)
            net_demand_samples = np.maximum(demand_samples - current_inventory, 0)
            
            # Find optimal order quantity for this period
            if len(net_demand_samples[net_demand_samples > 0]) > 0:
                optimal_qty = self.newsvendor_optimal_quantity(net_demand_samples)
            else:
                optimal_qty = 0
                
            optimal_orders.append(optimal_qty)
            
            # Calculate expected cost
            total_inventory = current_inventory + optimal_qty
            cost = self.calculate_expected_cost(optimal_qty, demand_samples)
            expected_costs.append(cost)
            
            # Update inventory for next period (simplified simulation)
            expected_demand = demand_predictions[t]['mean']
            current_inventory = max(0, total_inventory - expected_demand)
            inventory_levels.append(current_inventory)
        
        return {
            'optimal_orders': optimal_orders,
            'expected_costs': expected_costs,
            'inventory_levels': inventory_levels[:-1]  # Remove last element
        }

# Generate synthetic demand data
def generate_demand_data(n_periods=200):
    """Generate realistic demand data with trend, seasonality, and randomness"""
    np.random.seed(42)
    
    # Base demand level
    base_demand = 100
    
    # Trend component
    trend = np.linspace(0, 20, n_periods)
    
    # Seasonal components
    weekly_season = 30 * np.sin(2 * np.pi * np.arange(n_periods) / 7)
    monthly_season = 15 * np.sin(2 * np.pi * np.arange(n_periods) / 30)
    
    # Random noise
    noise = np.random.normal(0, 10, n_periods)
    
    # Combine components
    demand = base_demand + trend + weekly_season + monthly_season + noise
    demand = np.maximum(demand, 0)  # Ensure non-negative demand
    
    return demand

# Main execution
print("🚀 Starting Integrated Demand Forecasting and Inventory Optimization")
print("=" * 70)

# Step 1: Generate and prepare data
demand_data = generate_demand_data(150)
train_size = int(0.8 * len(demand_data))
train_data = demand_data[:train_size]
test_data = demand_data[train_size:]

print(f"📊 Data Summary:")
print(f"   Total periods: {len(demand_data)}")
print(f"   Training periods: {len(train_data)}")
print(f"   Test periods: {len(test_data)}")
print(f"   Average demand: {np.mean(demand_data):.2f}")
print(f"   Demand std: {np.std(demand_data):.2f}")

# Step 2: Train demand forecasting model
print("\n🤖 Training Machine Learning Demand Forecaster...")
forecaster = DemandForecaster()
forecaster.fit(train_data)

# Step 3: Generate demand predictions
print("🔮 Generating Demand Predictions with Uncertainty...")
forecast_horizon = 30
demand_predictions = forecaster.predict_with_uncertainty(
    train_data, 
    horizon=forecast_horizon,
    n_simulations=1000
)

# Step 4: Optimize inventory
print("⚙️ Optimizing Inventory Decisions...")
optimizer = InventoryOptimizer(
    holding_cost=2.0,     # $2 per unit per period
    stockout_cost=15.0,   # $15 penalty per unit shortage
    order_cost=1.0        # $1 per unit ordered
)

optimization_results = optimizer.optimize_inventory(demand_predictions)

# Step 5: Performance Analysis
print("\n📈 Performance Analysis:")
print("-" * 30)

# Forecasting accuracy (using available test data)
test_predictions = []
for i in range(min(len(test_data), forecast_horizon)):
    test_predictions.append(demand_predictions[i]['mean'])

if len(test_predictions) > 0:
    mae = mean_absolute_error(test_data[:len(test_predictions)], test_predictions)
    rmse = np.sqrt(mean_squared_error(test_data[:len(test_predictions)], test_predictions))
    mape = np.mean(np.abs((test_data[:len(test_predictions)] - test_predictions) / test_data[:len(test_predictions)])) * 100
    
    print(f"Forecasting Performance:")
    print(f"  MAE: {mae:.2f}")
    print(f"  RMSE: {rmse:.2f}")
    print(f"  MAPE: {mape:.2f}%")

# Inventory optimization results
total_expected_cost = sum(optimization_results['expected_costs'])
avg_order_quantity = np.mean(optimization_results['optimal_orders'])
avg_inventory_level = np.mean(optimization_results['inventory_levels'])

print(f"\nInventory Optimization Results:")
print(f"  Total expected cost: ${total_expected_cost:.2f}")
print(f"  Average order quantity: {avg_order_quantity:.2f} units")
print(f"  Average inventory level: {avg_inventory_level:.2f} units")

print("\n✅ Analysis Complete! Generating visualizations...")

# Create comprehensive visualizations
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
fig.suptitle('Integrated Demand Forecasting and Inventory Optimization Results', fontsize=16, fontweight='bold')

# Plot 1: Historical demand and forecasts
ax1 = axes[0, 0]
time_historical = range(len(demand_data))
time_forecast = range(len(train_data), len(train_data) + forecast_horizon)

ax1.plot(time_historical, demand_data, 'b-', alpha=0.7, label='Historical Demand', linewidth=2)
ax1.axvline(x=len(train_data), color='red', linestyle='--', alpha=0.7, label='Train/Test Split')

# Plot forecast with uncertainty bands
forecast_means = [pred['mean'] for pred in demand_predictions]
forecast_lower_80 = [pred['lower_80'] for pred in demand_predictions]
forecast_upper_80 = [pred['upper_80'] for pred in demand_predictions]
forecast_lower_95 = [pred['lower_95'] for pred in demand_predictions]
forecast_upper_95 = [pred['upper_95'] for pred in demand_predictions]

ax1.plot(time_forecast, forecast_means, 'g-', linewidth=2, label='Forecast Mean')
ax1.fill_between(time_forecast, forecast_lower_80, forecast_upper_80, alpha=0.3, color='green', label='80% Confidence')
ax1.fill_between(time_forecast, forecast_lower_95, forecast_upper_95, alpha=0.1, color='green', label='95% Confidence')

ax1.set_title('Demand Forecasting with Uncertainty Quantification')
ax1.set_xlabel('Time Period')
ax1.set_ylabel('Demand')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2: Optimal order quantities and inventory levels
ax2 = axes[0, 1]
forecast_periods = range(forecast_horizon)

ax2_twin = ax2.twinx()
bars = ax2.bar(forecast_periods, optimization_results['optimal_orders'], alpha=0.6, color='orange', label='Optimal Orders')
line = ax2_twin.plot(forecast_periods, optimization_results['inventory_levels'], 'r-', linewidth=2, marker='o', label='Inventory Level')

ax2.set_title('Optimal Ordering Policy and Inventory Levels')
ax2.set_xlabel('Time Period')
ax2.set_ylabel('Order Quantity', color='orange')
ax2_twin.set_ylabel('Inventory Level', color='red')
ax2.grid(True, alpha=0.3)

# Create combined legend
lines1, labels1 = ax2.get_legend_handles_labels()
lines2, labels2 = ax2_twin.get_legend_handles_labels()
ax2.legend(lines1 + lines2, labels1 + labels2, loc='upper right')

# Plot 3: Cost breakdown analysis
ax3 = axes[1, 0]
cost_components = []
periods = range(forecast_horizon)

for i, pred in enumerate(demand_predictions):
    order_qty = optimization_results['optimal_orders'][i]
    demand_samples = pred['samples']
    
    # Calculate cost components
    inventory_after_order = optimization_results['inventory_levels'][i] + order_qty
    inventory_levels = inventory_after_order - demand_samples
    
    holding_cost = np.mean(np.maximum(inventory_levels, 0)) * optimizer.holding_cost
    stockout_cost = np.mean(np.maximum(-inventory_levels, 0)) * optimizer.stockout_cost
    order_cost = order_qty * optimizer.order_cost
    
    cost_components.append([holding_cost, stockout_cost, order_cost])

cost_components = np.array(cost_components)

ax3.stackplot(periods, 
              cost_components[:, 0], 
              cost_components[:, 1], 
              cost_components[:, 2],
              labels=['Holding Cost', 'Stockout Cost', 'Order Cost'],
              alpha=0.7)

ax3.set_title('Cost Component Analysis Over Time')
ax3.set_xlabel('Time Period')
ax3.set_ylabel('Expected Cost ($)')
ax3.legend(loc='upper left')
ax3.grid(True, alpha=0.3)

# Plot 4: Demand distribution and optimal policy
ax4 = axes[1, 1]
# Show distribution for first few periods
sample_period = 5
demand_samples = demand_predictions[sample_period]['samples']
optimal_order = optimization_results['optimal_orders'][sample_period]

ax4.hist(demand_samples, bins=30, alpha=0.7, density=True, color='skyblue', label='Demand Distribution')
ax4.axvline(x=np.mean(demand_samples), color='blue', linestyle='-', linewidth=2, label=f'Mean Demand: {np.mean(demand_samples):.1f}')
ax4.axvline(x=optimal_order, color='red', linestyle='--', linewidth=2, label=f'Optimal Order: {optimal_order:.1f}')

ax4.set_title(f'Demand Distribution and Optimal Policy (Period {sample_period+1})')
ax4.set_xlabel('Demand')
ax4.set_ylabel('Probability Density')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Additional analysis: Service level and fill rate
print("\n📊 Additional Performance Metrics:")
print("-" * 40)

service_levels = []
fill_rates = []

for i, pred in enumerate(demand_predictions):
    order_qty = optimization_results['optimal_orders'][i]
    current_inv = optimization_results['inventory_levels'][i]
    total_available = current_inv + order_qty
    demand_samples = pred['samples']
    
    # Service level (probability of no stockout)
    service_level = np.mean(demand_samples <= total_available)
    service_levels.append(service_level)
    
    # Fill rate (expected fraction of demand satisfied)
    satisfied_demand = np.minimum(demand_samples, total_available)
    fill_rate = np.mean(satisfied_demand / demand_samples)
    fill_rates.append(fill_rate)

avg_service_level = np.mean(service_levels)
avg_fill_rate = np.mean(fill_rates)

print(f"Average Service Level: {avg_service_level:.1%}")
print(f"Average Fill Rate: {avg_fill_rate:.1%}")
print(f"Total Inventory Investment: ${np.sum(optimization_results['optimal_orders']):.2f}")

# Summary insights
print("\n🎯 Key Insights:")
print("-" * 20)
print("• The ML model captures both trend and seasonal patterns in demand")
print("• Uncertainty quantification enables robust inventory decisions")
print("• The newsvendor model balances holding and stockout costs optimally")
print(f"• Service level of {avg_service_level:.1%} achieved with optimized policy")
print("• Integration of ML and OR provides superior performance than using either alone")

print("\n🏁 Analysis completed successfully!")

def create_features(self, data, lookback=7):

residual_std = np.std(demand_data) * (1 + 0.1 * h)
pred_samples = np.random.normal(base_pred, residual_std, n_simulations)

def newsvendor_optimal_quantity(self, demand_samples):
    critical_ratio = self.stockout_cost / (self.stockout_cost + self.holding_cost)
    optimal_quantity = np.percentile(demand_samples, critical_ratio * 100)

def calculate_expected_cost(self, order_quantity, demand_samples):
    inventory_levels = order_quantity - demand_samples
    holding_costs = np.maximum(inventory_levels, 0) * self.holding_cost
    stockout_costs = np.maximum(-inventory_levels, 0) * self.stockout_cost
    order_costs = order_quantity * self.order_cost

!pip install pulp

import pulp
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Employees and days
employees = ["Alice", "Bob", "Charlie", "Dana", "Eva"]
days = ["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"]

# Availability: 1 means available, 0 means unavailable
availability = {
    "Alice":   [1, 1, 1, 1, 1, 0, 0],
    "Bob":     [1, 1, 1, 1, 1, 1, 1],
    "Charlie": [1, 1, 1, 0, 0, 0, 0],
    "Dana":    [1, 1, 1, 1, 1, 1, 0],
    "Eva":     [0, 0, 1, 1, 1, 1, 1]
}

# Create model
model = pulp.LpProblem("Shift_Scheduling", pulp.LpMinimize)

# Decision variables
x = pulp.LpVariable.dicts("shift", (employees, days), cat='Binary')

# Objective: minimize total number of shifts
model += pulp.lpSum([x[e][d] for e in employees for d in days])

# Constraint 1: At least 2 staff per day
for d in days:
    model += pulp.lpSum([x[e][d] for e in employees]) >= 2

# Constraint 2: Max 5 shifts per employee
for e in employees:
    model += pulp.lpSum([x[e][d] for d in days]) <= 5

# Constraint 3: Availability
for e in employees:
    for i, d in enumerate(days):
        if availability[e][i] == 0:
            model += x[e][d] == 0

# Solve
model.solve()
print(f"Status: {pulp.LpStatus[model.status]}")

# Build schedule table
schedule = pd.DataFrame(0, index=employees, columns=days)
for e in employees:
    for d in days:
        schedule.loc[e, d] = int(pulp.value(x[e][d]))

# Plot
plt.figure(figsize=(10, 5))
sns.heatmap(schedule, annot=True, cbar=False, cmap="YlGnBu", linewidths=.5)
plt.title("Weekly Shift Schedule (1 = Working, 0 = Off)")
plt.xlabel("Day")
plt.ylabel("Employee")
plt.show()

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

# Parameters
mu1, mu2 = 0.1, 0.2  # Expected returns
sigma1, sigma2 = 0.15, 0.25  # Standard deviations
rho = 0.3  # Correlation
alpha = 0.05  # Risk tolerance (probability of loss)
w0 = 0.5  # Initial weight for improvement problem

# Portfolio expected return
def portfolio_return(w, mu1, mu2):
    return w * mu1 + (1 - w) * mu2

# Portfolio standard deviation
def portfolio_std(w, sigma1, sigma2, rho):
    return np.sqrt(w**2 * sigma1**2 + (1 - w)**2 * sigma2**2 + 2 * w * (1 - w) * rho * sigma1 * sigma2)

# Risk constraint: P(R_p < 0) <= alpha
def risk_constraint(w, mu1, mu2, sigma1, sigma2, rho, alpha):
    mu_p = portfolio_return(w, mu1, mu2)
    sigma_p = portfolio_std(w, sigma1, sigma2, rho)
    return mu_p + portfolio_std(w, sigma1, sigma2, rho) * norm.ppf(alpha)

# Objective function to maximize return (minimize negative return)
def objective(w):
    return -portfolio_return(w[0], mu1, mu2)

# Constraints and bounds
constraints = [
    {'type': 'ineq', 'fun': lambda w: risk_constraint(w[0], mu1, mu2, sigma1, sigma2, rho, alpha)}
]
bounds = [(0, 1)]  # w between 0 and 1

# Optimization
initial_guess = [w0]
result = minimize(objective, initial_guess, method='SLSQP', bounds=bounds, constraints=constraints)

# Extract results
optimal_w = result.x[0]
optimal_return = portfolio_return(optimal_w, mu1, mu2)
optimal_std = portfolio_std(optimal_w, sigma1, sigma2, rho)
initial_return = portfolio_return(w0, mu1, mu2)
initial_std = portfolio_std(w0, sigma1, sigma2, rho)

# Generate data for efficient frontier
w_values = np.linspace(0, 1, 100)
returns = [portfolio_return(w, mu1, mu2) for w in w_values]
stds = [portfolio_std(w, sigma1, sigma2, rho) for w in w_values]
feasible = [risk_constraint(w, mu1, mu2, sigma1, sigma2, rho, alpha) >= 0 for w in w_values]

# Plotting
plt.figure(figsize=(10, 6))
plt.scatter(stds, returns, c=feasible, cmap='viridis', label='Portfolios')
plt.scatter(initial_std, initial_return, c='red', s=100, label='Initial Portfolio (w=0.5)')
plt.scatter(optimal_std, optimal_return, c='blue', s=100, label=f'Optimal Portfolio (w={optimal_w:.2f})')
plt.plot(stds, returns, 'k--', label='Efficient Frontier')
plt.axhline(0, color='gray', linestyle='--')
plt.title('Portfolio Optimization with Risk Constraint')
plt.xlabel('Standard Deviation (Risk)')
plt.ylabel('Expected Return')
plt.legend()
plt.grid(True)
plt.colorbar(label='Feasible (Green) / Infeasible (Yellow)')
plt.savefig('portfolio_optimization.png')
plt.show()

# Print results
print(f"Optimal weight for asset 1: {optimal_w:.4f}")
print(f"Optimal expected return: {optimal_return:.4f}")
print(f"Optimal standard deviation: {optimal_std:.4f}")
print(f"Initial expected return: {initial_return:.4f}")
print(f"Initial standard deviation: {initial_std:.4f}")