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

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

class TradePortfolioOptimizer:
    """
    A class to optimize trade partner portfolio allocation using Modern Portfolio Theory
    """
    
    def __init__(self, countries, expected_returns, covariance_matrix):
        """
        Initialize the optimizer with trade partner data
        
        Parameters:
        countries (list): List of trade partner country names
        expected_returns (np.array): Expected returns from each trade partner
        covariance_matrix (np.array): Covariance matrix of returns
        """
        self.countries = countries
        self.expected_returns = np.array(expected_returns)
        self.cov_matrix = np.array(covariance_matrix)
        self.n_assets = len(countries)
        
    def portfolio_performance(self, weights):
        """
        Calculate portfolio return and risk
        
        Parameters:
        weights (np.array): Portfolio weights
        
        Returns:
        tuple: (portfolio_return, portfolio_risk)
        """
        portfolio_return = np.sum(weights * self.expected_returns)
        portfolio_risk = np.sqrt(np.dot(weights.T, np.dot(self.cov_matrix, weights)))
        return portfolio_return, portfolio_risk
    
    def minimize_risk(self, target_return):
        """
        Minimize portfolio risk for a given target return
        
        Parameters:
        target_return (float): Desired portfolio return
        
        Returns:
        dict: Optimization result
        """
        # Objective function: minimize portfolio variance
        def objective(weights):
            return np.dot(weights.T, np.dot(self.cov_matrix, weights))
        
        # Constraints
        constraints = [
            {'type': 'eq', 'fun': lambda x: np.sum(x) - 1},  # weights sum to 1
            {'type': 'eq', 'fun': lambda x: np.sum(x * self.expected_returns) - target_return}  # target return
        ]
        
        # Bounds: each weight between 0 and 1
        bounds = tuple((0, 1) for _ in range(self.n_assets))
        
        # Initial guess: equal weights
        initial_guess = np.array([1/self.n_assets] * self.n_assets)
        
        # Optimize
        result = minimize(objective, initial_guess, method='SLSQP', 
                         bounds=bounds, constraints=constraints)
        
        return result
    
    def efficient_frontier(self, num_points=100):
        """
        Generate the efficient frontier
        
        Parameters:
        num_points (int): Number of points on the frontier
        
        Returns:
        tuple: (returns, risks, weights_array)
        """
        # Range of target returns
        min_return = np.min(self.expected_returns)
        max_return = np.max(self.expected_returns)
        target_returns = np.linspace(min_return, max_return, num_points)
        
        returns = []
        risks = []
        weights_array = []
        
        for target_return in target_returns:
            try:
                result = self.minimize_risk(target_return)
                if result.success:
                    weights = result.x
                    ret, risk = self.portfolio_performance(weights)
                    returns.append(ret)
                    risks.append(risk)
                    weights_array.append(weights)
            except:
                continue
                
        return np.array(returns), np.array(risks), np.array(weights_array)
    
    def max_sharpe_ratio(self, risk_free_rate=0.02):
        """
        Find portfolio with maximum Sharpe ratio
        
        Parameters:
        risk_free_rate (float): Risk-free rate
        
        Returns:
        dict: Optimization result
        """
        def negative_sharpe(weights):
            ret, risk = self.portfolio_performance(weights)
            return -(ret - risk_free_rate) / risk
        
        constraints = [{'type': 'eq', 'fun': lambda x: np.sum(x) - 1}]
        bounds = tuple((0, 1) for _ in range(self.n_assets))
        initial_guess = np.array([1/self.n_assets] * self.n_assets)
        
        result = minimize(negative_sharpe, initial_guess, method='SLSQP',
                         bounds=bounds, constraints=constraints)
        
        return result

# Example: Trade portfolio for a hypothetical country
def main():
    # Define trade partners and their characteristics
    countries = ['USA', 'China', 'Germany', 'Japan', 'UK', 'India', 'Brazil']
    
    # Expected returns (annual trade growth benefits as percentage)
    expected_returns = [0.08, 0.12, 0.06, 0.05, 0.07, 0.15, 0.11]
    
    # Covariance matrix (simplified for demonstration)
    # In reality, this would be calculated from historical trade data
    correlation_matrix = np.array([
        [1.00, 0.65, 0.70, 0.60, 0.75, 0.30, 0.40],  # USA
        [0.65, 1.00, 0.55, 0.50, 0.45, 0.60, 0.35],  # China
        [0.70, 0.55, 1.00, 0.65, 0.80, 0.25, 0.30],  # Germany
        [0.60, 0.50, 0.65, 1.00, 0.55, 0.35, 0.25],  # Japan
        [0.75, 0.45, 0.80, 0.55, 1.00, 0.20, 0.35],  # UK
        [0.30, 0.60, 0.25, 0.35, 0.20, 1.00, 0.45],  # India
        [0.40, 0.35, 0.30, 0.25, 0.35, 0.45, 1.00]   # Brazil
    ])
    
    # Standard deviations (trade volatility)
    std_devs = [0.15, 0.25, 0.12, 0.10, 0.14, 0.30, 0.22]
    
    # Calculate covariance matrix
    covariance_matrix = np.outer(std_devs, std_devs) * correlation_matrix
    
    # Create optimizer
    optimizer = TradePortfolioOptimizer(countries, expected_returns, covariance_matrix)
    
    # Generate efficient frontier
    print("Generating efficient frontier...")
    frontier_returns, frontier_risks, frontier_weights = optimizer.efficient_frontier()
    
    # Find maximum Sharpe ratio portfolio
    print("Finding maximum Sharpe ratio portfolio...")
    max_sharpe_result = optimizer.max_sharpe_ratio()
    max_sharpe_weights = max_sharpe_result.x
    max_sharpe_return, max_sharpe_risk = optimizer.portfolio_performance(max_sharpe_weights)
    
    # Find minimum variance portfolio
    print("Finding minimum variance portfolio...")
    min_var_result = optimizer.minimize_risk(np.min(frontier_returns))
    min_var_weights = min_var_result.x
    min_var_return, min_var_risk = optimizer.portfolio_performance(min_var_weights)
    
    # Create visualizations
    fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
    
    # 1. Efficient Frontier
    ax1.plot(frontier_risks, frontier_returns, 'b-', linewidth=2, label='Efficient Frontier')
    ax1.scatter(max_sharpe_risk, max_sharpe_return, marker='*', color='red', s=200, label='Max Sharpe Ratio')
    ax1.scatter(min_var_risk, min_var_return, marker='*', color='green', s=200, label='Min Variance')
    
    # Plot individual assets
    individual_risks = [np.sqrt(covariance_matrix[i,i]) for i in range(len(countries))]
    ax1.scatter(individual_risks, expected_returns, alpha=0.7, s=100, label='Individual Countries')
    
    for i, country in enumerate(countries):
        ax1.annotate(country, (individual_risks[i], expected_returns[i]), 
                    xytext=(5, 5), textcoords='offset points', fontsize=8)
    
    ax1.set_xlabel('Risk (Standard Deviation)')
    ax1.set_ylabel('Expected Return')
    ax1.set_title('Efficient Frontier - Trade Partner Portfolio')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. Maximum Sharpe Ratio Portfolio Allocation
    ax2.pie(max_sharpe_weights, labels=countries, autopct='%1.1f%%', startangle=90)
    ax2.set_title('Maximum Sharpe Ratio Portfolio\nAllocation')
    
    # 3. Minimum Variance Portfolio Allocation
    ax3.pie(min_var_weights, labels=countries, autopct='%1.1f%%', startangle=90)
    ax3.set_title('Minimum Variance Portfolio\nAllocation')
    
    # 4. Risk-Return Comparison
    portfolios = ['Max Sharpe', 'Min Variance']
    returns_comp = [max_sharpe_return, min_var_return]
    risks_comp = [max_sharpe_risk, min_var_risk]
    
    x = np.arange(len(portfolios))
    width = 0.35
    
    ax4_twin = ax4.twinx()
    bars1 = ax4.bar(x - width/2, returns_comp, width, label='Expected Return', alpha=0.8, color='blue')
    bars2 = ax4_twin.bar(x + width/2, risks_comp, width, label='Risk', alpha=0.8, color='red')
    
    ax4.set_xlabel('Portfolio Type')
    ax4.set_ylabel('Expected Return', color='blue')
    ax4_twin.set_ylabel('Risk (Std Dev)', color='red')
    ax4.set_title('Portfolio Comparison: Risk vs Return')
    ax4.set_xticks(x)
    ax4.set_xticklabels(portfolios)
    ax4.grid(True, alpha=0.3)
    
    # Add value labels on bars
    for bar in bars1:
        height = bar.get_height()
        ax4.annotate(f'{height:.3f}', xy=(bar.get_x() + bar.get_width()/2, height),
                    xytext=(0, 3), textcoords="offset points", ha='center', va='bottom')
    
    for bar in bars2:
        height = bar.get_height()
        ax4_twin.annotate(f'{height:.3f}', xy=(bar.get_x() + bar.get_width()/2, height),
                         xytext=(0, 3), textcoords="offset points", ha='center', va='bottom')
    
    plt.tight_layout()
    plt.show()
    
    # Print detailed results
    print("\n" + "="*60)
    print("TRADE PARTNER PORTFOLIO OPTIMIZATION RESULTS")
    print("="*60)
    
    print(f"\nMAXIMUM SHARPE RATIO PORTFOLIO:")
    print(f"Expected Return: {max_sharpe_return:.4f} ({max_sharpe_return*100:.2f}%)")
    print(f"Risk (Std Dev): {max_sharpe_risk:.4f} ({max_sharpe_risk*100:.2f}%)")
    print(f"Sharpe Ratio: {(max_sharpe_return - 0.02) / max_sharpe_risk:.4f}")
    print("\nAllocation:")
    for i, country in enumerate(countries):
        if max_sharpe_weights[i] > 0.01:  # Only show significant allocations
            print(f"  {country}: {max_sharpe_weights[i]*100:.2f}%")
    
    print(f"\nMINIMUM VARIANCE PORTFOLIO:")
    print(f"Expected Return: {min_var_return:.4f} ({min_var_return*100:.2f}%)")
    print(f"Risk (Std Dev): {min_var_risk:.4f} ({min_var_risk*100:.2f}%)")
    print("\nAllocation:")
    for i, country in enumerate(countries):
        if min_var_weights[i] > 0.01:  # Only show significant allocations
            print(f"  {country}: {min_var_weights[i]*100:.2f}%")
    
    # Risk-Return Analysis Table
    print(f"\nINDIVIDUAL COUNTRY ANALYSIS:")
    print("-" * 50)
    df_analysis = pd.DataFrame({
        'Country': countries,
        'Expected Return': [f"{r*100:.2f}%" for r in expected_returns],
        'Risk (Std Dev)': [f"{np.sqrt(covariance_matrix[i,i])*100:.2f}%" for i in range(len(countries))],
        'Return/Risk Ratio': [expected_returns[i]/np.sqrt(covariance_matrix[i,i]) for i in range(len(countries))]
    })
    print(df_analysis.to_string(index=False))
    
    return optimizer, max_sharpe_weights, min_var_weights

# Run the analysis
if __name__ == "__main__":
    optimizer, max_sharpe_weights, min_var_weights = main()

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize, differential_evolution
from scipy.constants import hbar, e, epsilon_0, pi
import warnings
warnings.filterwarnings('ignore')

# Physical constants and material parameters for CdSe
class QuantumDotParameters:
    def __init__(self):
        # CdSe material parameters
        self.E_g_bulk = 1.74  # eV, bulk bandgap
        self.m_e_star = 0.13  # effective electron mass (in units of m_0)
        self.m_h_star = 0.45  # effective hole mass (in units of m_0)
        self.epsilon_r = 10.0  # relative dielectric constant
        self.m_0 = 9.109e-31  # electron rest mass (kg)
        
        # Convert to convenient units
        self.hbar_eV = hbar / e  # hbar in eV·s
        
    def effective_bandgap(self, radius_nm):
        """
        Calculate effective bandgap using quantum confinement model
        
        Parameters:
        radius_nm: quantum dot radius in nanometers
        
        Returns:
        Effective bandgap in eV
        """
        R = radius_nm * 1e-9  # convert to meters
        
        # Confinement energy term
        confinement = (self.hbar_eV**2 * pi**2) / (2 * R**2) * \
                     (1/(self.m_e_star * self.m_0) + 1/(self.m_h_star * self.m_0)) * \
                     (1/(e * 1e-18))  # Convert to eV
        
        # Coulomb interaction term (binding energy)
        coulomb = (1.8 * e**2) / (4 * pi * self.epsilon_r * epsilon_0 * R) / e
        
        return self.E_g_bulk + confinement - coulomb
    
    def emission_wavelength(self, radius_nm):
        """Convert bandgap to emission wavelength in nm"""
        E_gap = self.effective_bandgap(radius_nm)
        # E = hc/λ, so λ = hc/E
        wavelength_nm = (hbar * 3e8 * 2 * pi) / (E_gap * e) * 1e9
        return wavelength_nm
    
    def quantum_yield(self, radius_nm, surface_quality=0.8):
        """
        Model quantum yield based on size and surface quality
        
        Parameters:
        radius_nm: quantum dot radius
        surface_quality: factor (0-1) representing surface passivation quality
        """
        # Empirical model: QY decreases with smaller size due to surface effects
        # and has optimal range around 2-4 nm for CdSe
        size_factor = np.exp(-(radius_nm - 3.0)**2 / 2.0)  # Gaussian centered at 3nm
        surface_factor = surface_quality
        
        # Base quantum yield with size and surface effects
        qy = 0.95 * size_factor * surface_factor
        
        # Add some realistic constraints
        if radius_nm < 1.0:  # Very small dots have poor QY
            qy *= 0.1
        elif radius_nm > 6.0:  # Very large dots approach bulk behavior
            qy *= 0.7
            
        return min(qy, 1.0)

# Initialize parameters
qd_params = QuantumDotParameters()

# Define optimization objectives
class QuantumDotOptimizer:
    def __init__(self, params, target_wavelength=550, weight_qy=0.5, weight_wavelength=0.5):
        self.params = params
        self.target_wavelength = target_wavelength  # target emission wavelength in nm
        self.weight_qy = weight_qy
        self.weight_wavelength = weight_wavelength
        
    def objective_function(self, x):
        """
        Multi-objective optimization function
        x[0]: radius in nm
        x[1]: surface quality factor (0-1)
        """
        radius, surface_quality = x
        
        # Calculate properties
        wavelength = self.params.emission_wavelength(radius)
        qy = self.params.quantum_yield(radius, surface_quality)
        
        # Wavelength deviation penalty
        wavelength_penalty = abs(wavelength - self.target_wavelength) / self.target_wavelength
        
        # Objective: minimize wavelength deviation, maximize QY
        objective = self.weight_wavelength * wavelength_penalty - self.weight_qy * qy
        
        return objective
    
    def optimize(self, method='differential_evolution'):
        """Perform optimization"""
        bounds = [(1.0, 6.0),    # radius bounds in nm
                 (0.5, 1.0)]     # surface quality bounds
        
        if method == 'differential_evolution':
            result = differential_evolution(
                self.objective_function, 
                bounds, 
                maxiter=1000,
                popsize=15,
                seed=42
            )
        else:
            # Alternative: local optimization
            x0 = [3.0, 0.8]  # initial guess
            result = minimize(
                self.objective_function, 
                x0, 
                bounds=bounds,
                method='L-BFGS-B'
            )
            
        return result

# Perform optimization for different target wavelengths
target_wavelengths = [450, 500, 550, 600, 650]  # Blue to Red
optimization_results = {}

print("Quantum Dot Optimization Results")
print("="*50)

for target_wl in target_wavelengths:
    optimizer = QuantumDotOptimizer(qd_params, target_wavelength=target_wl)
    result = optimizer.optimize()
    
    optimal_radius = result.x[0]
    optimal_surface_quality = result.x[1]
    
    # Calculate final properties
    final_wavelength = qd_params.emission_wavelength(optimal_radius)
    final_qy = qd_params.quantum_yield(optimal_radius, optimal_surface_quality)
    final_bandgap = qd_params.effective_bandgap(optimal_radius)
    
    optimization_results[target_wl] = {
        'radius': optimal_radius,
        'surface_quality': optimal_surface_quality,
        'wavelength': final_wavelength,
        'quantum_yield': final_qy,
        'bandgap': final_bandgap,
        'success': result.success
    }
    
    print(f"\nTarget wavelength: {target_wl} nm")
    print(f"Optimal radius: {optimal_radius:.2f} nm")
    print(f"Optimal surface quality: {optimal_surface_quality:.3f}")
    print(f"Achieved wavelength: {final_wavelength:.1f} nm")
    print(f"Quantum yield: {final_qy:.3f}")
    print(f"Effective bandgap: {final_bandgap:.3f} eV")
    print(f"Optimization success: {result.success}")

# Create comprehensive visualizations
plt.style.use('seaborn-v0_8')
fig = plt.figure(figsize=(16, 12))

# Plot 1: Size-dependent properties
ax1 = plt.subplot(2, 3, 1)
radius_range = np.linspace(1, 6, 100)
wavelengths = [qd_params.emission_wavelength(r) for r in radius_range]
bandgaps = [qd_params.effective_bandgap(r) for r in radius_range]

ax1.plot(radius_range, wavelengths, 'b-', linewidth=2, label='Emission wavelength')
ax1.set_xlabel('Radius (nm)')
ax1.set_ylabel('Wavelength (nm)', color='b')
ax1.tick_params(axis='y', labelcolor='b')
ax1.grid(True, alpha=0.3)

ax1_twin = ax1.twinx()
ax1_twin.plot(radius_range, bandgaps, 'r--', linewidth=2, label='Bandgap')
ax1_twin.set_ylabel('Bandgap (eV)', color='r')
ax1_twin.tick_params(axis='y', labelcolor='r')

ax1.set_title('Size-Dependent Optical Properties')

# Plot 2: Quantum yield vs size and surface quality
ax2 = plt.subplot(2, 3, 2)
surface_qualities = [0.6, 0.8, 1.0]
for sq in surface_qualities:
    qy_values = [qd_params.quantum_yield(r, sq) for r in radius_range]
    ax2.plot(radius_range, qy_values, linewidth=2, label=f'Surface quality: {sq}')

ax2.set_xlabel('Radius (nm)')
ax2.set_ylabel('Quantum Yield')
ax2.set_title('Quantum Yield vs Size and Surface Quality')
ax2.legend()
ax2.grid(True, alpha=0.3)

# Plot 3: Optimization results
ax3 = plt.subplot(2, 3, 3)
target_wls = list(optimization_results.keys())
optimal_radii = [optimization_results[wl]['radius'] for wl in target_wls]
achieved_wls = [optimization_results[wl]['wavelength'] for wl in target_wls]

ax3.scatter(target_wls, achieved_wls, c=optimal_radii, s=100, cmap='viridis', alpha=0.8)
ax3.plot([400, 700], [400, 700], 'k--', alpha=0.5, label='Perfect match')
ax3.set_xlabel('Target Wavelength (nm)')
ax3.set_ylabel('Achieved Wavelength (nm)')
ax3.set_title('Optimization Accuracy')
ax3.legend()
ax3.grid(True, alpha=0.3)

# Add colorbar for radius
cbar = plt.colorbar(ax3.collections[0], ax=ax3)
cbar.set_label('Optimal Radius (nm)')

# Plot 4: Quantum efficiency map
ax4 = plt.subplot(2, 3, 4)
radius_grid = np.linspace(1, 6, 50)
surface_grid = np.linspace(0.5, 1.0, 50)
R, S = np.meshgrid(radius_grid, surface_grid)
QY = np.zeros_like(R)

for i in range(len(radius_grid)):
    for j in range(len(surface_grid)):
        QY[j, i] = qd_params.quantum_yield(R[j, i], S[j, i])

contour = ax4.contourf(R, S, QY, levels=20, cmap='plasma')
ax4.set_xlabel('Radius (nm)')
ax4.set_ylabel('Surface Quality')
ax4.set_title('Quantum Yield Optimization Map')
plt.colorbar(contour, ax=ax4, label='Quantum Yield')

# Plot 5: Spectral tuning range
ax5 = plt.subplot(2, 3, 5)
colors = plt.cm.rainbow(np.linspace(0, 1, len(target_wls)))
for i, (target_wl, color) in enumerate(zip(target_wls, colors)):
    result = optimization_results[target_wl]
    ax5.bar(i, result['quantum_yield'], color=color, alpha=0.7, 
            label=f'{target_wl}nm (R={result["radius"]:.1f}nm)')

ax5.set_xlabel('Target Wavelength')
ax5.set_ylabel('Quantum Yield')
ax5.set_title('Achievable Performance Across Spectrum')
ax5.set_xticks(range(len(target_wls)))
ax5.set_xticklabels([f'{wl}nm' for wl in target_wls])
ax5.grid(True, alpha=0.3)

# Plot 6: Bandgap engineering
ax6 = plt.subplot(2, 3, 6)
radii_for_bg = [optimization_results[wl]['radius'] for wl in target_wls]
bandgaps_optimized = [optimization_results[wl]['bandgap'] for wl in target_wls]

ax6.scatter(radii_for_bg, bandgaps_optimized, c=target_wls, s=100, cmap='rainbow', alpha=0.8)
ax6.set_xlabel('Optimal Radius (nm)')
ax6.set_ylabel('Effective Bandgap (eV)')
ax6.set_title('Bandgap Engineering Results')
ax6.grid(True, alpha=0.3)

# Add colorbar for wavelength
cbar2 = plt.colorbar(ax6.collections[0], ax=ax6)
cbar2.set_label('Target Wavelength (nm)')

plt.tight_layout()
plt.show()

# Summary analysis
print("\n" + "="*50)
print("OPTIMIZATION SUMMARY AND ANALYSIS")
print("="*50)

print(f"\nSpectral Range: {min(achieved_wls):.1f} - {max(achieved_wls):.1f} nm")
print(f"Size Range: {min(optimal_radii):.2f} - {max(optimal_radii):.2f} nm")

avg_qy = np.mean([result['quantum_yield'] for result in optimization_results.values()])
print(f"Average Quantum Yield: {avg_qy:.3f}")

wavelength_accuracy = np.mean([abs(result['wavelength'] - target) / target * 100 
                              for target, result in optimization_results.items()])
print(f"Average Wavelength Accuracy: {wavelength_accuracy:.1f}% deviation")

print("\nKey Insights:")
print("1. Smaller quantum dots (1-2 nm) emit in blue region but have lower QY")
print("2. Optimal size range for high QY is 2.5-4 nm")
print("3. Surface quality is crucial for maintaining high quantum efficiency")
print("4. Bandgap can be tuned from ~3.5 eV (blue) to ~2.2 eV (red)")
print("5. Trade-off exists between wavelength precision and quantum yield optimization")

objective = self.weight_wavelength * wavelength_penalty - self.weight_qy * qy

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import seaborn as sns
from scipy.optimize import dual_annealing
import time

class IsingModel2D:
    """
    2D Ising Model implementation with simulated annealing optimization
    """
    
    def __init__(self, size=20, J=1.0, h=0.0):
        """
        Initialize the 2D Ising model
        
        Parameters:
        - size: Grid size (size x size)
        - J: Coupling constant (positive = ferromagnetic, negative = antiferromagnetic)
        - h: External magnetic field
        """
        self.size = size
        self.J = J
        self.h = h
        self.spins = np.random.choice([-1, 1], size=(size, size))
        self.energy_history = []
        self.magnetization_history = []
        
    def calculate_energy(self, spins=None):
        """
        Calculate the total energy of the spin configuration
        
        E = -J * sum(s_i * s_j) - h * sum(s_i)
        """
        if spins is None:
            spins = self.spins
            
        # Interaction energy (nearest neighbors)
        # Using periodic boundary conditions
        interaction_energy = 0
        
        # Horizontal interactions
        horizontal = np.sum(spins[:, :-1] * spins[:, 1:])  # Adjacent columns
        horizontal += np.sum(spins[:, -1] * spins[:, 0])   # Periodic boundary
        
        # Vertical interactions
        vertical = np.sum(spins[:-1, :] * spins[1:, :])    # Adjacent rows
        vertical += np.sum(spins[-1, :] * spins[0, :])     # Periodic boundary
        
        interaction_energy = -(self.J * (horizontal + vertical))
        
        # External field energy
        field_energy = -self.h * np.sum(spins)
        
        return interaction_energy + field_energy
    
    def calculate_magnetization(self, spins=None):
        """
        Calculate the magnetization (average spin)
        """
        if spins is None:
            spins = self.spins
        return np.mean(spins)
    
    def metropolis_step(self, temperature):
        """
        Perform one Monte Carlo step using Metropolis algorithm
        """
        # Choose a random spin
        i, j = np.random.randint(0, self.size, 2)
        
        # Calculate energy change if we flip this spin
        current_spin = self.spins[i, j]
        
        # Get neighbors with periodic boundary conditions
        neighbors = [
            self.spins[(i-1) % self.size, j],
            self.spins[(i+1) % self.size, j],
            self.spins[i, (j-1) % self.size],
            self.spins[i, (j+1) % self.size]
        ]
        
        # Energy change from flipping the spin
        delta_E = 2 * current_spin * (self.J * sum(neighbors) + self.h)
        
        # Accept or reject the flip
        if delta_E <= 0 or np.random.random() < np.exp(-delta_E / temperature):
            self.spins[i, j] = -current_spin
            return True
        return False
    
    def simulated_annealing(self, initial_temp=10.0, final_temp=0.1, 
                          cooling_rate=0.95, steps_per_temp=1000, verbose=True):
        """
        Optimize spin configuration using simulated annealing
        """
        temperature = initial_temp
        accepted_moves = 0
        total_moves = 0
        
        self.energy_history = []
        self.magnetization_history = []
        self.temperature_history = []
        
        print("Starting simulated annealing optimization...")
        print(f"Initial energy: {self.calculate_energy():.2f}")
        print(f"Initial magnetization: {self.calculate_magnetization():.3f}")
        print("-" * 50)
        
        while temperature > final_temp:
            for _ in range(steps_per_temp):
                if self.metropolis_step(temperature):
                    accepted_moves += 1
                total_moves += 1
            
            # Record current state
            current_energy = self.calculate_energy()
            current_magnetization = self.calculate_magnetization()
            
            self.energy_history.append(current_energy)
            self.magnetization_history.append(current_magnetization)
            self.temperature_history.append(temperature)
            
            if verbose and len(self.energy_history) % 10 == 0:
                acceptance_rate = accepted_moves / total_moves if total_moves > 0 else 0
                print(f"T={temperature:.3f}, E={current_energy:.2f}, "
                      f"M={current_magnetization:.3f}, Accept={acceptance_rate:.2f}")
            
            # Cool down
            temperature *= cooling_rate
        
        final_energy = self.calculate_energy()
        final_magnetization = self.calculate_magnetization()
        
        print("-" * 50)
        print(f"Optimization completed!")
        print(f"Final energy: {final_energy:.2f}")
        print(f"Final magnetization: {final_magnetization:.3f}")
        print(f"Total acceptance rate: {accepted_moves/total_moves:.3f}")
        
        return {
            'final_energy': final_energy,
            'final_magnetization': final_magnetization,
            'acceptance_rate': accepted_moves/total_moves,
            'total_steps': total_moves
        }

def compare_different_parameters():
    """
    Compare optimization results for different J and h values
    """
    # Different parameter combinations
    params = [
        {'J': 1.0, 'h': 0.0, 'name': 'Ferromagnetic (J>0, h=0)'},
        {'J': -1.0, 'h': 0.0, 'name': 'Antiferromagnetic (J<0, h=0)'},
        {'J': 1.0, 'h': 0.5, 'name': 'Ferromagnetic with field (J>0, h>0)'},
        {'J': 1.0, 'h': -0.5, 'name': 'Ferromagnetic with opposite field (J>0, h<0)'}
    ]
    
    results = []
    models = []
    
    for param in params:
        print(f"\n{'='*60}")
        print(f"Testing: {param['name']}")
        print(f"{'='*60}")
        
        model = IsingModel2D(size=15, J=param['J'], h=param['h'])
        result = model.simulated_annealing(
            initial_temp=5.0, 
            final_temp=0.01, 
            cooling_rate=0.9, 
            steps_per_temp=500,
            verbose=False
        )
        
        results.append({**result, **param})
        models.append(model)
    
    return results, models

def plot_results(results, models):
    """
    Create comprehensive plots of the optimization results
    """
    fig = plt.figure(figsize=(20, 16))
    
    # 1. Final spin configurations
    for i, (model, result) in enumerate(zip(models, results)):
        ax = plt.subplot(4, 4, i+1)
        
        # Custom colormap for spins: red for -1, blue for +1
        colors = ['red', 'blue']
        cmap = ListedColormap(colors)
        
        im = ax.imshow(model.spins, cmap=cmap, vmin=-1, vmax=1)
        ax.set_title(f"{result['name']}\nE={result['final_energy']:.1f}, M={result['final_magnetization']:.2f}", 
                    fontsize=10)
        ax.set_xticks([])
        ax.set_yticks([])
        
        # Add colorbar for the first plot
        if i == 0:
            cbar = plt.colorbar(im, ax=ax, shrink=0.8)
            cbar.set_ticks([-1, 1])
            cbar.set_ticklabels(['↓ (-1)', '↑ (+1)'])
    
    # 2. Energy evolution during optimization
    ax_energy = plt.subplot(4, 4, (5, 8))
    for i, (model, result) in enumerate(zip(models, results)):
        if len(model.energy_history) > 0:
            ax_energy.plot(model.energy_history, label=result['name'], linewidth=2)
    
    ax_energy.set_xlabel('Cooling Steps')
    ax_energy.set_ylabel('Energy')
    ax_energy.set_title('Energy Evolution During Annealing')
    ax_energy.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    ax_energy.grid(True, alpha=0.3)
    
    # 3. Magnetization evolution
    ax_mag = plt.subplot(4, 4, (9, 12))
    for i, (model, result) in enumerate(zip(models, results)):
        if len(model.magnetization_history) > 0:
            ax_mag.plot(model.magnetization_history, label=result['name'], linewidth=2)
    
    ax_mag.set_xlabel('Cooling Steps')
    ax_mag.set_ylabel('Magnetization')
    ax_mag.set_title('Magnetization Evolution During Annealing')
    ax_mag.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    ax_mag.grid(True, alpha=0.3)
    ax_mag.axhline(y=0, color='black', linestyle='--', alpha=0.5)
    
    # 4. Temperature schedule
    ax_temp = plt.subplot(4, 4, (13, 16))
    if len(models[0].temperature_history) > 0:
        ax_temp.semilogy(models[0].temperature_history, 'g-', linewidth=3, label='Temperature')
        ax_temp.set_xlabel('Cooling Steps')
        ax_temp.set_ylabel('Temperature (log scale)')
        ax_temp.set_title('Temperature Schedule')
        ax_temp.grid(True, alpha=0.3)
        ax_temp.legend()
    
    plt.tight_layout()
    plt.show()
    
    # Summary statistics
    print("\n" + "="*80)
    print("OPTIMIZATION SUMMARY")
    print("="*80)
    print(f"{'Configuration':<35} {'Final Energy':<15} {'Magnetization':<15} {'Accept Rate':<12}")
    print("-"*80)
    for result in results:
        print(f"{result['name']:<35} {result['final_energy']:<15.2f} "
              f"{result['final_magnetization']:<15.3f} {result['acceptance_rate']:<12.3f}")

# Run the analysis
print("🔬 SPIN SYSTEM OPTIMIZATION ANALYSIS")
print("🧲 2D Ising Model with Simulated Annealing")
print("="*60)

# Perform optimization for different parameter sets
results, models = compare_different_parameters()

# Create comprehensive visualization
plot_results(results, models)

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

# Physical constants for Rb-87
class Rb87Constants:
    def __init__(self):
        self.Gamma = 2 * np.pi * 6.067e6  # Natural linewidth (rad/s)
        self.lambda_laser = 780.241e-9    # Transition wavelength (m)
        self.k = 2 * np.pi / self.lambda_laser  # Wave vector
        self.mass = 1.443e-25             # Mass of Rb-87 (kg)
        self.mu_B = 9.274e-24             # Bohr magneton (J/T)
        self.g_J = 2.0                    # Landé g-factor
        self.hbar = 1.055e-34             # Reduced Planck constant

# Initialize constants
rb87 = Rb87Constants()

class MOTSimulation:
    def __init__(self, constants):
        self.const = constants
        
    def scattering_rate(self, detuning, saturation_param):
        """
        Calculate photon scattering rate
        
        Parameters:
        detuning: Laser detuning from resonance (Hz)
        saturation_param: Saturation parameter s0
        
        Returns:
        Scattering rate (photons/s)
        """
        delta_normalized = 2 * detuning / self.const.Gamma
        rate = self.const.Gamma * saturation_param / (1 + saturation_param + delta_normalized**2)
        return rate
    
    def radiation_pressure_force(self, velocity, detuning, saturation_param, magnetic_gradient, position):
        """
        Calculate the total radiation pressure force on an atom
        
        Parameters:
        velocity: Atomic velocity (m/s)
        detuning: Base laser detuning (Hz)
        saturation_param: Saturation parameter
        magnetic_gradient: Magnetic field gradient (T/m)
        position: Atomic position (m)
        
        Returns:
        Force (N)
        """
        # Doppler shift
        doppler_shift = self.const.k * velocity
        
        # Zeeman shift due to magnetic field
        zeeman_shift = self.const.mu_B * self.const.g_J * magnetic_gradient * position / self.const.hbar
        
        # Effective detuning
        eff_detuning = detuning + doppler_shift - zeeman_shift
        
        # Scattering rate
        scatter_rate = self.scattering_rate(eff_detuning, saturation_param)
        
        # Force (assuming counter-propagating beams)
        force = -self.const.hbar * self.const.k * scatter_rate * np.sign(velocity)
        
        return force
    
    def capture_velocity(self, detuning, saturation_param, magnetic_gradient):
        """
        Calculate maximum capture velocity for the MOT
        """
        # Simplified model: maximum velocity that can be slowed to zero
        max_accel = abs(self.radiation_pressure_force(100, detuning, saturation_param, magnetic_gradient, 0.001))
        capture_distance = 0.01  # 1 cm typical MOT size
        
        # Using v^2 = 2*a*d
        v_capture = np.sqrt(2 * max_accel * capture_distance / self.const.mass)
        return v_capture
    
    def trap_depth(self, detuning, saturation_param, magnetic_gradient):
        """
        Calculate effective trap depth (simplified model)
        """
        # Maximum restoring force at trap edge
        trap_radius = 0.005  # 5 mm
        max_force = abs(self.radiation_pressure_force(0, detuning, saturation_param, 
                                                    magnetic_gradient, trap_radius))
        
        # Potential energy at trap edge
        trap_depth = max_force * trap_radius
        return trap_depth
    
    def trapped_atom_number(self, detuning, saturation_param, magnetic_gradient):
        """
        Estimate number of trapped atoms (simplified model)
        
        This combines capture velocity, trap depth, and loading rate
        """
        v_cap = self.capture_velocity(detuning, saturation_param, magnetic_gradient)
        depth = self.trap_depth(detuning, saturation_param, magnetic_gradient)
        
        # Simplified loading model
        loading_rate = v_cap * saturation_param / (1 + saturation_param)
        trap_lifetime = depth * 1e15  # Simplified lifetime model
        
        # Steady-state atom number
        N_atoms = loading_rate * trap_lifetime
        
        # Add some realistic constraints
        if detuning > -0.5 * self.const.Gamma or detuning < -5 * self.const.Gamma:
            N_atoms *= 0.1  # Poor performance outside optimal range
            
        if saturation_param > 10:
            N_atoms *= 0.5  # Heating at high intensity
            
        return N_atoms

# Create simulation instance
mot_sim = MOTSimulation(rb87)

def objective_function(params):
    """
    Objective function to maximize (we minimize the negative)
    params = [detuning (in units of Gamma), saturation_parameter, magnetic_gradient (T/m)]
    """
    detuning_gamma, sat_param, mag_grad = params
    detuning_hz = detuning_gamma * rb87.Gamma / (2 * np.pi)
    
    # Calculate negative atom number (for minimization)
    atom_number = mot_sim.trapped_atom_number(detuning_hz, sat_param, mag_grad)
    return -atom_number

# Define parameter ranges for optimization
bounds = [
    (-5, -0.5),    # Detuning in units of Gamma (red-detuned)
    (0.1, 10),     # Saturation parameter
    (0.01, 0.5)    # Magnetic gradient (T/m)
]

# Initial guess
initial_params = [-2.5, 3.0, 0.1]

print("Starting MOT optimization...")
print(f"Initial parameters: Δ/Γ = {initial_params[0]:.2f}, s₀ = {initial_params[1]:.2f}, B' = {initial_params[2]:.3f} T/m")

# Perform optimization
result = minimize(objective_function, initial_params, bounds=bounds, method='L-BFGS-B')

optimal_params = result.x
optimal_detuning_gamma = optimal_params[0]
optimal_sat_param = optimal_params[1]
optimal_mag_grad = optimal_params[2]

print("\nOptimization Results:")
print(f"Optimal detuning: Δ/Γ = {optimal_detuning_gamma:.3f}")
print(f"Optimal saturation parameter: s₀ = {optimal_sat_param:.3f}")
print(f"Optimal magnetic gradient: B' = {optimal_mag_grad:.3f} T/m")
print(f"Maximum trapped atoms: {-result.fun:.2e}")

# Calculate performance metrics at optimal parameters
optimal_detuning_hz = optimal_detuning_gamma * rb87.Gamma / (2 * np.pi)
v_capture = mot_sim.capture_velocity(optimal_detuning_hz, optimal_sat_param, optimal_mag_grad)
trap_depth = mot_sim.trap_depth(optimal_detuning_hz, optimal_sat_param, optimal_mag_grad)

print(f"\nPerformance at optimal parameters:")
print(f"Capture velocity: {v_capture:.1f} m/s")
print(f"Trap depth: {trap_depth:.2e} J")
print(f"Temperature equivalent: {trap_depth/1.38e-23*1e6:.1f} µK")

# Create detailed parameter sweep for visualization
detuning_range = np.linspace(-4, -1, 50)
sat_param_range = np.linspace(0.5, 8, 50)
mag_grad_fixed = optimal_mag_grad

print("\nGenerating parameter sweep data...")

# 2D parameter sweep
D, S = np.meshgrid(detuning_range, sat_param_range)
atom_numbers = np.zeros_like(D)

for i in range(len(sat_param_range)):
    for j in range(len(detuning_range)):
        detuning_hz = D[i,j] * rb87.Gamma / (2 * np.pi)
        atom_numbers[i,j] = mot_sim.trapped_atom_number(detuning_hz, S[i,j], mag_grad_fixed)

print("Parameter sweep completed. Creating visualizations...")

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

# Plot 1: 2D contour plot of atom number vs detuning and saturation parameter
ax1 = plt.subplot(2, 3, 1)
contour = plt.contourf(D, S, atom_numbers, levels=50, cmap='viridis')
plt.colorbar(contour, label='Trapped Atoms')
plt.plot(optimal_detuning_gamma, optimal_sat_param, 'r*', markersize=15, label='Optimal Point')
plt.xlabel('Detuning (Δ/Γ)')
plt.ylabel('Saturation Parameter (s₀)')
plt.title('MOT Performance Map\n(Fixed B\' = {:.3f} T/m)'.format(optimal_mag_grad))
plt.legend()
plt.grid(True, alpha=0.3)

# Plot 2: Cross-section at optimal saturation parameter
ax2 = plt.subplot(2, 3, 2)
fixed_sat_atoms = [mot_sim.trapped_atom_number(d * rb87.Gamma / (2 * np.pi), optimal_sat_param, optimal_mag_grad) 
                   for d in detuning_range]
plt.plot(detuning_range, fixed_sat_atoms, 'b-', linewidth=2, label=f's₀ = {optimal_sat_param:.2f}')
plt.axvline(optimal_detuning_gamma, color='r', linestyle='--', label='Optimal Δ/Γ')
plt.xlabel('Detuning (Δ/Γ)')
plt.ylabel('Trapped Atoms')
plt.title('Atom Number vs Detuning')
plt.legend()
plt.grid(True, alpha=0.3)

# Plot 3: Cross-section at optimal detuning
ax3 = plt.subplot(2, 3, 3)
fixed_det_atoms = [mot_sim.trapped_atom_number(optimal_detuning_hz, s, optimal_mag_grad) 
                   for s in sat_param_range]
plt.plot(sat_param_range, fixed_det_atoms, 'g-', linewidth=2, label=f'Δ/Γ = {optimal_detuning_gamma:.2f}')
plt.axvline(optimal_sat_param, color='r', linestyle='--', label='Optimal s₀')
plt.xlabel('Saturation Parameter (s₀)')
plt.ylabel('Trapped Atoms')
plt.title('Atom Number vs Intensity')
plt.legend()
plt.grid(True, alpha=0.3)

# Plot 4: Magnetic gradient dependence
ax4 = plt.subplot(2, 3, 4)
mag_grad_range = np.linspace(0.01, 0.3, 100)
mag_grad_atoms = [mot_sim.trapped_atom_number(optimal_detuning_hz, optimal_sat_param, bg) 
                  for bg in mag_grad_range]
plt.plot(mag_grad_range, mag_grad_atoms, 'm-', linewidth=2)
plt.axvline(optimal_mag_grad, color='r', linestyle='--', label='Optimal B\'')
plt.xlabel('Magnetic Gradient (T/m)')
plt.ylabel('Trapped Atoms')
plt.title('Atom Number vs Magnetic Gradient')
plt.legend()
plt.grid(True, alpha=0.3)

# Plot 5: Scattering rate vs detuning for different intensities
ax5 = plt.subplot(2, 3, 5)
detuning_fine = np.linspace(-5, 0, 200)
for s_val in [0.5, 1.0, 2.0, 5.0]:
    scatter_rates = [mot_sim.scattering_rate(d * rb87.Gamma / (2 * np.pi), s_val) 
                    for d in detuning_fine]
    plt.plot(detuning_fine, np.array(scatter_rates) / rb87.Gamma, 
             label=f's₀ = {s_val}', linewidth=2)

plt.axvline(optimal_detuning_gamma, color='r', linestyle='--', alpha=0.7, label='Optimal Δ/Γ')
plt.xlabel('Detuning (Δ/Γ)')
plt.ylabel('Scattering Rate (Γ)')
plt.title('Scattering Rate vs Detuning')
plt.legend()
plt.grid(True, alpha=0.3)

# Plot 6: 3D surface plot (smaller)
ax6 = plt.subplot(2, 3, 6, projection='3d')
D_coarse, S_coarse = np.meshgrid(np.linspace(-4, -1, 20), np.linspace(0.5, 8, 20))
atoms_coarse = np.zeros_like(D_coarse)

for i in range(D_coarse.shape[0]):
    for j in range(D_coarse.shape[1]):
        detuning_hz = D_coarse[i,j] * rb87.Gamma / (2 * np.pi)
        atoms_coarse[i,j] = mot_sim.trapped_atom_number(detuning_hz, S_coarse[i,j], optimal_mag_grad)

surf = ax6.plot_surface(D_coarse, S_coarse, atoms_coarse, cmap='viridis', alpha=0.8)
ax6.scatter([optimal_detuning_gamma], [optimal_sat_param], [-result.fun], 
           color='red', s=100, label='Optimal Point')
ax6.set_xlabel('Detuning (Δ/Γ)')
ax6.set_ylabel('Saturation Parameter (s₀)')
ax6.set_zlabel('Trapped Atoms')
ax6.set_title('3D Performance Surface')

plt.tight_layout()
plt.show()

print("\nAnalysis Summary:")
print("="*50)
print(f"The optimization reveals that the MOT performs best with:")
print(f"• Red detuning of {abs(optimal_detuning_gamma):.2f}Γ")
print(f"• Moderate laser intensity (s₀ = {optimal_sat_param:.2f})")
print(f"• Magnetic gradient of {optimal_mag_grad:.3f} T/m")
print(f"\nThis configuration provides:")
print(f"• Capture velocity: {v_capture:.1f} m/s")
print(f"• Effective temperature: {trap_depth/1.38e-23*1e6:.0f} µK")
print(f"• Maximum trapped atoms: {-result.fun:.1e}")

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

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

class QuantumControlOptimizer:
    """
    A class for optimizing quantum state control of a two-level system (qubit).
    
    The system Hamiltonian is H(t) = (ω₀/2)σz + (u(t)/2)σx
    where u(t) is the control field we want to optimize.
    """
    
    def __init__(self, omega0=1.0, T=5.0, N=100, lambda_reg=0.01):
        """
        Initialize the quantum control optimizer.
        
        Parameters:
        - omega0: Natural frequency of the qubit
        - T: Total evolution time
        - N: Number of time steps
        - lambda_reg: Regularization parameter for control amplitude
        """
        self.omega0 = omega0
        self.T = T
        self.N = N
        self.dt = T / N
        self.lambda_reg = lambda_reg
        
        # Pauli matrices
        self.sigma_x = np.array([[0, 1], [1, 0]], dtype=complex)
        self.sigma_z = np.array([[1, 0], [0, -1]], dtype=complex)
        
        # Initial and target states
        self.psi_initial = np.array([1, 0], dtype=complex)  # |0⟩
        self.psi_target = np.array([0, 1], dtype=complex)   # |1⟩
        
        # Time grid
        self.time_grid = np.linspace(0, T, N+1)
        
    def hamiltonian(self, u_t):
        """
        Construct the Hamiltonian for a given control field value.
        
        H(t) = (ω₀/2)σz + (u(t)/2)σx
        """
        return (self.omega0/2) * self.sigma_z + (u_t/2) * self.sigma_x
    
    def time_evolution_operator(self, control_sequence):
        """
        Compute the time evolution operator for the entire control sequence.
        
        U = T exp[-i ∫₀ᵀ H(t)dt] ≈ ∏ᵢ exp[-i H(tᵢ) Δt]
        """
        U = np.eye(2, dtype=complex)
        
        for i in range(len(control_sequence)):
            H_t = self.hamiltonian(control_sequence[i])
            U_step = expm(-1j * H_t * self.dt)
            U = U_step @ U
            
        return U
    
    def fidelity(self, control_sequence):
        """
        Compute the fidelity between the evolved state and target state.
        
        F = |⟨ψ_target|U|ψ_initial⟩|²
        """
        U = self.time_evolution_operator(control_sequence)
        psi_final = U @ self.psi_initial
        
        # Fidelity is the squared overlap
        overlap = np.vdot(self.psi_target, psi_final)
        return np.abs(overlap)**2
    
    def objective_function(self, control_sequence):
        """
        Objective function to minimize: -Fidelity + λ∫|u(t)|²dt
        
        We want to maximize fidelity while minimizing control effort.
        """
        fid = self.fidelity(control_sequence)
        control_effort = np.sum(control_sequence**2) * self.dt
        
        return -fid + self.lambda_reg * control_effort
    
    def optimize_control(self, initial_guess=None, method='BFGS'):
        """
        Optimize the control sequence using scipy.optimize.minimize.
        """
        if initial_guess is None:
            # Start with a simple sinusoidal guess
            initial_guess = 2 * np.sin(2 * np.pi * self.time_grid[:-1] / self.T)
        
        # Optimization bounds (optional)
        bounds = [(-10, 10) for _ in range(self.N)]
        
        result = minimize(
            self.objective_function,
            initial_guess,
            method=method,
            bounds=bounds,
            options={'maxiter': 1000, 'disp': True}
        )
        
        return result
    
    def simulate_evolution(self, control_sequence):
        """
        Simulate the quantum state evolution under the control sequence.
        Returns the state vector at each time step.
        """
        states = []
        psi = self.psi_initial.copy()
        states.append(psi.copy())
        
        for i in range(len(control_sequence)):
            H_t = self.hamiltonian(control_sequence[i])
            U_step = expm(-1j * H_t * self.dt)
            psi = U_step @ psi
            states.append(psi.copy())
            
        return np.array(states)

# Initialize the optimizer
optimizer = QuantumControlOptimizer(omega0=2.0, T=3.0, N=150, lambda_reg=0.001)

print("=== Quantum State Control Optimization ===")
print(f"System parameters:")
print(f"- Natural frequency ω₀: {optimizer.omega0}")
print(f"- Evolution time T: {optimizer.T}")
print(f"- Time steps N: {optimizer.N}")
print(f"- Regularization λ: {optimizer.lambda_reg}")
print(f"- Initial state: |0⟩")
print(f"- Target state: |1⟩")
print()

# Initial guess: simple sine wave
initial_control = 3 * np.sin(4 * np.pi * optimizer.time_grid[:-1] / optimizer.T)

print("Testing initial guess...")
initial_fidelity = optimizer.fidelity(initial_control)
print(f"Initial fidelity: {initial_fidelity:.6f}")
print()

# Optimize the control sequence
print("Starting optimization...")
result = optimizer.optimize_control(initial_guess=initial_control)

optimal_control = result.x
optimal_fidelity = optimizer.fidelity(optimal_control)

print(f"\nOptimization completed!")
print(f"- Success: {result.success}")
print(f"- Final fidelity: {optimal_fidelity:.6f}")
print(f"- Objective function value: {result.fun:.6f}")
print(f"- Number of iterations: {result.nit}")

# Simulate the evolution with optimal control
optimal_states = optimizer.simulate_evolution(optimal_control)

# Also simulate with initial guess for comparison
initial_states = optimizer.simulate_evolution(initial_control)

print(f"\nFidelity comparison:")
print(f"- Initial guess: {initial_fidelity:.6f}")
print(f"- Optimized control: {optimal_fidelity:.6f}")
print(f"- Improvement: {optimal_fidelity - initial_fidelity:.6f}")

# Create comprehensive plots
fig = plt.figure(figsize=(16, 12))

# Plot 1: Control sequences comparison
ax1 = plt.subplot(2, 3, 1)
plt.plot(optimizer.time_grid[:-1], initial_control, 'b--', linewidth=2, 
         label='Initial guess', alpha=0.7)
plt.plot(optimizer.time_grid[:-1], optimal_control, 'r-', linewidth=2, 
         label='Optimized control')
plt.xlabel('Time')
plt.ylabel('Control amplitude u(t)')
plt.title('Control Field Comparison')
plt.legend()
plt.grid(True, alpha=0.3)

# Plot 2: State populations evolution (initial guess)
ax2 = plt.subplot(2, 3, 2)
prob_0_initial = np.abs(initial_states[:, 0])**2
prob_1_initial = np.abs(initial_states[:, 1])**2

plt.plot(optimizer.time_grid, prob_0_initial, 'b-', linewidth=2, label='|0⟩ population')
plt.plot(optimizer.time_grid, prob_1_initial, 'r-', linewidth=2, label='|1⟩ population')
plt.xlabel('Time')
plt.ylabel('Population')
plt.title('State Evolution (Initial Guess)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(0, 1)

# Plot 3: State populations evolution (optimized)
ax3 = plt.subplot(2, 3, 3)
prob_0_optimal = np.abs(optimal_states[:, 0])**2
prob_1_optimal = np.abs(optimal_states[:, 1])**2

plt.plot(optimizer.time_grid, prob_0_optimal, 'b-', linewidth=2, label='|0⟩ population')
plt.plot(optimizer.time_grid, prob_1_optimal, 'r-', linewidth=2, label='|1⟩ population')
plt.xlabel('Time')
plt.ylabel('Population')
plt.title('State Evolution (Optimized)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(0, 1)

# Plot 4: Bloch sphere representation (initial guess)
ax4 = plt.subplot(2, 3, 4, projection='3d')

# Calculate Bloch vector components for initial guess
x_initial = 2 * np.real(initial_states[:, 0] * np.conj(initial_states[:, 1]))
y_initial = -2 * np.imag(initial_states[:, 0] * np.conj(initial_states[:, 1]))
z_initial = np.abs(initial_states[:, 0])**2 - np.abs(initial_states[:, 1])**2

# Plot trajectory
ax4.plot(x_initial, y_initial, z_initial, 'b-', linewidth=2, alpha=0.7)
ax4.scatter(x_initial[0], y_initial[0], z_initial[0], color='green', s=100, label='Start')
ax4.scatter(x_initial[-1], y_initial[-1], z_initial[-1], color='red', s=100, label='End')

# Draw unit sphere
u = np.linspace(0, 2 * np.pi, 50)
v = np.linspace(0, np.pi, 50)
x_sphere = np.outer(np.cos(u), np.sin(v))
y_sphere = np.outer(np.sin(u), np.sin(v))
z_sphere = np.outer(np.ones(np.size(u)), np.cos(v))
ax4.plot_surface(x_sphere, y_sphere, z_sphere, alpha=0.1, color='gray')

ax4.set_xlabel('X')
ax4.set_ylabel('Y')
ax4.set_zlabel('Z')
ax4.set_title('Bloch Sphere (Initial)')
ax4.legend()

# Plot 5: Bloch sphere representation (optimized)
ax5 = plt.subplot(2, 3, 5, projection='3d')

# Calculate Bloch vector components for optimized control
x_optimal = 2 * np.real(optimal_states[:, 0] * np.conj(optimal_states[:, 1]))
y_optimal = -2 * np.imag(optimal_states[:, 0] * np.conj(optimal_states[:, 1]))
z_optimal = np.abs(optimal_states[:, 0])**2 - np.abs(optimal_states[:, 1])**2

# Plot trajectory
ax5.plot(x_optimal, y_optimal, z_optimal, 'r-', linewidth=2)
ax5.scatter(x_optimal[0], y_optimal[0], z_optimal[0], color='green', s=100, label='Start')
ax5.scatter(x_optimal[-1], y_optimal[-1], z_optimal[-1], color='red', s=100, label='End')

# Draw unit sphere
ax5.plot_surface(x_sphere, y_sphere, z_sphere, alpha=0.1, color='gray')

ax5.set_xlabel('X')
ax5.set_ylabel('Y')
ax5.set_zlabel('Z')
ax5.set_title('Bloch Sphere (Optimized)')
ax5.legend()

# Plot 6: Fidelity comparison
ax6 = plt.subplot(2, 3, 6)
categories = ['Initial Guess', 'Optimized']
fidelities = [initial_fidelity, optimal_fidelity]
colors = ['skyblue', 'salmon']

bars = plt.bar(categories, fidelities, color=colors, alpha=0.8)
plt.ylabel('Fidelity')
plt.title('Fidelity Comparison')
plt.ylim(0, 1)

# Add value labels on bars
for bar, fid in zip(bars, fidelities):
    height = bar.get_height()
    plt.text(bar.get_x() + bar.get_width()/2., height + 0.01,
             f'{fid:.4f}', ha='center', va='bottom', fontweight='bold')

plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Additional analysis: Control effort comparison
initial_effort = np.sum(initial_control**2) * optimizer.dt
optimal_effort = np.sum(optimal_control**2) * optimizer.dt

print(f"\nControl effort analysis:")
print(f"- Initial guess effort: {initial_effort:.4f}")
print(f"- Optimized control effort: {optimal_effort:.4f}")
print(f"- Effort ratio (opt/initial): {optimal_effort/initial_effort:.4f}")

# Frequency analysis of control pulses
from scipy.fft import fft, fftfreq

fig2, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

# FFT of initial control
freqs = fftfreq(optimizer.N, optimizer.dt)
initial_fft = np.abs(fft(initial_control))
optimal_fft = np.abs(fft(optimal_control))

ax1.plot(freqs[:optimizer.N//2], initial_fft[:optimizer.N//2], 'b-', 
         linewidth=2, label='Initial guess')
ax1.set_xlabel('Frequency')
ax1.set_ylabel('Amplitude')
ax1.set_title('Frequency Spectrum (Initial)')
ax1.grid(True, alpha=0.3)

ax2.plot(freqs[:optimizer.N//2], optimal_fft[:optimizer.N//2], 'r-', 
         linewidth=2, label='Optimized')
ax2.set_xlabel('Frequency')
ax2.set_ylabel('Amplitude')
ax2.set_title('Frequency Spectrum (Optimized)')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("\n=== Analysis Summary ===")
print(f"The optimization successfully improved the fidelity from {initial_fidelity:.4f} to {optimal_fidelity:.4f}")
print(f"The optimal control achieves {optimal_fidelity*100:.2f}% fidelity in transferring |0⟩ → |1⟩")
print(f"Final |1⟩ population: {prob_1_optimal[-1]:.4f}")
print(f"The control effort was {'reduced' if optimal_effort < initial_effort else 'increased'} by optimization")