import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize, differential_evolution
from scipy.special import comb
import seaborn as sns
from mpl_toolkits.mplot3d import Axes3D
import warnings
warnings.filterwarnings('ignore')

class QuantumSecretSharingOptimizer:
    def __init__(self, alpha=1.0, beta=2.0, gamma=0.5, p_compromise=0.1, max_resource=100):
        """
        Initialize the quantum secret sharing optimizer
        
        Parameters:
        alpha: Weight for quantum overhead cost
        beta: Weight for security benefit  
        gamma: Weight for information entropy
        p_compromise: Probability of party compromise
        max_resource: Maximum resource budget
        """
        self.alpha = alpha
        self.beta = beta  
        self.gamma = gamma
        self.p_compromise = p_compromise
        self.max_resource = max_resource
        
    def information_entropy(self, k, n):
        """
        Calculate information entropy H(k,n) for the quantum secret sharing scheme
        
        H(k,n) = -∑_{i=k}^{n} C(n,i) * p^i * (1-p)^{n-i} * log2(C(n,i) * p^i * (1-p)^{n-i})
        """
        entropy = 0.0
        p = self.p_compromise
        
        for i in range(k, n+1):
            if i <= n:  # Ensure valid binomial coefficient
                prob_i = comb(n, i, exact=False) * (p**i) * ((1-p)**(n-i))
                if prob_i > 1e-15:  # Avoid log(0)
                    entropy -= prob_i * np.log2(prob_i)
        
        return entropy
    
    def cost_function(self, params):
        """
        Cost function to minimize: C(k,n) = α(n/k) + β(1-(k-1)/(n-1)) + γH(k,n)
        """
        k, n = int(params[0]), int(params[1])
        
        # Constraint checks
        if k < 2 or n < k or n > 20 or k * n > self.max_resource:
            return 1e6  # Large penalty for invalid parameters
            
        # Quantum overhead cost
        overhead_cost = self.alpha * (n / k)
        
        # Security benefit (higher k relative to n is more secure)
        if n > 1:
            security_benefit = self.beta * (1 - (k-1)/(n-1))
        else:
            security_benefit = self.beta
            
        # Information entropy cost
        entropy_cost = self.gamma * self.information_entropy(k, n)
        
        total_cost = overhead_cost - security_benefit + entropy_cost
        
        return total_cost
    
    def optimize_parameters(self):
        """
        Find optimal k and n parameters using differential evolution
        """
        # Define bounds: k in [2, 10], n in [k, 20]
        bounds = [(2, 10), (2, 20)]
        
        # Use differential evolution for global optimization
        result = differential_evolution(
            self.cost_function,
            bounds,
            seed=42,
            maxiter=1000,
            popsize=15,
            atol=1e-6,
            tol=1e-6
        )
        
        optimal_k = int(result.x[0])
        optimal_n = int(result.x[1])
        optimal_cost = result.fun
        
        return optimal_k, optimal_n, optimal_cost, result
    
    def analyze_parameter_space(self):
        """
        Analyze the cost function across different k and n values
        """
        k_range = range(2, 11)
        n_range = range(2, 21)
        
        cost_matrix = np.zeros((len(k_range), len(n_range)))
        
        for i, k in enumerate(k_range):
            for j, n in enumerate(n_range):
                if n >= k:  # Valid constraint
                    cost_matrix[i, j] = self.cost_function([k, n])
                else:
                    cost_matrix[i, j] = np.nan
                    
        return k_range, n_range, cost_matrix
    
    def security_analysis(self, k, n):
        """
        Analyze security properties of the (k,n) scheme
        """
        # Probability that fewer than k parties are compromised (secure)
        p = self.p_compromise
        secure_prob = sum(comb(n, i, exact=False) * (p**i) * ((1-p)**(n-i)) 
                         for i in range(k))
        
        # Expected number of compromised parties
        expected_compromised = n * p
        
        # Security margin
        security_margin = k - expected_compromised
        
        return {
            'secure_probability': secure_prob,
            'expected_compromised': expected_compromised,
            'security_margin': security_margin,
            'quantum_overhead': n / k,
            'information_entropy': self.information_entropy(k, n)
        }

def plot_optimization_results(optimizer):
    """
    Create comprehensive visualizations of the optimization results
    """
    # Get optimal parameters
    opt_k, opt_n, opt_cost, result = optimizer.optimize_parameters()
    
    # Analyze parameter space
    k_range, n_range, cost_matrix = optimizer.analyze_parameter_space()
    
    # Create figure with subplots
    fig = plt.figure(figsize=(20, 15))
    
    # 1. Cost function heatmap
    ax1 = plt.subplot(2, 3, 1)
    mask = np.isnan(cost_matrix)
    sns.heatmap(cost_matrix, 
                xticklabels=n_range, 
                yticklabels=k_range,
                annot=True, 
                fmt='.2f', 
                cmap='viridis_r',
                mask=mask,
                cbar_kws={'label': 'Cost Function Value'})
    plt.title('Cost Function Heatmap\n$C(k,n) = \\alpha\\frac{n}{k} + \\beta(1-\\frac{k-1}{n-1}) + \\gamma H(k,n)$')
    plt.xlabel('Number of Parties (n)')
    plt.ylabel('Threshold (k)')
    
    # Mark optimal point
    opt_k_idx = list(k_range).index(opt_k)
    opt_n_idx = list(n_range).index(opt_n)
    plt.scatter(opt_n_idx + 0.5, opt_k_idx + 0.5, 
               color='red', s=200, marker='*', 
               label=f'Optimal: k={opt_k}, n={opt_n}')
    plt.legend()
    
    # 2. 3D surface plot
    ax2 = plt.subplot(2, 3, 2, projection='3d')
    K, N = np.meshgrid(k_range, n_range)
    # Transpose cost_matrix to match meshgrid orientation
    cost_surface = cost_matrix.T
    
    # Create mask for valid combinations
    valid_mask = N >= K
    K_valid = K[valid_mask]
    N_valid = N[valid_mask]
    cost_valid = cost_surface[valid_mask]
    
    scatter = ax2.scatter(K_valid, N_valid, cost_valid, 
                         c=cost_valid, cmap='viridis_r', s=50)
    ax2.scatter([opt_k], [opt_n], [opt_cost], 
               color='red', s=200, marker='*')
    
    ax2.set_xlabel('Threshold (k)')
    ax2.set_ylabel('Number of Parties (n)')
    ax2.set_zlabel('Cost Function Value')
    ax2.set_title('3D Cost Function Surface')
    plt.colorbar(scatter, ax=ax2, shrink=0.8)
    
    # 3. Security analysis
    ax3 = plt.subplot(2, 3, 3)
    security_data = []
    k_values = []
    n_values = []
    
    for k in range(2, 11):
        for n in range(k, min(k+10, 21)):
            security = optimizer.security_analysis(k, n)
            security_data.append(security['secure_probability'])
            k_values.append(k)
            n_values.append(n)
    
    scatter = ax3.scatter(k_values, n_values, c=security_data, 
                         cmap='RdYlGn', s=60, alpha=0.7)
    ax3.scatter([opt_k], [opt_n], color='red', s=200, marker='*', 
               label=f'Optimal: k={opt_k}, n={opt_n}')
    
    plt.colorbar(scatter, ax=ax3, label='Security Probability')
    ax3.set_xlabel('Threshold (k)')
    ax3.set_ylabel('Number of Parties (n)')
    ax3.set_title('Security Probability Analysis')
    ax3.legend()
    
    # 4. Quantum overhead vs Security trade-off
    ax4 = plt.subplot(2, 3, 4)
    overhead_data = []
    security_data = []
    
    for k in range(2, 11):
        for n in range(k, min(k+8, 16)):
            security = optimizer.security_analysis(k, n)
            overhead_data.append(security['quantum_overhead'])
            security_data.append(security['secure_probability'])
    
    scatter = ax4.scatter(overhead_data, security_data, 
                         c=range(len(overhead_data)), cmap='plasma', alpha=0.7)
    
    # Mark optimal point
    opt_security = optimizer.security_analysis(opt_k, opt_n)
    ax4.scatter([opt_security['quantum_overhead']], 
               [opt_security['secure_probability']], 
               color='red', s=200, marker='*', 
               label=f'Optimal: ({opt_security["quantum_overhead"]:.2f}, {opt_security["secure_probability"]:.3f})')
    
    ax4.set_xlabel('Quantum Overhead (n/k)')
    ax4.set_ylabel('Security Probability')
    ax4.set_title('Quantum Overhead vs Security Trade-off')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 5. Information entropy analysis
    ax5 = plt.subplot(2, 3, 5)
    entropy_matrix = np.zeros((len(k_range), len(n_range)))
    
    for i, k in enumerate(k_range):
        for j, n in enumerate(n_range):
            if n >= k:
                entropy_matrix[i, j] = optimizer.information_entropy(k, n)
            else:
                entropy_matrix[i, j] = np.nan
    
    mask = np.isnan(entropy_matrix)
    sns.heatmap(entropy_matrix, 
                xticklabels=n_range, 
                yticklabels=k_range,
                annot=True, 
                fmt='.3f', 
                cmap='Blues',
                mask=mask,
                cbar_kws={'label': 'Information Entropy'})
    
    plt.scatter(opt_n_idx + 0.5, opt_k_idx + 0.5, 
               color='red', s=200, marker='*')
    plt.title('Information Entropy $H(k,n)$')
    plt.xlabel('Number of Parties (n)')
    plt.ylabel('Threshold (k)')
    
    # 6. Cost components breakdown
    ax6 = plt.subplot(2, 3, 6)
    
    # Calculate cost components for optimal solution
    overhead_cost = optimizer.alpha * (opt_n / opt_k)
    security_benefit = optimizer.beta * (1 - (opt_k-1)/(opt_n-1)) if opt_n > 1 else optimizer.beta
    entropy_cost = optimizer.gamma * optimizer.information_entropy(opt_k, opt_n)
    
    components = ['Quantum\nOverhead', 'Security\nBenefit', 'Information\nEntropy']
    values = [overhead_cost, -security_benefit, entropy_cost]  # Negative for benefit
    colors = ['red', 'green', 'blue']
    
    bars = ax6.bar(components, values, color=colors, alpha=0.7)
    ax6.axhline(y=0, color='black', linestyle='-', linewidth=0.8)
    ax6.set_ylabel('Cost Component Value')
    ax6.set_title(f'Cost Components Breakdown\nOptimal: k={opt_k}, n={opt_n}')
    
    # Add value labels on bars
    for bar, value in zip(bars, values):
        height = bar.get_height()
        ax6.text(bar.get_x() + bar.get_width()/2., height + (0.1 if height > 0 else -0.1),
                f'{value:.3f}', ha='center', va='bottom' if height > 0 else 'top')
    
    plt.tight_layout()
    plt.show()
    
    return opt_k, opt_n, opt_cost

# Run the optimization
print("=== Quantum Secret Sharing Optimization ===\n")

# Initialize optimizer with specific parameters
optimizer = QuantumSecretSharingOptimizer(
    alpha=1.0,    # Weight for quantum overhead
    beta=2.0,     # Weight for security benefit
    gamma=0.5,    # Weight for information entropy
    p_compromise=0.1,  # 10% probability of party compromise
    max_resource=100   # Resource budget constraint
)

print("Optimizer Parameters:")
print(f"α (overhead weight): {optimizer.alpha}")
print(f"β (security weight): {optimizer.beta}")  
print(f"γ (entropy weight): {optimizer.gamma}")
print(f"Compromise probability: {optimizer.p_compromise}")
print(f"Resource budget: {optimizer.max_resource}\n")

# Perform optimization
opt_k, opt_n, opt_cost = plot_optimization_results(optimizer)

print("=== Optimization Results ===")
print(f"Optimal threshold (k): {opt_k}")
print(f"Optimal number of parties (n): {opt_n}")
print(f"Optimal cost: {opt_cost:.4f}\n")

# Detailed analysis of optimal solution
optimal_security = optimizer.security_analysis(opt_k, opt_n)

print("=== Security Analysis of Optimal Solution ===")
print(f"Security probability: {optimal_security['secure_probability']:.4f}")
print(f"Expected compromised parties: {optimal_security['expected_compromised']:.2f}")
print(f"Security margin: {optimal_security['security_margin']:.2f}")
print(f"Quantum overhead (n/k): {optimal_security['quantum_overhead']:.2f}")
print(f"Information entropy: {optimal_security['information_entropy']:.4f}")

print("\n=== Cost Function Components ===")
overhead = optimizer.alpha * (opt_n / opt_k)
security = optimizer.beta * (1 - (opt_k-1)/(opt_n-1))
entropy = optimizer.gamma * optimal_security['information_entropy']
total = overhead - security + entropy

print(f"Quantum overhead cost: α(n/k) = {overhead:.4f}")
print(f"Security benefit: β(1-(k-1)/(n-1)) = {security:.4f}")
print(f"Information entropy cost: γH(k,n) = {entropy:.4f}")
print(f"Total cost: {total:.4f}")

import numpy as np
import matplotlib.pyplot as plt
from itertools import product, combinations
import seaborn as sns
from scipy.optimize import minimize
from mpl_toolkits.mplot3d import Axes3D

# Set up the plotting style
plt.style.use('seaborn-v0_8')
sns.set_palette("husl")

class SteaneCode:
    """
    Implementation of the 7-qubit Steane quantum error correction code.
    This code can correct any single qubit error (X, Y, or Z errors).
    """
    
    def __init__(self):
        # Generator matrix for the Steane code (stabilizer generators)
        # These define the parity check constraints
        self.stabilizers_x = np.array([
            [1, 1, 1, 1, 0, 0, 0],  # X₁X₂X₃X₄
            [1, 1, 0, 0, 1, 1, 0],  # X₁X₂X₅X₆  
            [1, 0, 1, 0, 1, 0, 1]   # X₁X₃X₅X₇
        ])
        
        self.stabilizers_z = np.array([
            [1, 1, 1, 1, 0, 0, 0],  # Z₁Z₂Z₃Z₄
            [1, 1, 0, 0, 1, 1, 0],  # Z₁Z₂Z₅Z₆
            [1, 0, 1, 0, 1, 0, 1]   # Z₁Z₃Z₅Z₇
        ])
        
        # Logical operators
        self.logical_x = np.array([1, 1, 1, 1, 1, 1, 1])  # X̄ = X₁X₂X₃X₄X₅X₆X₇
        self.logical_z = np.array([1, 1, 1, 1, 1, 1, 1])  # Z̄ = Z₁Z₂Z₃Z₄Z₅Z₆Z₇
    
    def syndrome(self, error_pattern, error_type='X'):
        """
        Calculate the syndrome for a given error pattern.
        The syndrome uniquely identifies single qubit errors.
        """
        if error_type == 'X':
            stabilizers = self.stabilizers_z  # Z stabilizers detect X errors
        else:  # error_type == 'Z'
            stabilizers = self.stabilizers_x  # X stabilizers detect Z errors
            
        return (stabilizers @ error_pattern) % 2
    
    def decode_syndrome(self, syndrome, error_type='X'):
        """
        Decode the syndrome to find the error location.
        Returns the qubit index where the error occurred (0-6).
        """
        # Create syndrome lookup table for single qubit errors
        syndrome_table = {}
        for i in range(7):
            error = np.zeros(7, dtype=int)
            error[i] = 1
            syn = self.syndrome(error, error_type)
            syndrome_table[tuple(syn)] = i
        
        return syndrome_table.get(tuple(syndrome), -1)  # -1 if no single error found

class QuantumErrorSimulator:
    """
    Simulator for quantum errors and error correction performance analysis.
    """
    
    def __init__(self, code):
        self.code = code
        
    def generate_random_errors(self, num_errors, error_types=['X', 'Z', 'Y']):
        """
        Generate random quantum errors for simulation.
        Y errors are treated as both X and Z errors occurring simultaneously.
        """
        errors = []
        for _ in range(num_errors):
            error_type = np.random.choice(error_types)
            qubit_pos = np.random.randint(0, 7)
            
            if error_type == 'Y':
                # Y error = X error + Z error
                x_error = np.zeros(7, dtype=int)
                z_error = np.zeros(7, dtype=int)
                x_error[qubit_pos] = 1
                z_error[qubit_pos] = 1
                errors.append(('Y', x_error, z_error))
            else:
                error = np.zeros(7, dtype=int)
                error[qubit_pos] = 1
                if error_type == 'X':
                    errors.append(('X', error, np.zeros(7, dtype=int)))
                else:  # error_type == 'Z'
                    errors.append(('Z', np.zeros(7, dtype=int), error))
        
        return errors
    
    def test_error_correction(self, errors):
        """
        Test the error correction capability for a list of errors.
        Returns success rate and detailed results.
        """
        results = []
        successful_corrections = 0
        
        for error_type, x_error, z_error in errors:
            # Calculate syndromes
            x_syndrome = self.code.syndrome(x_error, 'X') if np.any(x_error) else np.zeros(3, dtype=int)
            z_syndrome = self.code.syndrome(z_error, 'Z') if np.any(z_error) else np.zeros(3, dtype=int)
            
            # Decode syndromes
            x_correction = self.code.decode_syndrome(x_syndrome, 'X') if np.any(x_syndrome) else -1
            z_correction = self.code.decode_syndrome(z_syndrome, 'Z') if np.any(z_syndrome) else -1
            
            # Check if correction is successful
            x_success = (x_correction == np.argmax(x_error)) if np.any(x_error) else True
            z_success = (z_correction == np.argmax(z_error)) if np.any(z_error) else True
            
            overall_success = x_success and z_success
            if overall_success:
                successful_corrections += 1
            
            results.append({
                'error_type': error_type,
                'x_error_pos': np.argmax(x_error) if np.any(x_error) else -1,
                'z_error_pos': np.argmax(z_error) if np.any(z_error) else -1,
                'x_syndrome': x_syndrome,
                'z_syndrome': z_syndrome,
                'x_correction': x_correction,
                'z_correction': z_correction,
                'success': overall_success
            })
        
        success_rate = successful_corrections / len(errors)
        return success_rate, results

def optimize_threshold_analysis(steane_code, error_rates):
    """
    Analyze the error threshold - the maximum error rate below which 
    error correction provides a net benefit.
    """
    simulator = QuantumErrorSimulator(steane_code)
    success_rates = []
    
    for error_rate in error_rates:
        # Simulate errors with given probability
        num_trials = 1000
        total_success = 0
        
        for _ in range(num_trials):
            # Each qubit has error_rate probability of having an error
            errors = []
            for qubit in range(7):
                if np.random.random() < error_rate:
                    error_type = np.random.choice(['X', 'Z', 'Y'])
                    if error_type == 'Y':
                        x_error = np.zeros(7, dtype=int)
                        z_error = np.zeros(7, dtype=int)
                        x_error[qubit] = 1
                        z_error[qubit] = 1
                        errors.append(('Y', x_error, z_error))
                    else:
                        error = np.zeros(7, dtype=int)
                        error[qubit] = 1
                        if error_type == 'X':
                            errors.append(('X', error, np.zeros(7, dtype=int)))
                        else:
                            errors.append(('Z', np.zeros(7, dtype=int), error))
            
            if errors:  # Only test if there are errors
                success_rate, _ = simulator.test_error_correction(errors)
                total_success += success_rate
            else:
                total_success += 1  # No errors = perfect success
                
        success_rates.append(total_success / num_trials)
    
    return success_rates

def distance_analysis():
    """
    Analyze the distance properties of the Steane code.
    The distance is the minimum weight of non-trivial logical operators.
    """
    steane = SteaneCode()
    
    # Calculate minimum distance
    # For the Steane code, the distance should be 3
    logical_x_weight = np.sum(steane.logical_x)
    logical_z_weight = np.sum(steane.logical_z)
    
    return min(logical_x_weight, logical_z_weight)

# Main analysis and optimization
def main_analysis():
    """
    Main function to run comprehensive analysis of the Steane code optimization.
    """
    print("=== Quantum Error Correction Code Optimization Analysis ===\n")
    
    # Initialize the Steane code
    steane_code = SteaneCode()
    simulator = QuantumErrorSimulator(steane_code)
    
    # 1. Basic error correction testing
    print("1. Testing single qubit error correction:")
    single_errors = simulator.generate_random_errors(100, ['X', 'Z', 'Y'])
    success_rate, results = simulator.test_error_correction(single_errors)
    print(f"   Success rate for single qubit errors: {success_rate:.2%}")
    
    # 2. Distance analysis
    distance = distance_analysis()
    print(f"2. Code distance: {distance}")
    print(f"   This means the code can correct up to {(distance-1)//2} errors")
    
    # 3. Threshold analysis
    print("\n3. Performing threshold analysis...")
    error_rates = np.logspace(-4, -1, 20)  # From 0.01% to 10%
    success_rates = optimize_threshold_analysis(steane_code, error_rates)
    
    return {
        'steane_code': steane_code,
        'single_error_success_rate': success_rate,
        'single_error_results': results,
        'distance': distance,
        'error_rates': error_rates,
        'success_rates': success_rates
    }

# Run the analysis
results = main_analysis()

# Create comprehensive visualizations
def create_comprehensive_plots(results):
    """
    Create detailed plots for the quantum error correction analysis.
    """
    fig = plt.figure(figsize=(20, 15))
    
    # Plot 1: Error correction success rate vs error rate (Threshold analysis)
    ax1 = plt.subplot(2, 3, 1)
    plt.semilogx(results['error_rates'], results['success_rates'], 'b-', linewidth=3, marker='o')
    plt.axhline(y=0.5, color='r', linestyle='--', alpha=0.7, label='50% threshold')
    plt.xlabel('Physical Error Rate', fontsize=12)
    plt.ylabel('Logical Success Rate', fontsize=12)
    plt.title('Error Threshold Analysis\nSteane Code Performance', fontsize=14, fontweight='bold')
    plt.grid(True, alpha=0.3)
    plt.legend()
    
    # Find approximate threshold
    threshold_idx = np.where(np.array(results['success_rates']) < 0.5)[0]
    if len(threshold_idx) > 0:
        threshold = results['error_rates'][threshold_idx[0]]
        plt.axvline(x=threshold, color='g', linestyle=':', alpha=0.7, 
                   label=f'Threshold ≈ {threshold:.1e}')
        plt.legend()
    
    # Plot 2: Syndrome pattern analysis
    ax2 = plt.subplot(2, 3, 2)
    steane = results['steane_code']
    
    # Create syndrome lookup table for visualization
    syndrome_patterns = []
    error_positions = []
    for i in range(7):
        x_error = np.zeros(7, dtype=int)
        x_error[i] = 1
        x_syndrome = steane.syndrome(x_error, 'X')
        syndrome_patterns.append(x_syndrome)
        error_positions.append(i)
    
    syndrome_matrix = np.array(syndrome_patterns).T
    im = plt.imshow(syndrome_matrix, cmap='RdYlBu', aspect='auto')
    plt.colorbar(im, ax=ax2)
    plt.xlabel('Qubit Position', fontsize=12)
    plt.ylabel('Syndrome Bit', fontsize=12)
    plt.title('X-Error Syndrome Patterns\nSteane Code', fontsize=14, fontweight='bold')
    plt.xticks(range(7), [f'Q{i}' for i in range(7)])
    plt.yticks(range(3), ['S₀', 'S₁', 'S₂'])
    
    # Plot 3: Error type distribution analysis
    ax3 = plt.subplot(2, 3, 3)
    error_types = [r['error_type'] for r in results['single_error_results']]
    success_by_type = {}
    for error_type in ['X', 'Z', 'Y']:
        type_results = [r for r in results['single_error_results'] if r['error_type'] == error_type]
        if type_results:
            success_by_type[error_type] = np.mean([r['success'] for r in type_results])
        else:
            success_by_type[error_type] = 0
    
    bars = plt.bar(success_by_type.keys(), success_by_type.values(), 
                   color=['red', 'blue', 'green'], alpha=0.7)
    plt.ylabel('Correction Success Rate', fontsize=12)
    plt.xlabel('Error Type', fontsize=12)
    plt.title('Success Rate by Error Type\n(Single Qubit Errors)', fontsize=14, fontweight='bold')
    plt.ylim(0, 1.1)
    
    # Add value labels on bars
    for bar, value in zip(bars, success_by_type.values()):
        plt.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.02, 
                f'{value:.2%}', ha='center', va='bottom', fontsize=11)
    
    # Plot 4: Code parameters visualization
    ax4 = plt.subplot(2, 3, 4)
    params = {
        'Physical Qubits (n)': 7,
        'Logical Qubits (k)': 1, 
        'Distance (d)': results['distance'],
        'Correctable Errors': (results['distance']-1)//2
    }
    
    y_pos = np.arange(len(params))
    values = list(params.values())
    bars = plt.barh(y_pos, values, color='skyblue', alpha=0.8)
    plt.yticks(y_pos, list(params.keys()))
    plt.xlabel('Value', fontsize=12)
    plt.title('Steane Code Parameters\n[[7,1,3]] Code', fontsize=14, fontweight='bold')
    
    # Add value labels
    for i, (bar, value) in enumerate(zip(bars, values)):
        plt.text(bar.get_width() + 0.1, bar.get_y() + bar.get_height()/2, 
                str(value), ha='left', va='center', fontsize=11)
    
    # Plot 5: Stabilizer generator matrix visualization
    ax5 = plt.subplot(2, 3, 5)
    stabilizers = np.vstack([steane.stabilizers_x, steane.stabilizers_z])
    im = plt.imshow(stabilizers, cmap='RdBu', aspect='auto')
    plt.colorbar(im, ax=ax5)
    plt.xlabel('Qubit Position', fontsize=12)
    plt.ylabel('Stabilizer Generator', fontsize=12)
    plt.title('Stabilizer Generators\nX-type (top) and Z-type (bottom)', fontsize=14, fontweight='bold')
    plt.xticks(range(7), [f'Q{i}' for i in range(7)])
    plt.yticks(range(6), ['X₁', 'X₂', 'X₃', 'Z₁', 'Z₂', 'Z₃'])
    
    # Plot 6: Error correction efficiency analysis
    ax6 = plt.subplot(2, 3, 6)
    
    # Calculate encoding efficiency metrics
    encoding_rate = 1/7  # k/n = 1/7
    overhead = 7/1  # n/k = 7/1
    
    # Create a comparison with other codes (theoretical)
    codes = ['Steane [7,1,3]', 'Repetition [3,1,3]', 'Shor [9,1,3]', 'Surface [17,1,5]']
    rates = [1/7, 1/3, 1/9, 1/17]
    distances = [3, 3, 3, 5]
    
    scatter = plt.scatter(rates, distances, s=200, alpha=0.7, 
                         c=['red', 'blue', 'green', 'orange'])
    plt.xlabel('Encoding Rate (k/n)', fontsize=12)
    plt.ylabel('Code Distance', fontsize=12)
    plt.title('Code Efficiency Comparison\nRate vs Distance Trade-off', fontsize=14, fontweight='bold')
    
    # Add labels for each point
    for i, code in enumerate(codes):
        plt.annotate(code, (rates[i], distances[i]), 
                    xytext=(5, 5), textcoords='offset points', fontsize=10)
    
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()
    
    # Additional 3D visualization of error correction landscape
    fig2 = plt.figure(figsize=(12, 8))
    ax = fig2.add_subplot(111, projection='3d')
    
    # Create a 3D surface showing error correction capability
    x_errors = np.arange(8)  # 0-7 possible X errors
    z_errors = np.arange(8)  # 0-7 possible Z errors
    X, Z = np.meshgrid(x_errors, z_errors)
    
    # Calculate correction success for combinations of errors
    success_surface = np.zeros_like(X, dtype=float)
    for i, x_err_count in enumerate(x_errors):
        for j, z_err_count in enumerate(z_errors):
            total_errors = x_err_count + z_err_count
            if total_errors == 0:
                success_surface[j, i] = 1.0
            elif total_errors == 1:
                success_surface[j, i] = 1.0  # Single errors always correctable
            elif total_errors == 2:
                success_surface[j, i] = 0.5  # Some two-error patterns correctable
            else:
                success_surface[j, i] = 0.1  # Multiple errors rarely correctable
    
    surf = ax.plot_surface(X, Z, success_surface, cmap='viridis', alpha=0.8)
    ax.set_xlabel('Number of X Errors')
    ax.set_ylabel('Number of Z Errors') 
    ax.set_zlabel('Correction Success Probability')
    ax.set_title('3D Error Correction Landscape\nSteane Code Performance')
    
    fig2.colorbar(surf)
    plt.show()

# Generate all plots
create_comprehensive_plots(results)

# Print detailed analysis results
print("\n=== Detailed Analysis Results ===")
print(f"Code Parameters: [[7,1,{results['distance']}]]")
print(f"Encoding Rate: {1/7:.3f}")
print(f"Single Error Correction Success: {results['single_error_success_rate']:.2%}")

# Syndrome table analysis
print("\nSyndrome Analysis:")
steane = results['steane_code']
print("X-Error Syndrome Lookup Table:")
for i in range(7):
    x_error = np.zeros(7, dtype=int)
    x_error[i] = 1
    syndrome = steane.syndrome(x_error, 'X')
    print(f"  Qubit {i}: Syndrome {syndrome} → Binary: {np.binary_repr(syndrome[0], 1)}{np.binary_repr(syndrome[1], 1)}{np.binary_repr(syndrome[2], 1)}")

print(f"\nOptimization Insights:")
print(f"• The Steane code achieves optimal parameters for single error correction")
print(f"• Distance-3 provides the minimum overhead for correcting 1 error")
print(f"• Symmetric design enables correction of all Pauli errors (X, Y, Z)")
print(f"• Threshold analysis shows practical error rates where QEC provides benefit")

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar, minimize
from scipy.linalg import logm
import seaborn as sns

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

class QuantumDepolarizingChannel:
    """
    Implementation of a quantum depolarizing channel for capacity optimization.
    """
    
    def __init__(self, p):
        """
        Initialize the depolarizing channel.
        
        Parameters:
        p (float): Depolarization parameter between 0 and 1
        """
        self.p = p
        
    def apply_channel(self, rho):
        """
        Apply the depolarizing channel to a density matrix.
        
        The depolarizing channel is defined as:
        Φ(ρ) = (1-p)ρ + p(I/2)
        
        Parameters:
        rho (numpy.ndarray): Input density matrix
        
        Returns:
        numpy.ndarray: Output density matrix after channel application
        """
        identity = np.eye(2) / 2  # Maximally mixed state
        return (1 - self.p) * rho + self.p * identity
    
    def von_neumann_entropy(self, rho):
        """
        Calculate the von Neumann entropy of a density matrix.
        
        S(ρ) = -Tr(ρ log₂(ρ))
        
        Parameters:
        rho (numpy.ndarray): Density matrix
        
        Returns:
        float: von Neumann entropy
        """
        # Add small regularization to avoid log(0)
        eigenvals = np.real(np.linalg.eigvals(rho))
        eigenvals = eigenvals[eigenvals > 1e-12]
        
        if len(eigenvals) == 0:
            return 0.0
            
        return -np.sum(eigenvals * np.log2(eigenvals))
    
    def mutual_information(self, input_probs, input_states):
        """
        Calculate mutual information I(X:Y) = H(Y) - H(Y|X)
        
        Parameters:
        input_probs (numpy.ndarray): Probabilities of input states
        input_states (list): List of input density matrices
        
        Returns:
        float: Mutual information
        """
        # Calculate average output state
        avg_output = np.zeros((2, 2), dtype=complex)
        for i, state in enumerate(input_states):
            output_state = self.apply_channel(state)
            avg_output += input_probs[i] * output_state
        
        # H(Y) - entropy of average output
        h_y = self.von_neumann_entropy(avg_output)
        
        # H(Y|X) - conditional entropy
        h_y_given_x = 0.0
        for i, state in enumerate(input_states):
            output_state = self.apply_channel(state)
            h_y_given_x += input_probs[i] * self.von_neumann_entropy(output_state)
        
        return h_y - h_y_given_x
    
    def theoretical_capacity(self):
        """
        Calculate the theoretical capacity of the depolarizing channel.
        
        Returns:
        float: Theoretical channel capacity
        """
        if self.p == 0:
            return 1.0
        elif self.p == 1:
            return 0.0
        else:
            term1 = 1 + (1 - self.p) * np.log2(1 - self.p)
            term2 = (self.p / 3) * np.log2(self.p / 3) if self.p > 0 else 0
            return term1 + term2

def create_qubit_state(theta, phi):
    """
    Create a general qubit state |ψ⟩ = cos(θ/2)|0⟩ + e^(iφ)sin(θ/2)|1⟩
    
    Parameters:
    theta (float): Polar angle
    phi (float): Azimuthal angle
    
    Returns:
    numpy.ndarray: Density matrix of the qubit state
    """
    psi = np.array([np.cos(theta/2), np.exp(1j * phi) * np.sin(theta/2)])
    return np.outer(psi, np.conj(psi))

def optimize_capacity_single_input(p_values):
    """
    Optimize capacity for single input state across different depolarization parameters.
    
    Parameters:
    p_values (numpy.ndarray): Array of depolarization parameters
    
    Returns:
    tuple: Arrays of optimized capacities and optimal parameters
    """
    optimized_capacities = []
    optimal_thetas = []
    optimal_phis = []
    
    for p in p_values:
        channel = QuantumDepolarizingChannel(p)
        
        def objective(params):
            theta, phi = params
            state = create_qubit_state(theta, phi)
            return -channel.mutual_information([1.0], [state])
        
        # Optimize over theta and phi
        result = minimize(objective, [np.pi/2, 0], 
                         bounds=[(0, np.pi), (0, 2*np.pi)],
                         method='L-BFGS-B')
        
        optimized_capacities.append(-result.fun)
        optimal_thetas.append(result.x[0])
        optimal_phis.append(result.x[1])
    
    return np.array(optimized_capacities), np.array(optimal_thetas), np.array(optimal_phis)

def optimize_capacity_two_inputs(p_values):
    """
    Optimize capacity using two input states with optimal probabilities.
    
    Parameters:
    p_values (numpy.ndarray): Array of depolarization parameters
    
    Returns:
    tuple: Arrays of optimized capacities and optimal parameters
    """
    optimized_capacities = []
    optimal_params = []
    
    for p in p_values:
        channel = QuantumDepolarizingChannel(p)
        
        def objective(params):
            prob1, theta1, phi1, theta2, phi2 = params
            prob2 = 1 - prob1
            
            state1 = create_qubit_state(theta1, phi1)
            state2 = create_qubit_state(theta2, phi2)
            
            return -channel.mutual_information([prob1, prob2], [state1, state2])
        
        # Initial guess: equal probability, orthogonal states
        initial_guess = [0.5, 0, 0, np.pi, 0]
        bounds = [(0.01, 0.99), (0, np.pi), (0, 2*np.pi), (0, np.pi), (0, 2*np.pi)]
        
        result = minimize(objective, initial_guess, bounds=bounds, method='L-BFGS-B')
        
        optimized_capacities.append(-result.fun)
        optimal_params.append(result.x)
    
    return np.array(optimized_capacities), optimal_params

# Main analysis
print("Quantum Communication Capacity Optimization Analysis")
print("=" * 55)

# Define range of depolarization parameters
p_values = np.linspace(0.001, 0.999, 50)

# Calculate theoretical capacities
theoretical_capacities = []
for p in p_values:
    channel = QuantumDepolarizingChannel(p)
    theoretical_capacities.append(channel.theoretical_capacity())

theoretical_capacities = np.array(theoretical_capacities)

# Optimize with single input state
print("Optimizing with single input state...")
single_capacities, single_thetas, single_phis = optimize_capacity_single_input(p_values)

# Optimize with two input states
print("Optimizing with two input states...")
two_capacities, two_params = optimize_capacity_two_inputs(p_values)

# Create comprehensive visualization
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
fig.suptitle('Quantum Depolarizing Channel Capacity Optimization', fontsize=16, fontweight='bold')

# Plot 1: Capacity comparison
axes[0, 0].plot(p_values, theoretical_capacities, 'k-', linewidth=3, label='Theoretical Capacity', alpha=0.8)
axes[0, 0].plot(p_values, single_capacities, 'r--', linewidth=2, label='Single Input Optimization', alpha=0.7)
axes[0, 0].plot(p_values, two_capacities, 'b:', linewidth=2, label='Two Input Optimization', alpha=0.7)
axes[0, 0].set_xlabel('Depolarization Parameter (p)')
axes[0, 0].set_ylabel('Channel Capacity (bits)')
axes[0, 0].set_title('Channel Capacity vs Depolarization Parameter')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Plot 2: Optimal single input parameters
ax2 = axes[0, 1]
ax2_twin = ax2.twinx()

line1 = ax2.plot(p_values, single_thetas, 'g-', linewidth=2, label='θ (polar angle)')
line2 = ax2_twin.plot(p_values, single_phis, 'orange', linestyle='--', linewidth=2, label='φ (azimuthal angle)')

ax2.set_xlabel('Depolarization Parameter (p)')
ax2.set_ylabel('θ (radians)', color='g')
ax2_twin.set_ylabel('φ (radians)', color='orange')
ax2.set_title('Optimal Single Input State Parameters')

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

# Plot 3: Capacity difference analysis
capacity_diff = two_capacities - single_capacities
axes[1, 0].plot(p_values, capacity_diff, 'm-', linewidth=2)
axes[1, 0].fill_between(p_values, 0, capacity_diff, alpha=0.3, color='magenta')
axes[1, 0].set_xlabel('Depolarization Parameter (p)')
axes[1, 0].set_ylabel('Capacity Difference (bits)')
axes[1, 0].set_title('Advantage of Two Input States over Single Input')
axes[1, 0].grid(True, alpha=0.3)

# Plot 4: Error analysis
theoretical_error_single = np.abs(theoretical_capacities - single_capacities)
theoretical_error_two = np.abs(theoretical_capacities - two_capacities)

axes[1, 1].semilogy(p_values, theoretical_error_single, 'r-', linewidth=2, label='Single Input Error')
axes[1, 1].semilogy(p_values, theoretical_error_two, 'b-', linewidth=2, label='Two Input Error')
axes[1, 1].set_xlabel('Depolarization Parameter (p)')
axes[1, 1].set_ylabel('Absolute Error (log scale)')
axes[1, 1].set_title('Optimization Error vs Theoretical Values')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Detailed analysis at specific points
print("\nDetailed Analysis at Key Points:")
print("-" * 40)

key_points = [0.1, 0.3, 0.5, 0.7, 0.9]
for p in key_points:
    idx = np.argmin(np.abs(p_values - p))
    print(f"\nDepolarization p = {p:.1f}:")
    print(f"  Theoretical Capacity: {theoretical_capacities[idx]:.4f} bits")
    print(f"  Single Input Capacity: {single_capacities[idx]:.4f} bits")
    print(f"  Two Input Capacity: {two_capacities[idx]:.4f} bits")
    print(f"  Optimal θ: {single_thetas[idx]:.3f} rad ({single_thetas[idx]*180/np.pi:.1f}°)")
    print(f"  Optimal φ: {single_phis[idx]:.3f} rad ({single_phis[idx]*180/np.pi:.1f}°)")

# Performance summary
print(f"\nPerformance Summary:")
print(f"  Maximum theoretical capacity: {np.max(theoretical_capacities):.4f} bits at p = {p_values[np.argmax(theoretical_capacities)]:.3f}")
print(f"  Average single input error: {np.mean(theoretical_error_single):.6f} bits")
print(f"  Average two input error: {np.mean(theoretical_error_two):.6f} bits")
print(f"  Maximum advantage of two inputs: {np.max(capacity_diff):.6f} bits")

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar, minimize
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns

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

class QKDOptimizer:
    """
    Quantum Key Distribution Rate Optimizer for BB84 Protocol
    """
    
    def __init__(self, distance=50, wavelength=1550e-9, fiber_loss=0.2):
        """
        Initialize QKD system parameters
        
        Parameters:
        - distance: transmission distance in km
        - wavelength: photon wavelength in meters
        - fiber_loss: fiber attenuation in dB/km
        """
        self.distance = distance
        self.wavelength = wavelength
        self.fiber_loss = fiber_loss
        
        # Physical constants
        self.c = 3e8  # speed of light
        self.h = 6.626e-34  # Planck constant
        
        # System parameters
        self.detector_efficiency = 0.1  # 10% detector efficiency
        self.dark_count_rate = 1e-6  # dark counts per pulse
        self.pulse_rate = 1e9  # 1 GHz pulse rate
        
    def channel_transmission(self, distance):
        """Calculate channel transmission efficiency"""
        loss_db = self.fiber_loss * distance
        return 10**(-loss_db/10)
    
    def binary_entropy(self, x):
        """Binary entropy function h(x)"""
        if x <= 0 or x >= 1:
            return 0
        return -x * np.log2(x) - (1-x) * np.log2(1-x)
    
    def qber_calculation(self, distance, detector_eff):
        """
        Calculate Quantum Bit Error Rate (QBER)
        Includes intrinsic errors and detector dark counts
        """
        eta = self.channel_transmission(distance)
        
        # Signal rate
        signal_rate = self.pulse_rate * eta * detector_eff
        
        # Dark count contribution to error rate
        dark_error_rate = self.dark_count_rate * self.pulse_rate
        
        # Intrinsic QBER (assumed 2% for perfect system)
        intrinsic_qber = 0.02
        
        # Total QBER calculation
        if signal_rate > 0:
            qber = (intrinsic_qber * signal_rate + 0.5 * dark_error_rate) / (signal_rate + dark_error_rate)
        else:
            qber = 0.5  # No signal case
            
        return min(qber, 0.5)  # QBER cannot exceed 50%
    
    def error_correction_efficiency(self, qber):
        """
        Error correction efficiency factor f(e)
        Uses typical CASCADE protocol efficiency
        """
        if qber <= 0:
            return 1.0
        return 1.22  # Typical CASCADE efficiency
    
    def key_rate(self, distance, detector_eff):
        """
        Calculate secure key generation rate
        
        R = μ * η * η_d * [1 - h(e) - f(e) * h(e)]
        
        Where:
        - μ: pulse rate
        - η: channel transmission
        - η_d: detector efficiency  
        - e: QBER
        - h(e): binary entropy
        - f(e): error correction efficiency
        """
        eta = self.channel_transmission(distance)
        qber = self.qber_calculation(distance, detector_eff)
        
        # Calculate quantum bit rate
        q_rate = self.pulse_rate * eta * detector_eff
        
        # Check if QBER is below theoretical limit (11% for BB84)
        if qber > 0.11:
            return 0  # No secure key can be generated
        
        # Calculate key rate components
        h_e = self.binary_entropy(qber)
        f_e = self.error_correction_efficiency(qber)
        
        # Secure key rate formula
        rate = q_rate * (1 - h_e - f_e * h_e)
        
        return max(0, rate)  # Rate cannot be negative
    
    def optimize_detector_efficiency(self, distance):
        """Optimize detector efficiency for maximum key rate at given distance"""
        
        def objective(detector_eff):
            return -self.key_rate(distance, detector_eff[0])
        
        # Detector efficiency must be between 0.01% and 90%
        bounds = [(0.0001, 0.9)]
        
        result = minimize(objective, x0=[0.1], bounds=bounds, method='L-BFGS-B')
        
        optimal_eff = result.x[0]
        max_rate = -result.fun
        
        return optimal_eff, max_rate
    
    def analyze_system_performance(self):
        """Comprehensive analysis of QKD system performance"""
        
        # Distance range for analysis
        distances = np.linspace(1, 200, 100)
        
        # Detector efficiency range
        detector_effs = np.linspace(0.01, 0.5, 50)
        
        # Calculate key rates for different parameters
        key_rates_distance = []
        qber_values = []
        optimal_efficiencies = []
        max_rates = []
        
        print("Analyzing QKD system performance...")
        print("=" * 50)
        
        for dist in distances:
            rate = self.key_rate(dist, self.detector_efficiency)
            qber = self.qber_calculation(dist, self.detector_efficiency)
            
            key_rates_distance.append(rate)
            qber_values.append(qber)
            
            # Find optimal detector efficiency for this distance
            opt_eff, max_rate = self.optimize_detector_efficiency(dist)
            optimal_efficiencies.append(opt_eff)
            max_rates.append(max_rate)
        
        # Create 2D parameter space analysis
        Distance, DetectorEff = np.meshgrid(distances[:50], detector_effs)
        KeyRates2D = np.zeros_like(Distance)
        
        for i, det_eff in enumerate(detector_effs):
            for j, dist in enumerate(distances[:50]):
                KeyRates2D[i, j] = self.key_rate(dist, det_eff)
        
        return {
            'distances': distances,
            'key_rates_distance': key_rates_distance,
            'qber_values': qber_values,
            'optimal_efficiencies': optimal_efficiencies,
            'max_rates': max_rates,
            'detector_effs': detector_effs,
            'Distance': Distance,
            'DetectorEff': DetectorEff,
            'KeyRates2D': KeyRates2D
        }

# Initialize QKD system
qkd = QKDOptimizer(distance=50, fiber_loss=0.2)

# Perform comprehensive analysis
results = qkd.analyze_system_performance()

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

# Plot 1: Key rate vs Distance
ax1 = plt.subplot(2, 3, 1)
plt.plot(results['distances'], np.array(results['key_rates_distance']) * 1e-6, 'b-', linewidth=2, label='Standard System')
plt.plot(results['distances'], np.array(results['max_rates']) * 1e-6, 'r--', linewidth=2, label='Optimized System')
plt.xlabel('Distance (km)')
plt.ylabel('Key Rate (Mbps)')
plt.title('QKD Key Rate vs Transmission Distance')
plt.legend()
plt.grid(True, alpha=0.3)
plt.xlim(0, 200)

# Plot 2: QBER vs Distance
ax2 = plt.subplot(2, 3, 2)
plt.plot(results['distances'], np.array(results['qber_values']) * 100, 'g-', linewidth=2)
plt.axhline(y=11, color='r', linestyle='--', alpha=0.7, label='Security Threshold (11%)')
plt.xlabel('Distance (km)')
plt.ylabel('QBER (%)')
plt.title('Quantum Bit Error Rate vs Distance')
plt.legend()
plt.grid(True, alpha=0.3)
plt.xlim(0, 200)

# Plot 3: Optimal Detector Efficiency vs Distance
ax3 = plt.subplot(2, 3, 3)
plt.plot(results['distances'], np.array(results['optimal_efficiencies']) * 100, 'm-', linewidth=2)
plt.xlabel('Distance (km)')
plt.ylabel('Optimal Detector Efficiency (%)')
plt.title('Optimal Detector Efficiency vs Distance')
plt.grid(True, alpha=0.3)
plt.xlim(0, 200)

# Plot 4: 3D Parameter Space
ax4 = plt.subplot(2, 3, 4, projection='3d')
surf = ax4.plot_surface(results['Distance'], results['DetectorEff'] * 100, 
                       results['KeyRates2D'] * 1e-6, cmap='viridis', alpha=0.8)
ax4.set_xlabel('Distance (km)')
ax4.set_ylabel('Detector Efficiency (%)')
ax4.set_zlabel('Key Rate (Mbps)')
ax4.set_title('3D Parameter Space Analysis')
plt.colorbar(surf, ax=ax4, shrink=0.5)

# Plot 5: Heatmap of parameter space
ax5 = plt.subplot(2, 3, 5)
heatmap_data = results['KeyRates2D'] * 1e-6
im = ax5.imshow(heatmap_data, aspect='auto', origin='lower', cmap='plasma')
ax5.set_xlabel('Distance Index')
ax5.set_ylabel('Detector Efficiency Index')
ax5.set_title('Key Rate Heatmap (Mbps)')
plt.colorbar(im, ax=ax5)

# Plot 6: Optimization comparison
ax6 = plt.subplot(2, 3, 6)
improvement = (np.array(results['max_rates']) - np.array(results['key_rates_distance'])) / np.array(results['key_rates_distance']) * 100
valid_improvement = improvement[np.array(results['key_rates_distance']) > 0]
valid_distances = results['distances'][np.array(results['key_rates_distance']) > 0]

plt.plot(valid_distances, valid_improvement, 'purple', linewidth=2)
plt.xlabel('Distance (km)')
plt.ylabel('Rate Improvement (%)')
plt.title('Optimization Benefit vs Distance')
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Print optimization results
print("\nQKD System Optimization Results")
print("=" * 50)

# Find maximum transmission distance
max_distance_idx = np.where(np.array(results['key_rates_distance']) > 0)[0]
if len(max_distance_idx) > 0:
    max_distance = results['distances'][max_distance_idx[-1]]
    print(f"Maximum transmission distance: {max_distance:.1f} km")
else:
    print("No secure transmission possible with current parameters")

# Find optimal operating point
max_rate_idx = np.argmax(results['max_rates'])
optimal_distance = results['distances'][max_rate_idx]
optimal_rate = results['max_rates'][max_rate_idx]
optimal_detector_eff = results['optimal_efficiencies'][max_rate_idx]

print(f"Optimal operating distance: {optimal_distance:.1f} km")
print(f"Maximum achievable key rate: {optimal_rate*1e-6:.3f} Mbps")
print(f"Optimal detector efficiency: {optimal_detector_eff*100:.1f}%")

# Calculate specific examples
print(f"\nSpecific Examples:")
print("-" * 30)

test_distances = [10, 50, 100, 150]
for dist in test_distances:
    standard_rate = qkd.key_rate(dist, qkd.detector_efficiency)
    opt_eff, opt_rate = qkd.optimize_detector_efficiency(dist)
    qber = qkd.qber_calculation(dist, qkd.detector_efficiency)
    
    if standard_rate > 0:
        improvement = (opt_rate - standard_rate) / standard_rate * 100
        print(f"Distance {dist} km:")
        print(f"  Standard rate: {standard_rate*1e-6:.4f} Mbps")
        print(f"  Optimized rate: {opt_rate*1e-6:.4f} Mbps")
        print(f"  Improvement: {improvement:.1f}%")
        print(f"  QBER: {qber*100:.2f}%")
        print(f"  Optimal detector efficiency: {opt_eff*100:.1f}%")
    else:
        print(f"Distance {dist} km: No secure transmission possible")
    print()