import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from matplotlib.patches import Rectangle
import seaborn as sns

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

class AntennaArrayOptimizer:
    """
    Class for optimizing electromagnetic field distribution in antenna arrays
    """
    
    def __init__(self, num_elements=8, element_spacing=0.5, frequency=2.4e9):
        """
        Initialize the antenna array optimizer
        
        Parameters:
        - num_elements: Number of antenna elements
        - element_spacing: Spacing between elements in wavelengths
        - frequency: Operating frequency in Hz
        """
        self.num_elements = num_elements
        self.element_spacing = element_spacing  # in wavelengths
        self.frequency = frequency
        self.c = 3e8  # Speed of light
        self.wavelength = self.c / self.frequency
        self.k = 2 * np.pi / self.wavelength  # Wave number
        
        # Element positions (in wavelengths)
        self.element_positions = np.arange(num_elements) * element_spacing
        
    def array_factor(self, theta, amplitudes, phases):
        """
        Calculate the array factor for given amplitudes and phases
        
        Array factor formula:
        AF(θ) = Σ(n=0 to N-1) I_n * exp(j * (k * d_n * cos(θ) + φ_n))
        
        Parameters:
        - theta: Angle array in radians
        - amplitudes: Current amplitudes for each element
        - phases: Phase shifts for each element
        
        Returns:
        - Complex array factor
        """
        theta = np.atleast_1d(theta)
        AF = np.zeros_like(theta, dtype=complex)
        
        for n in range(self.num_elements):
            # Phase contribution from element position and applied phase
            phase_term = self.k * self.element_positions[n] * self.wavelength * np.cos(theta) + phases[n]
            AF += amplitudes[n] * np.exp(1j * phase_term)
            
        return AF
    
    def desired_pattern(self, theta, beam_direction=0, sidelobe_level=-20):
        """
        Define the desired radiation pattern
        
        Parameters:
        - theta: Angle array in radians
        - beam_direction: Main beam direction in radians
        - sidelobe_level: Desired sidelobe level in dB
        """
        # Simple desired pattern: high gain in beam direction, low elsewhere
        desired = np.ones_like(theta) * 10**(sidelobe_level/20)
        
        # Main beam region (±30 degrees around beam direction)
        beam_width = np.pi/6  # 30 degrees
        main_beam_mask = np.abs(theta - beam_direction) <= beam_width
        desired[main_beam_mask] = 1.0
        
        return desired
    
    def cost_function(self, x, theta_range, desired_pattern):
        """
        Cost function to minimize
        
        The optimization variables x contain:
        - x[0:N]: Amplitudes for each element
        - x[N:2N]: Phases for each element
        """
        N = self.num_elements
        amplitudes = x[:N]
        phases = x[N:2*N]
        
        # Calculate array factor
        AF = self.array_factor(theta_range, amplitudes, phases)
        actual_pattern = np.abs(AF)
        
        # Normalize patterns
        actual_pattern = actual_pattern / np.max(actual_pattern)
        desired_norm = desired_pattern / np.max(desired_pattern)
        
        # Cost function: Mean squared error + regularization
        mse = np.mean((actual_pattern - desired_norm)**2)
        
        # Add penalty for very small amplitudes to avoid numerical issues
        amplitude_penalty = np.sum(1.0 / (amplitudes + 1e-6))
        
        return mse + 0.001 * amplitude_penalty
    
    def optimize_pattern(self, beam_direction=0, sidelobe_level=-20):
        """
        Optimize the antenna array pattern
        """
        # Define angle range for optimization
        theta_range = np.linspace(0, 2*np.pi, 360)
        desired = self.desired_pattern(theta_range, beam_direction, sidelobe_level)
        
        # Initial guess: uniform amplitude, zero phase
        x0 = np.ones(2 * self.num_elements)
        x0[self.num_elements:] = 0  # Zero initial phases
        
        # Constraints: amplitudes should be positive
        bounds = [(0.1, 2.0) for _ in range(self.num_elements)] + \
                 [(-np.pi, np.pi) for _ in range(self.num_elements)]
        
        # Perform optimization
        result = minimize(
            self.cost_function,
            x0,
            args=(theta_range, desired),
            method='L-BFGS-B',
            bounds=bounds,
            options={'maxiter': 1000}
        )
        
        return result, theta_range, desired
    
    def plot_results(self, result, theta_range, desired_pattern, beam_direction=0):
        """
        Create comprehensive plots of the optimization results
        """
        # Extract optimized parameters
        N = self.num_elements
        opt_amplitudes = result.x[:N]
        opt_phases = result.x[N:2*N]
        
        # Calculate patterns
        AF_initial = self.array_factor(theta_range, np.ones(N), np.zeros(N))
        AF_optimized = self.array_factor(theta_range, opt_amplitudes, opt_phases)
        
        # Create figure with subplots
        fig = plt.figure(figsize=(16, 12))
        
        # 1. Radiation Pattern Comparison (Linear Scale)
        ax1 = plt.subplot(2, 3, 1)
        theta_deg = np.degrees(theta_range)
        
        plt.plot(theta_deg, np.abs(AF_initial)/np.max(np.abs(AF_initial)), 
                'b--', label='Initial (Uniform)', linewidth=2)
        plt.plot(theta_deg, np.abs(AF_optimized)/np.max(np.abs(AF_optimized)), 
                'r-', label='Optimized', linewidth=2)
        plt.plot(theta_deg, desired_pattern/np.max(desired_pattern), 
                'g:', label='Desired', linewidth=2, alpha=0.7)
        
        plt.axvline(x=np.degrees(beam_direction), color='k', linestyle=':', alpha=0.5, label='Beam Direction')
        plt.xlabel('Angle (degrees)')
        plt.ylabel('Normalized Amplitude')
        plt.title('Radiation Pattern Comparison (Linear)')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.xlim(0, 360)
        
        # 2. Radiation Pattern in dB
        ax2 = plt.subplot(2, 3, 2)
        initial_db = 20 * np.log10(np.abs(AF_initial)/np.max(np.abs(AF_initial)) + 1e-10)
        optimized_db = 20 * np.log10(np.abs(AF_optimized)/np.max(np.abs(AF_optimized)) + 1e-10)
        
        plt.plot(theta_deg, initial_db, 'b--', label='Initial', linewidth=2)
        plt.plot(theta_deg, optimized_db, 'r-', label='Optimized', linewidth=2)
        plt.axvline(x=np.degrees(beam_direction), color='k', linestyle=':', alpha=0.5)
        
        plt.xlabel('Angle (degrees)')
        plt.ylabel('Gain (dB)')
        plt.title('Radiation Pattern (dB Scale)')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.ylim(-50, 5)
        plt.xlim(0, 360)
        
        # 3. Polar Plot
        ax3 = plt.subplot(2, 3, 3, projection='polar')
        optimized_norm = np.abs(AF_optimized)/np.max(np.abs(AF_optimized))
        
        ax3.plot(theta_range, optimized_norm, 'r-', linewidth=2, label='Optimized')
        ax3.plot(theta_range, np.abs(AF_initial)/np.max(np.abs(AF_initial)), 
                'b--', linewidth=2, alpha=0.7, label='Initial')
        
        ax3.set_title('Polar Radiation Pattern')
        ax3.legend(loc='upper right', bbox_to_anchor=(1.3, 1.0))
        ax3.grid(True)
        
        # 4. Element Amplitudes
        ax4 = plt.subplot(2, 3, 4)
        element_numbers = np.arange(1, N+1)
        
        plt.bar(element_numbers, opt_amplitudes, alpha=0.7, color='orange', 
                edgecolor='black', linewidth=1.5)
        plt.xlabel('Element Number')
        plt.ylabel('Amplitude')
        plt.title('Optimized Element Amplitudes')
        plt.grid(True, alpha=0.3)
        
        # Add values on top of bars
        for i, v in enumerate(opt_amplitudes):
            plt.text(i+1, v+0.02, f'{v:.2f}', ha='center', va='bottom', fontweight='bold')
        
        # 5. Element Phases
        ax5 = plt.subplot(2, 3, 5)
        phases_deg = np.degrees(opt_phases)
        
        plt.bar(element_numbers, phases_deg, alpha=0.7, color='lightblue', 
                edgecolor='black', linewidth=1.5)
        plt.xlabel('Element Number')
        plt.ylabel('Phase (degrees)')
        plt.title('Optimized Element Phases')
        plt.grid(True, alpha=0.3)
        
        # Add values on top of bars
        for i, v in enumerate(phases_deg):
            plt.text(i+1, v+2, f'{v:.1f}°', ha='center', va='bottom', fontweight='bold')
        
        # 6. Convergence Information
        ax6 = plt.subplot(2, 3, 6)
        ax6.axis('off')
        
        # Display optimization results
        info_text = f"""
        Optimization Results:
        
        • Number of Elements: {N}
        • Element Spacing: {self.element_spacing}λ
        • Frequency: {self.frequency/1e9:.1f} GHz
        • Beam Direction: {np.degrees(beam_direction):.1f}°
        
        • Final Cost: {result.fun:.6f}
        • Iterations: {result.nit}
        • Success: {result.success}
        
        Performance Metrics:
        • Max Amplitude: {np.max(opt_amplitudes):.3f}
        • Min Amplitude: {np.min(opt_amplitudes):.3f}
        • Phase Range: {np.degrees(np.max(opt_phases) - np.min(opt_phases)):.1f}°
        """
        
        ax6.text(0.1, 0.9, info_text, transform=ax6.transAxes, fontsize=11,
                verticalalignment='top', fontfamily='monospace',
                bbox=dict(boxstyle='round', facecolor='lightgray', alpha=0.8))
        
        plt.tight_layout()
        return fig

# Example usage and demonstration
def demonstrate_optimization():
    """
    Demonstrate the electromagnetic field optimization with different scenarios
    """
    print("🔬 Electromagnetic Field Distribution Optimization Demo")
    print("=" * 60)
    
    # Create optimizer instance
    optimizer = AntennaArrayOptimizer(num_elements=8, element_spacing=0.5, frequency=2.4e9)
    
    # Scenario 1: Main beam at 0 degrees (broadside)
    print("\n📡 Scenario 1: Broadside Array (0°)")
    result1, theta_range, desired1 = optimizer.optimize_pattern(beam_direction=0, sidelobe_level=-20)
    fig1 = optimizer.plot_results(result1, theta_range, desired1, beam_direction=0)
    plt.suptitle('Broadside Array Optimization (0° beam)', fontsize=16, fontweight='bold')
    plt.show()
    
    # Scenario 2: Steered beam at 30 degrees
    print("\n📡 Scenario 2: Steered Array (30°)")
    beam_angle = np.pi/6  # 30 degrees
    result2, theta_range, desired2 = optimizer.optimize_pattern(beam_direction=beam_angle, sidelobe_level=-25)
    fig2 = optimizer.plot_results(result2, theta_range, desired2, beam_direction=beam_angle)
    plt.suptitle('Steered Array Optimization (30° beam)', fontsize=16, fontweight='bold')
    plt.show()
    
    # Performance comparison
    print("\n📊 Performance Analysis:")
    print("-" * 40)
    
    N = optimizer.num_elements
    opt_amps1 = result1.x[:N]
    opt_phases1 = result1.x[N:2*N]
    opt_amps2 = result2.x[:N]
    opt_phases2 = result2.x[N:2*N]
    
    print(f"Broadside Array:")
    print(f"  • Final cost: {result1.fun:.6f}")
    print(f"  • Amplitude variation: {np.std(opt_amps1):.3f}")
    print(f"  • Phase variation: {np.degrees(np.std(opt_phases1)):.1f}°")
    
    print(f"\nSteered Array:")
    print(f"  • Final cost: {result2.fun:.6f}")
    print(f"  • Amplitude variation: {np.std(opt_amps2):.3f}")
    print(f"  • Phase variation: {np.degrees(np.std(opt_phases2)):.1f}°")
    
    # Calculate directivity improvement
    AF_initial = optimizer.array_factor(theta_range, np.ones(N), np.zeros(N))
    AF_opt1 = optimizer.array_factor(theta_range, opt_amps1, opt_phases1)
    
    directivity_initial = np.max(np.abs(AF_initial)**2) / np.mean(np.abs(AF_initial)**2)
    directivity_opt = np.max(np.abs(AF_opt1)**2) / np.mean(np.abs(AF_opt1)**2)
    
    print(f"\nDirectivity Improvement:")
    print(f"  • Initial: {10*np.log10(directivity_initial):.1f} dB")
    print(f"  • Optimized: {10*np.log10(directivity_opt):.1f} dB")
    print(f"  • Improvement: {10*np.log10(directivity_opt/directivity_initial):.1f} dB")

if __name__ == "__main__":
    demonstrate_optimization()

import numpy as np
import matplotlib.pyplot as plt
import random
from matplotlib.animation import FuncAnimation
import seaborn as sns

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

class SimulatedAnnealingTSP:
    """
    Simulated Annealing implementation for Traveling Salesman Problem
    Demonstrates phase transition optimization behavior
    """
    
    def __init__(self, cities, initial_temp=1000, cooling_rate=0.995, min_temp=0.01):
        """
        Initialize the SA algorithm
        
        Parameters:
        cities: array of city coordinates [(x1,y1), (x2,y2), ...]
        initial_temp: Starting temperature T_0
        cooling_rate: Cooling factor α (T_new = α * T_old)
        min_temp: Minimum temperature to stop
        """
        self.cities = np.array(cities)
        self.n_cities = len(cities)
        self.initial_temp = initial_temp
        self.cooling_rate = cooling_rate
        self.min_temp = min_temp
        
        # Distance matrix calculation
        self.distance_matrix = self._calculate_distance_matrix()
        
        # Tracking variables
        self.temperatures = []
        self.energies = []
        self.best_energies = []
        self.acceptance_rates = []
        
    def _calculate_distance_matrix(self):
        """Calculate Euclidean distance matrix between all cities"""
        n = self.n_cities
        dist_matrix = np.zeros((n, n))
        for i in range(n):
            for j in range(n):
                if i != j:
                    dist_matrix[i][j] = np.sqrt(
                        (self.cities[i][0] - self.cities[j][0])**2 + 
                        (self.cities[i][1] - self.cities[j][1])**2
                    )
        return dist_matrix
    
    def _calculate_total_distance(self, tour):
        """
        Calculate total tour distance
        Energy function: E(tour) = Σ d(city_i, city_(i+1))
        """
        total_distance = 0
        for i in range(len(tour)):
            from_city = tour[i]
            to_city = tour[(i + 1) % len(tour)]
            total_distance += self.distance_matrix[from_city][to_city]
        return total_distance
    
    def _generate_neighbor(self, tour):
        """
        Generate neighbor solution using 2-opt swap
        This is a local search move that reverses a segment of the tour
        """
        new_tour = tour.copy()
        i, j = sorted(random.sample(range(len(tour)), 2))
        # Reverse the segment between i and j
        new_tour[i:j+1] = new_tour[i:j+1][::-1]
        return new_tour
    
    def _acceptance_probability(self, current_energy, new_energy, temperature):
        """
        Calculate acceptance probability using Boltzmann distribution
        P(accept) = exp(-ΔE / T) if ΔE > 0, else 1
        """
        if new_energy < current_energy:
            return 1.0
        if temperature == 0:
            return 0.0
        return np.exp(-(new_energy - current_energy) / temperature)
    
    def solve(self, max_iterations=50000):
        """
        Main simulated annealing algorithm
        Demonstrates phase transition from exploration to exploitation
        """
        # Initialize with random tour
        current_tour = list(range(self.n_cities))
        random.shuffle(current_tour)
        current_energy = self._calculate_total_distance(current_tour)
        
        # Best solution tracking
        best_tour = current_tour.copy()
        best_energy = current_energy
        
        # Temperature initialization
        temperature = self.initial_temp
        
        # Statistics tracking
        accepted_moves = 0
        total_moves = 0
        
        print("Starting Simulated Annealing optimization...")
        print(f"Initial energy: {current_energy:.2f}")
        
        iteration = 0
        while temperature > self.min_temp and iteration < max_iterations:
            # Generate neighbor solution
            new_tour = self._generate_neighbor(current_tour)
            new_energy = self._calculate_total_distance(new_tour)
            
            # Calculate acceptance probability
            accept_prob = self._acceptance_probability(current_energy, new_energy, temperature)
            
            # Accept or reject the move
            if random.random() < accept_prob:
                current_tour = new_tour
                current_energy = new_energy
                accepted_moves += 1
                
                # Update best solution
                if current_energy < best_energy:
                    best_tour = current_tour.copy()
                    best_energy = current_energy
            
            total_moves += 1
            
            # Record statistics every 100 iterations
            if iteration % 100 == 0:
                self.temperatures.append(temperature)
                self.energies.append(current_energy)
                self.best_energies.append(best_energy)
                
                # Calculate acceptance rate for last 100 moves
                if total_moves >= 100:
                    acceptance_rate = accepted_moves / min(100, total_moves)
                    self.acceptance_rates.append(acceptance_rate)
                else:
                    self.acceptance_rates.append(accepted_moves / total_moves if total_moves > 0 else 0)
                
                accepted_moves = 0  # Reset counter
            
            # Cool down the temperature (exponential cooling schedule)
            temperature *= self.cooling_rate
            iteration += 1
        
        print(f"Optimization completed after {iteration} iterations")
        print(f"Final temperature: {temperature:.6f}")
        print(f"Best energy found: {best_energy:.2f}")
        
        return best_tour, best_energy
    
    def visualize_results(self, best_tour):
        """
        Create comprehensive visualization of the optimization process
        """
        fig = plt.figure(figsize=(16, 12))
        
        # 1. Temperature decay (Phase transition control parameter)
        ax1 = plt.subplot(2, 3, 1)
        plt.plot(self.temperatures, 'b-', linewidth=2, alpha=0.8)
        plt.title('Temperature Decay\n(Phase Transition Control)', fontsize=12, fontweight='bold')
        plt.xlabel('Iteration (×100)')
        plt.ylabel('Temperature T')
        plt.yscale('log')
        plt.grid(True, alpha=0.3)
        
        # Add phase regions
        high_temp_region = len(self.temperatures) // 3
        plt.axvspan(0, high_temp_region, alpha=0.2, color='red', label='High T: Exploration')
        plt.axvspan(high_temp_region, len(self.temperatures), alpha=0.2, color='blue', label='Low T: Exploitation')
        plt.legend()
        
        # 2. Energy evolution (Objective function)
        ax2 = plt.subplot(2, 3, 2)
        plt.plot(self.energies, 'g-', alpha=0.6, label='Current Energy', linewidth=1)
        plt.plot(self.best_energies, 'r-', linewidth=2, label='Best Energy')
        plt.title('Energy Evolution\nE(tour) = Σ distances', fontsize=12, fontweight='bold')
        plt.xlabel('Iteration (×100)')
        plt.ylabel('Total Distance')
        plt.legend()
        plt.grid(True, alpha=0.3)
        
        # 3. Acceptance rate (Phase transition indicator)
        ax3 = plt.subplot(2, 3, 3)
        plt.plot(self.acceptance_rates, 'purple', linewidth=2)
        plt.title('Acceptance Rate\nP(accept) = exp(-ΔE/T)', fontsize=12, fontweight='bold')
        plt.xlabel('Iteration (×100)')
        plt.ylabel('Acceptance Rate')
        plt.ylim(0, 1)
        plt.grid(True, alpha=0.3)
        
        # Add horizontal lines for reference
        plt.axhline(y=0.5, color='orange', linestyle='--', alpha=0.7, label='50% acceptance')
        plt.axhline(y=0.1, color='red', linestyle='--', alpha=0.7, label='10% acceptance')
        plt.legend()
        
        # 4. Phase space plot (Temperature vs Energy)
        ax4 = plt.subplot(2, 3, 4)
        scatter = plt.scatter(self.temperatures, self.energies, 
                            c=range(len(self.temperatures)), 
                            cmap='viridis', alpha=0.7, s=20)
        plt.xlabel('Temperature T')
        plt.ylabel('Energy E')
        plt.title('Phase Space\n(T vs E trajectory)', fontsize=12, fontweight='bold')
        plt.xscale('log')
        cbar = plt.colorbar(scatter)
        cbar.set_label('Iteration')
        plt.grid(True, alpha=0.3)
        
        # 5. Best tour visualization
        ax5 = plt.subplot(2, 3, 5)
        
        # Plot cities
        city_x = [self.cities[i][0] for i in best_tour]
        city_y = [self.cities[i][1] for i in best_tour]
        
        # Close the tour
        city_x.append(city_x[0])
        city_y.append(city_y[0])
        
        plt.plot(city_x, city_y, 'bo-', linewidth=2, markersize=8, alpha=0.8)
        plt.scatter(self.cities[:, 0], self.cities[:, 1], c='red', s=100, zorder=5)
        
        # Add city labels
        for i, (x, y) in enumerate(self.cities):
            plt.annotate(str(i), (x, y), xytext=(5, 5), textcoords='offset points',
                        fontsize=8, fontweight='bold')
        
        plt.title(f'Optimal Tour\nTotal Distance: {self.best_energies[-1]:.2f}', 
                 fontsize=12, fontweight='bold')
        plt.xlabel('X Coordinate')
        plt.ylabel('Y Coordinate')
        plt.grid(True, alpha=0.3)
        plt.axis('equal')
        
        # 6. Convergence analysis
        ax6 = plt.subplot(2, 3, 6)
        
        # Calculate improvement rate
        improvements = []
        for i in range(1, len(self.best_energies)):
            improvement = (self.best_energies[i-1] - self.best_energies[i]) / self.best_energies[i-1] * 100
            improvements.append(max(0, improvement))
        
        plt.plot(improvements, 'orange', linewidth=2)
        plt.title('Convergence Analysis\nImprovement Rate (%)', fontsize=12, fontweight='bold')
        plt.xlabel('Iteration (×100)')
        plt.ylabel('Improvement Rate (%)')
        plt.grid(True, alpha=0.3)
        
        # Add phase transition markers
        if len(improvements) > 10:
            transition_point = len(improvements) // 3
            plt.axvline(x=transition_point, color='red', linestyle='--', alpha=0.7, 
                       label='Exploration→Exploitation')
            plt.legend()
        
        plt.tight_layout()
        plt.show()
        
        return fig

# Generate random cities for TSP instance
def generate_random_cities(n_cities, seed=42):
    """Generate n_cities random city coordinates"""
    np.random.seed(seed)
    random.seed(seed)
    
    cities = []
    for i in range(n_cities):
        x = np.random.uniform(0, 100)
        y = np.random.uniform(0, 100)
        cities.append((x, y))
    
    return cities

# Main execution
if __name__ == "__main__":
    print("=== Phase Transition Optimization: Simulated Annealing for TSP ===\n")
    
    # Problem setup
    n_cities = 20
    cities = generate_random_cities(n_cities)
    
    print(f"Generated {n_cities} random cities:")
    for i, (x, y) in enumerate(cities[:5]):  # Show first 5 cities
        print(f"City {i}: ({x:.2f}, {y:.2f})")
    print("...\n")
    
    # Initialize and run simulated annealing
    sa = SimulatedAnnealingTSP(
        cities=cities,
        initial_temp=1000.0,  # High initial temperature for exploration
        cooling_rate=0.995,   # Slow cooling for gradual phase transition
        min_temp=0.01         # Low final temperature for exploitation
    )
    
    # Solve the problem
    best_tour, best_distance = sa.solve(max_iterations=50000)
    
    # Display results
    print(f"\n=== RESULTS ===")
    print(f"Best tour found: {best_tour}")
    print(f"Best distance: {best_distance:.2f}")
    print(f"Number of temperature steps: {len(sa.temperatures)}")
    
    # Create comprehensive visualization
    print("\nGenerating visualization...")
    sa.visualize_results(best_tour)
    
    # Additional analysis: Phase transition characteristics
    print(f"\n=== PHASE TRANSITION ANALYSIS ===")
    
    # Calculate phase transition point (where acceptance rate drops significantly)
    if len(sa.acceptance_rates) > 10:
        acceptance_array = np.array(sa.acceptance_rates)
        # Find where acceptance rate first drops below 20%
        transition_indices = np.where(acceptance_array < 0.2)[0]
        if len(transition_indices) > 0:
            transition_point = transition_indices[0]
            transition_temp = sa.temperatures[transition_point]
            print(f"Phase transition occurs around iteration {transition_point * 100}")
            print(f"Transition temperature: {transition_temp:.2f}")
            print(f"At transition - Acceptance rate: {sa.acceptance_rates[transition_point]:.2%}")
    
    # Energy landscape analysis
    initial_energy = sa.energies[0] if sa.energies else 0
    final_energy = sa.best_energies[-1] if sa.best_energies else 0
    improvement = (initial_energy - final_energy) / initial_energy * 100
    
    print(f"Total improvement: {improvement:.1f}%")
    print(f"Final acceptance rate: {sa.acceptance_rates[-1]:.2%}")
    
    print("\n=== MATHEMATICAL FORMULATION ===")
    print("Objective function: f(π) = Σ d(πᵢ, π₍ᵢ₊₁₎ mod n)")
    print("Acceptance probability: P(ΔE) = exp(-ΔE/T)")
    print("Cooling schedule: T(k+1) = α × T(k), where α = 0.995")
    print("Phase transition: High T (exploration) → Low T (exploitation)")

def _acceptance_probability(self, current_energy, new_energy, temperature):
    if new_energy < current_energy:
        return 1.0  # Always accept improvements
    if temperature == 0:
        return 0.0  # Never accept worse solutions at T=0
    return np.exp(-(new_energy - current_energy) / temperature)

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

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

class CoolingSystemOptimizer:
    """
    A class to optimize data center cooling systems
    """
    
    def __init__(self):
        # System parameters
        self.alpha = 1.2  # Energy cost coefficient
        self.beta = 0.8   # Flow rate cost coefficient  
        self.gamma = 2.5  # Temperature deviation cost coefficient
        self.T_ambient = 25.0  # Ambient temperature (°C)
        
        # Constraints
        self.T_min = 18.0  # Minimum cooling temperature (°C)
        self.T_max = 28.0  # Maximum cooling temperature (°C)
        self.F_min = 10.0  # Minimum flow rate (m³/min)
        self.F_max = 100.0 # Maximum flow rate (m³/min)
        
    def server_power_factor(self, T_c):
        """
        Calculate server power factor based on cooling temperature
        Higher temperatures lead to higher power consumption
        """
        return 1.0 + 0.05 * np.exp((T_c - 20) / 10)
    
    def cooling_energy_cost(self, T_c):
        """
        Calculate energy cost due to server power consumption
        """
        return self.alpha * 1000 * self.server_power_factor(T_c)  # Base server power: 1000W
    
    def cooling_system_cost(self, T_c, F):
        """
        Calculate cooling system operational cost
        """
        flow_cost = self.beta * F**2
        temp_deviation_cost = self.gamma * (self.T_ambient - T_c)**2
        return flow_cost + temp_deviation_cost
    
    def total_cost(self, params):
        """
        Calculate total system cost
        params: [T_c, F] - cooling temperature and flow rate
        """
        T_c, F = params
        energy_cost = self.cooling_energy_cost(T_c)
        system_cost = self.cooling_system_cost(T_c, F)
        return energy_cost + system_cost
    
    def constraints(self):
        """
        Define optimization constraints
        """
        return [
            {'type': 'ineq', 'fun': lambda x: x[0] - self.T_min},  # T_c >= T_min
            {'type': 'ineq', 'fun': lambda x: self.T_max - x[0]},  # T_c <= T_max
            {'type': 'ineq', 'fun': lambda x: x[1] - self.F_min},  # F >= F_min
            {'type': 'ineq', 'fun': lambda x: self.F_max - x[1]}   # F <= F_max
        ]
    
    def optimize(self):
        """
        Perform optimization using scipy.optimize.minimize
        """
        # Initial guess: middle values
        x0 = [(self.T_min + self.T_max) / 2, (self.F_min + self.F_max) / 2]
        
        # Optimize
        result = minimize(
            self.total_cost,
            x0,
            method='SLSQP',
            constraints=self.constraints(),
            options={'disp': True}
        )
        
        return result
    
    def analyze_sensitivity(self, optimal_params):
        """
        Perform sensitivity analysis around optimal point
        """
        T_c_opt, F_opt = optimal_params
        
        # Temperature sensitivity
        T_range = np.linspace(self.T_min, self.T_max, 50)
        costs_T = [self.total_cost([T, F_opt]) for T in T_range]
        
        # Flow rate sensitivity  
        F_range = np.linspace(self.F_min, self.F_max, 50)
        costs_F = [self.total_cost([T_c_opt, F]) for F in F_range]
        
        return T_range, costs_T, F_range, costs_F
    
    def create_cost_surface(self):
        """
        Create 3D cost surface for visualization
        """
        T_grid = np.linspace(self.T_min, self.T_max, 30)
        F_grid = np.linspace(self.F_min, self.F_max, 30)
        T_mesh, F_mesh = np.meshgrid(T_grid, F_grid)
        
        Cost_mesh = np.zeros_like(T_mesh)
        for i in range(T_mesh.shape[0]):
            for j in range(T_mesh.shape[1]):
                Cost_mesh[i, j] = self.total_cost([T_mesh[i, j], F_mesh[i, j]])
        
        return T_mesh, F_mesh, Cost_mesh

# Initialize optimizer
optimizer = CoolingSystemOptimizer()

print("=== Cooling System Optimization ===")
print(f"Temperature range: {optimizer.T_min}°C - {optimizer.T_max}°C")
print(f"Flow rate range: {optimizer.F_min} - {optimizer.F_max} m³/min")
print(f"Ambient temperature: {optimizer.T_ambient}°C")
print()

# Perform optimization
result = optimizer.optimize()

print("\n=== Optimization Results ===")
if result.success:
    T_c_optimal, F_optimal = result.x
    print(f"Optimal cooling temperature: {T_c_optimal:.2f}°C")
    print(f"Optimal flow rate: {F_optimal:.2f} m³/min")
    print(f"Minimum total cost: ${result.fun:.2f}")
    
    # Cost breakdown
    energy_cost = optimizer.cooling_energy_cost(T_c_optimal)
    system_cost = optimizer.cooling_system_cost(T_c_optimal, F_optimal)
    print(f"Energy cost: ${energy_cost:.2f}")
    print(f"Cooling system cost: ${system_cost:.2f}")
else:
    print("Optimization failed!")
    print(result.message)

# Sensitivity analysis
T_range, costs_T, F_range, costs_F = optimizer.analyze_sensitivity(result.x)

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

# 1. Cost vs Temperature (Flow rate fixed at optimal)
ax1 = plt.subplot(2, 3, 1)
plt.plot(T_range, costs_T, 'b-', linewidth=2, label='Total Cost')
plt.axvline(x=T_c_optimal, color='red', linestyle='--', linewidth=2, label=f'Optimal T_c = {T_c_optimal:.1f}°C')
plt.xlabel('Cooling Temperature (°C)')
plt.ylabel('Total Cost ($)')
plt.title('Sensitivity Analysis: Cost vs Temperature')
plt.grid(True, alpha=0.3)
plt.legend()

# 2. Cost vs Flow Rate (Temperature fixed at optimal)
ax2 = plt.subplot(2, 3, 2)
plt.plot(F_range, costs_F, 'g-', linewidth=2, label='Total Cost')
plt.axvline(x=F_optimal, color='red', linestyle='--', linewidth=2, label=f'Optimal F = {F_optimal:.1f} m³/min')
plt.xlabel('Flow Rate (m³/min)')
plt.ylabel('Total Cost ($)')
plt.title('Sensitivity Analysis: Cost vs Flow Rate')
plt.grid(True, alpha=0.3)
plt.legend()

# 3. 3D Cost Surface
ax3 = plt.subplot(2, 3, 3, projection='3d')
T_mesh, F_mesh, Cost_mesh = optimizer.create_cost_surface()
surface = ax3.plot_surface(T_mesh, F_mesh, Cost_mesh, cmap='viridis', alpha=0.8)
ax3.scatter([T_c_optimal], [F_optimal], [result.fun], color='red', s=100, label='Optimal Point')
ax3.set_xlabel('Temperature (°C)')
ax3.set_ylabel('Flow Rate (m³/min)')
ax3.set_zlabel('Total Cost ($)')
ax3.set_title('3D Cost Surface')
plt.colorbar(surface, ax=ax3, shrink=0.5)

# 4. Cost Components Analysis
ax4 = plt.subplot(2, 3, 4)
T_analysis = np.linspace(optimizer.T_min, optimizer.T_max, 50)
energy_costs = [optimizer.cooling_energy_cost(T) for T in T_analysis]
system_costs = [optimizer.cooling_system_cost(T, F_optimal) for T in T_analysis]

plt.plot(T_analysis, energy_costs, 'b-', linewidth=2, label='Energy Cost')
plt.plot(T_analysis, system_costs, 'r-', linewidth=2, label='Cooling System Cost')
plt.plot(T_analysis, np.array(energy_costs) + np.array(system_costs), 'k--', linewidth=2, label='Total Cost')
plt.axvline(x=T_c_optimal, color='orange', linestyle=':', linewidth=2, label='Optimal Point')
plt.xlabel('Cooling Temperature (°C)')
plt.ylabel('Cost ($)')
plt.title('Cost Components vs Temperature')
plt.legend()
plt.grid(True, alpha=0.3)

# 5. Contour plot
ax5 = plt.subplot(2, 3, 5)
contour = plt.contour(T_mesh, F_mesh, Cost_mesh, levels=20, colors='black', alpha=0.4)
plt.contourf(T_mesh, F_mesh, Cost_mesh, levels=20, cmap='RdYlBu_r', alpha=0.8)
plt.colorbar(label='Total Cost ($)')
plt.plot(T_c_optimal, F_optimal, 'r*', markersize=15, label='Optimal Point')
plt.xlabel('Temperature (°C)')
plt.ylabel('Flow Rate (m³/min)')
plt.title('Cost Contour Map')
plt.legend()

# 6. Server power factor visualization
ax6 = plt.subplot(2, 3, 6)
T_power = np.linspace(15, 35, 100)
power_factors = [optimizer.server_power_factor(T) for T in T_power]
plt.plot(T_power, power_factors, 'purple', linewidth=2)
plt.axvline(x=T_c_optimal, color='red', linestyle='--', linewidth=2, label=f'Optimal T_c = {T_c_optimal:.1f}°C')
plt.xlabel('Temperature (°C)')
plt.ylabel('Server Power Factor')
plt.title('Server Power Factor vs Temperature')
plt.grid(True, alpha=0.3)
plt.legend()

plt.tight_layout()
plt.show()

# Print detailed analysis
print("\n=== Detailed Analysis ===")
print(f"At optimal point:")
print(f"  Server power factor: {optimizer.server_power_factor(T_c_optimal):.3f}")
print(f"  Temperature deviation from ambient: {abs(T_c_optimal - optimizer.T_ambient):.1f}°C")
print(f"  Cost breakdown:")
print(f"    - Energy cost: ${energy_cost:.2f} ({energy_cost/(energy_cost+system_cost)*100:.1f}%)")
print(f"    - System cost: ${system_cost:.2f} ({system_cost/(energy_cost+system_cost)*100:.1f}%)")

# Compare with suboptimal solutions
print(f"\n=== Comparison with Suboptimal Solutions ===")
scenarios = [
    ("High Temperature", optimizer.T_max, F_optimal),
    ("Low Temperature", optimizer.T_min, F_optimal),
    ("High Flow Rate", T_c_optimal, optimizer.F_max),
    ("Low Flow Rate", T_c_optimal, optimizer.F_min)
]

for name, T, F in scenarios:
    cost = optimizer.total_cost([T, F])
    savings = ((cost - result.fun) / cost) * 100
    print(f"{name}: T={T:.1f}°C, F={F:.1f} m³/min, Cost=${cost:.2f}, Savings={savings:.1f}%")

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

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

class MaxEntropyTwoLevel:
    """
    Class to analyze maximum entropy principle for a two-level system
    """
    
    def __init__(self, epsilon=1.0, k_B=1.0):
        """
        Initialize the two-level system
        
        Parameters:
        epsilon (float): Energy difference between levels
        k_B (float): Boltzmann constant (set to 1 for simplicity)
        """
        self.epsilon = epsilon
        self.k_B = k_B
    
    def entropy(self, p1):
        """
        Calculate entropy for given probability p1
        S = -k_B * [p0*ln(p0) + p1*ln(p1)]
        
        Parameters:
        p1 (float): Probability of excited state
        
        Returns:
        float: Entropy value
        """
        p0 = 1 - p1
        
        # Handle edge cases to avoid log(0)
        if p0 <= 0 or p1 <= 0:
            return -np.inf
        
        return -self.k_B * (p0 * np.log(p0) + p1 * np.log(p1))
    
    def average_energy(self, p1):
        """
        Calculate average energy
        <E> = p0*E0 + p1*E1 = p1*epsilon
        
        Parameters:
        p1 (float): Probability of excited state
        
        Returns:
        float: Average energy
        """
        return p1 * self.epsilon
    
    def analytical_solution(self, avg_energy):
        """
        Analytical solution using Lagrange multipliers
        The result is the canonical distribution: p1 = exp(-beta*epsilon) / Z
        where beta is determined by the constraint <E> = avg_energy
        
        Parameters:
        avg_energy (float): Target average energy
        
        Returns:
        tuple: (p0, p1, beta, temperature)
        """
        if avg_energy <= 0:
            return 1.0, 0.0, np.inf, 0.0
        elif avg_energy >= self.epsilon:
            return 0.0, 1.0, -np.inf, np.inf
        
        # From constraint: p1 = avg_energy / epsilon
        p1 = avg_energy / self.epsilon
        p0 = 1 - p1
        
        # Calculate beta from the canonical distribution
        if p1 > 0 and p0 > 0:
            beta = np.log(p0/p1) / self.epsilon
            temperature = 1.0 / (self.k_B * beta) if beta > 0 else np.inf
        else:
            beta = 0
            temperature = np.inf
            
        return p0, p1, beta, temperature
    
    def numerical_optimization(self, avg_energy):
        """
        Numerical solution using scipy.optimize
        Maximize entropy subject to constraints
        
        Parameters:
        avg_energy (float): Target average energy
        
        Returns:
        tuple: (p0, p1, max_entropy, optimization_result)
        """
        
        # Objective function (negative entropy for minimization)
        def objective(p):
            return -self.entropy(p[0])  # p[0] is p1
        
        # Constraint: average energy
        def energy_constraint(p):
            return self.average_energy(p[0]) - avg_energy
        
        # Constraint: normalization (p0 + p1 = 1 is automatically satisfied)
        constraints = [{'type': 'eq', 'fun': energy_constraint}]
        
        # Bounds: 0 <= p1 <= 1
        bounds = [(0, 1)]
        
        # Initial guess
        x0 = [avg_energy / self.epsilon]
        
        # Optimize
        result = minimize(objective, x0, method='SLSQP', 
                         bounds=bounds, constraints=constraints)
        
        p1_opt = result.x[0]
        p0_opt = 1 - p1_opt
        max_entropy = self.entropy(p1_opt)
        
        return p0_opt, p1_opt, max_entropy, result
    
    def plot_entropy_surface(self):
        """
        Plot entropy as a function of p1 and average energy
        """
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
        
        # Plot 1: Entropy vs p1
        p1_values = np.linspace(0.001, 0.999, 1000)
        entropy_values = [self.entropy(p1) for p1 in p1_values]
        
        ax1.plot(p1_values, entropy_values, 'b-', linewidth=2, label='S(p₁)')
        ax1.set_xlabel('p₁ (Probability of Excited State)', fontsize=12)
        ax1.set_ylabel('Entropy S', fontsize=12)
        ax1.set_title('Entropy vs Probability Distribution', fontsize=14)
        ax1.grid(True, alpha=0.3)
        ax1.legend()
        
        # Maximum entropy point (unconstrained)
        max_entropy_p1 = 0.5
        max_entropy_val = self.entropy(max_entropy_p1)
        ax1.plot(max_entropy_p1, max_entropy_val, 'ro', markersize=8, 
                label=f'Maximum (p₁={max_entropy_p1}, S={max_entropy_val:.3f})')
        ax1.legend()
        
        # Plot 2: Entropy vs Average Energy (with constraint)
        avg_energies = np.linspace(0.01, self.epsilon*0.99, 100)
        constrained_entropies = []
        p1_constrained = []
        
        for avg_e in avg_energies:
            p0, p1, _, _ = self.analytical_solution(avg_e)
            constrained_entropies.append(self.entropy(p1))
            p1_constrained.append(p1)
        
        ax2.plot(avg_energies, constrained_entropies, 'g-', linewidth=2, 
                label='Constrained Maximum Entropy')
        ax2.set_xlabel('Average Energy ⟨E⟩', fontsize=12)
        ax2.set_ylabel('Maximum Entropy S', fontsize=12)
        ax2.set_title('Maximum Entropy vs Average Energy Constraint', fontsize=14)
        ax2.grid(True, alpha=0.3)
        ax2.legend()
        
        plt.tight_layout()
        plt.show()
        
        return fig
    
    def plot_temperature_analysis(self):
        """
        Plot temperature dependence and compare with analytical solution
        """
        fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
        
        # Temperature range
        temperatures = np.logspace(-1, 1.5, 100)  # 0.1 to ~30
        betas = 1.0 / (self.k_B * temperatures)
        
        # Calculate probabilities using canonical distribution
        Z = 1 + np.exp(-betas * self.epsilon)  # Partition function
        p0_canonical = 1 / Z
        p1_canonical = np.exp(-betas * self.epsilon) / Z
        
        # Average energies and entropies
        avg_energies_canonical = p1_canonical * self.epsilon
        entropies_canonical = -self.k_B * (p0_canonical * np.log(p0_canonical) + 
                                          p1_canonical * np.log(p1_canonical))
        
        # Plot 1: Probabilities vs Temperature
        ax1.semilogx(temperatures, p0_canonical, 'b-', linewidth=2, label='p₀ (Ground state)')
        ax1.semilogx(temperatures, p1_canonical, 'r-', linewidth=2, label='p₁ (Excited state)')
        ax1.set_xlabel('Temperature T', fontsize=12)
        ax1.set_ylabel('Probability', fontsize=12)
        ax1.set_title('State Probabilities vs Temperature', fontsize=14)
        ax1.legend()
        ax1.grid(True, alpha=0.3)
        
        # Plot 2: Average Energy vs Temperature
        ax2.semilogx(temperatures, avg_energies_canonical, 'g-', linewidth=2)
        ax2.set_xlabel('Temperature T', fontsize=12)
        ax2.set_ylabel('Average Energy ⟨E⟩', fontsize=12)
        ax2.set_title('Average Energy vs Temperature', fontsize=14)
        ax2.grid(True, alpha=0.3)
        ax2.axhline(y=self.epsilon/2, color='k', linestyle='--', alpha=0.5, 
                   label='ε/2 (high-T limit)')
        ax2.legend()
        
        # Plot 3: Entropy vs Temperature
        ax3.semilogx(temperatures, entropies_canonical, 'purple', linewidth=2)
        ax3.set_xlabel('Temperature T', fontsize=12)
        ax3.set_ylabel('Entropy S', fontsize=12)
        ax3.set_title('Entropy vs Temperature', fontsize=14)
        ax3.grid(True, alpha=0.3)
        ax3.axhline(y=self.k_B*np.log(2), color='k', linestyle='--', alpha=0.5, 
                   label=f'k_B ln(2) = {self.k_B*np.log(2):.3f}')
        ax3.legend()
        
        # Plot 4: Phase diagram (Temperature vs Average Energy)
        ax4.plot(avg_energies_canonical, temperatures, 'orange', linewidth=2)
        ax4.set_xlabel('Average Energy ⟨E⟩', fontsize=12)
        ax4.set_ylabel('Temperature T', fontsize=12)
        ax4.set_title('Phase Diagram: Temperature vs Energy', fontsize=14)
        ax4.grid(True, alpha=0.3)
        ax4.set_yscale('log')
        
        plt.tight_layout()
        plt.show()
        
        return fig
    
    def compare_solutions(self, test_energies):
        """
        Compare analytical and numerical solutions
        
        Parameters:
        test_energies (array): Array of average energies to test
        
        Returns:
        dict: Comparison results
        """
        results = {
            'avg_energies': test_energies,
            'analytical': {'p0': [], 'p1': [], 'entropy': [], 'temperature': []},
            'numerical': {'p0': [], 'p1': [], 'entropy': []},
            'differences': {'p0': [], 'p1': [], 'entropy': []}
        }
        
        print("Comparison of Analytical vs Numerical Solutions")
        print("=" * 60)
        print(f"{'Avg Energy':>10} {'p₀ (Ana)':>10} {'p₀ (Num)':>10} {'p₁ (Ana)':>10} {'p₁ (Num)':>10} {'Temp':>10}")
        print("-" * 60)
        
        for avg_e in test_energies:
            # Analytical solution
            p0_ana, p1_ana, beta, temp = self.analytical_solution(avg_e)
            entropy_ana = self.entropy(p1_ana)
            
            # Numerical solution
            p0_num, p1_num, entropy_num, _ = self.numerical_optimization(avg_e)
            
            # Store results
            results['analytical']['p0'].append(p0_ana)
            results['analytical']['p1'].append(p1_ana)
            results['analytical']['entropy'].append(entropy_ana)
            results['analytical']['temperature'].append(temp)
            
            results['numerical']['p0'].append(p0_num)
            results['numerical']['p1'].append(p1_num)
            results['numerical']['entropy'].append(entropy_num)
            
            results['differences']['p0'].append(abs(p0_ana - p0_num))
            results['differences']['p1'].append(abs(p1_ana - p1_num))
            results['differences']['entropy'].append(abs(entropy_ana - entropy_num))
            
            print(f"{avg_e:10.3f} {p0_ana:10.4f} {p0_num:10.4f} {p1_ana:10.4f} {p1_num:10.4f} {temp:10.2f}")
        
        return results

# Example usage and analysis
def main():
    # Create system with epsilon = 2.0 (energy units)
    system = MaxEntropyTwoLevel(epsilon=2.0, k_B=1.0)
    
    print("Maximum Entropy Principle: Two-Level System Analysis")
    print("=" * 55)
    print(f"System parameters: ε = {system.epsilon}, k_B = {system.k_B}")
    print()
    
    # Test case: average energy = 0.8
    avg_energy_test = 0.8
    
    print(f"Case Study: Average Energy = {avg_energy_test}")
    print("-" * 30)
    
    # Analytical solution
    p0_ana, p1_ana, beta, temp = system.analytical_solution(avg_energy_test)
    entropy_ana = system.entropy(p1_ana)
    
    print("Analytical Solution (Lagrange multipliers):")
    print(f"  p₀ = {p0_ana:.4f}")
    print(f"  p₁ = {p1_ana:.4f}")
    print(f"  β = {beta:.4f}")
    print(f"  Temperature = {temp:.4f}")
    print(f"  Entropy = {entropy_ana:.4f}")
    print()
    
    # Numerical solution
    p0_num, p1_num, entropy_num, result = system.numerical_optimization(avg_energy_test)
    
    print("Numerical Solution (scipy.optimize):")
    print(f"  p₀ = {p0_num:.4f}")
    print(f"  p₁ = {p1_num:.4f}")
    print(f"  Entropy = {entropy_num:.4f}")
    print(f"  Optimization success: {result.success}")
    print()
    
    # Verification
    print("Verification:")
    print(f"  Normalization: p₀ + p₁ = {p0_ana + p1_ana:.6f}")
    print(f"  Average energy: ⟨E⟩ = {system.average_energy(p1_ana):.6f}")
    print(f"  Target energy: {avg_energy_test}")
    print()
    
    # Generate plots
    print("Generating visualization plots...")
    system.plot_entropy_surface()
    system.plot_temperature_analysis()
    
    # Compare solutions for different energies
    test_energies = np.linspace(0.1, 1.9, 10)
    comparison = system.compare_solutions(test_energies)
    
    print()
    print("Maximum differences between analytical and numerical solutions:")
    print(f"  Δp₀_max = {max(comparison['differences']['p0']):.2e}")
    print(f"  Δp₁_max = {max(comparison['differences']['p1']):.2e}")
    print(f"  ΔS_max = {max(comparison['differences']['entropy']):.2e}")

if __name__ == "__main__":
    main()