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

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

class NonlinearPortfolioOptimizer:
    def __init__(self, n_assets=5, n_simulations=10000):
        self.n_assets = n_assets
        self.n_simulations = n_simulations
        self.setup_market_data()
        
    def setup_market_data(self):
        """Generate synthetic market data with nonlinear correlations"""
        # Expected returns (annualized)
        self.mu = np.array([0.08, 0.12, 0.15, 0.10, 0.06])
        
        # Asset names
        self.asset_names = ['Stock A', 'Stock B', 'Stock C', 'Bond A', 'Bond B']
        
        # Base correlation matrix
        base_corr = np.array([
            [1.0, 0.3, 0.5, -0.2, -0.1],
            [0.3, 1.0, 0.7, -0.1, 0.0],
            [0.5, 0.7, 1.0, 0.1, 0.2],
            [-0.2, -0.1, 0.1, 1.0, 0.8],
            [-0.1, 0.0, 0.2, 0.8, 1.0]
        ])
        
        # Volatilities (annualized)
        self.sigma = np.array([0.20, 0.25, 0.30, 0.08, 0.06])
        
        # Create covariance matrix
        self.cov_matrix = np.outer(self.sigma, self.sigma) * base_corr
        
        # Transaction cost coefficients
        self.transaction_costs = np.array([0.002, 0.003, 0.004, 0.001, 0.001])
        
        # Generate nonlinear correlation scenario returns
        self.generate_nonlinear_returns()
        
    def generate_nonlinear_returns(self):
        """Generate returns with nonlinear correlation structure"""
        # Generate base multivariate normal returns
        base_returns = multivariate_normal.rvs(
            mean=self.mu/252,  # Daily returns
            cov=self.cov_matrix/252,
            size=self.n_simulations
        )
        
        # Add nonlinear correlation effects
        # Regime-dependent correlations based on market stress
        market_factor = base_returns[:, :3].mean(axis=1)  # Average stock performance
        stress_indicator = market_factor < np.percentile(market_factor, 25)
        
        # Increase correlations during stress periods
        stressed_returns = base_returns.copy()
        for i in range(self.n_simulations):
            if stress_indicator[i]:
                # Increase correlations between stocks during stress
                factor = 1.5
                stressed_returns[i, :3] = (
                    base_returns[i, :3] * 0.7 + 
                    market_factor[i] * 0.3 * factor
                )
        
        self.return_scenarios = stressed_returns
        
    def nonlinear_transaction_cost(self, weights):
        """Calculate nonlinear transaction costs"""
        return np.sum(self.transaction_costs * np.power(weights, 1.5))
    
    def portfolio_var(self, weights, confidence_level=0.05):
        """Calculate portfolio VaR using Monte Carlo simulation"""
        portfolio_returns = np.dot(self.return_scenarios, weights)
        var = -np.percentile(portfolio_returns, confidence_level * 100)
        return var
    
    def portfolio_expected_return(self, weights):
        """Calculate expected portfolio return"""
        return np.dot(weights, self.mu)
    
    def objective_function(self, weights):
        """Objective function: maximize expected return minus transaction costs"""
        expected_return = self.portfolio_expected_return(weights)
        transaction_cost = self.nonlinear_transaction_cost(weights)
        return -(expected_return - transaction_cost)  # Negative for minimization
    
    def var_constraint(self, weights, max_var):
        """VaR constraint function"""
        return max_var - self.portfolio_var(weights)
    
    def optimize_portfolio(self, max_var=0.02):
        """Solve the portfolio optimization problem"""
        # Initial guess (equal weights)
        x0 = np.ones(self.n_assets) / self.n_assets
        
        # Constraints
        constraints = [
            {'type': 'eq', 'fun': lambda x: np.sum(x) - 1.0},  # Sum to 1
            {'type': 'ineq', 'fun': lambda x: self.var_constraint(x, max_var)}  # VaR constraint
        ]
        
        # Bounds (long-only)
        bounds = [(0, 1) for _ in range(self.n_assets)]
        
        # Solve optimization
        result = optimize.minimize(
            self.objective_function,
            x0,
            method='SLSQP',
            bounds=bounds,
            constraints=constraints,
            options={'disp': True, 'maxiter': 1000}
        )
        
        return result
    
    def efficient_frontier(self, var_levels):
        """Generate efficient frontier for different VaR levels"""
        results = []
        for var_level in var_levels:
            try:
                result = self.optimize_portfolio(max_var=var_level)
                if result.success:
                    weights = result.x
                    expected_return = self.portfolio_expected_return(weights)
                    actual_var = self.portfolio_var(weights)
                    transaction_cost = self.nonlinear_transaction_cost(weights)
                    
                    results.append({
                        'var_limit': var_level,
                        'expected_return': expected_return,
                        'actual_var': actual_var,
                        'transaction_cost': transaction_cost,
                        'net_return': expected_return - transaction_cost,
                        'weights': weights
                    })
            except:
                continue
        
        return pd.DataFrame(results)

# Initialize optimizer
optimizer = NonlinearPortfolioOptimizer()

# Solve for a specific VaR constraint
print("=== Portfolio Optimization with Nonlinear Correlations ===")
print(f"Assets: {optimizer.asset_names}")
print(f"Expected Returns: {optimizer.mu}")
print(f"Volatilities: {optimizer.sigma}")
print(f"Transaction Cost Coefficients: {optimizer.transaction_costs}")

# Optimize with 2% VaR limit
result = optimizer.optimize_portfolio(max_var=0.02)

if result.success:
    optimal_weights = result.x
    print(f"\n=== Optimization Results ===")
    print(f"Optimization successful: {result.success}")
    print(f"Optimal weights:")
    for i, (asset, weight) in enumerate(zip(optimizer.asset_names, optimal_weights)):
        print(f"  {asset}: {weight:.4f} ({weight*100:.2f}%)")
    
    # Calculate portfolio metrics
    expected_return = optimizer.portfolio_expected_return(optimal_weights)
    portfolio_var = optimizer.portfolio_var(optimal_weights)
    transaction_cost = optimizer.nonlinear_transaction_cost(optimal_weights)
    net_return = expected_return - transaction_cost
    
    print(f"\n=== Portfolio Metrics ===")
    print(f"Expected Return: {expected_return:.4f} ({expected_return*100:.2f}%)")
    print(f"Transaction Cost: {transaction_cost:.4f} ({transaction_cost*100:.2f}%)")
    print(f"Net Expected Return: {net_return:.4f} ({net_return*100:.2f}%)")
    print(f"Portfolio VaR (5%): {portfolio_var:.4f} ({portfolio_var*100:.2f}%)")
else:
    print("Optimization failed!")
    print(result.message)

# Generate efficient frontier
print("\n=== Generating Efficient Frontier ===")
var_levels = np.linspace(0.01, 0.05, 20)
efficient_frontier_data = optimizer.efficient_frontier(var_levels)

print(f"Successfully computed {len(efficient_frontier_data)} frontier points")

# Create comprehensive visualizations
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Advanced Portfolio Optimization Analysis', fontsize=16, fontweight='bold')

# 1. Efficient Frontier
ax1 = axes[0, 0]
if len(efficient_frontier_data) > 0:
    ax1.scatter(efficient_frontier_data['actual_var'] * 100, 
               efficient_frontier_data['net_return'] * 100,
               c='blue', s=50, alpha=0.7, label='Efficient Frontier')
    ax1.scatter(portfolio_var * 100, net_return * 100, 
               c='red', s=100, marker='*', label='Optimal Portfolio')
ax1.set_xlabel('Portfolio VaR (%)')
ax1.set_ylabel('Net Expected Return (%)')
ax1.set_title('Efficient Frontier')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 2. Optimal Portfolio Allocation
ax2 = axes[0, 1]
if result.success:
    colors = plt.cm.Set3(np.linspace(0, 1, len(optimizer.asset_names)))
    wedges, texts, autotexts = ax2.pie(optimal_weights, labels=optimizer.asset_names, 
                                      autopct='%1.1f%%', colors=colors, startangle=90)
    ax2.set_title('Optimal Portfolio Allocation')

# 3. Asset Returns vs Transaction Costs
ax3 = axes[0, 2]
x_pos = np.arange(len(optimizer.asset_names))
width = 0.35
bars1 = ax3.bar(x_pos - width/2, optimizer.mu * 100, width, 
                label='Expected Return (%)', alpha=0.8, color='green')
bars2 = ax3.bar(x_pos + width/2, optimizer.transaction_costs * 100, width,
                label='Transaction Cost Coeff (%)', alpha=0.8, color='red')
ax3.set_xlabel('Assets')
ax3.set_ylabel('Percentage (%)')
ax3.set_title('Expected Returns vs Transaction Costs')
ax3.set_xticks(x_pos)
ax3.set_xticklabels(optimizer.asset_names, rotation=45)
ax3.legend()
ax3.grid(True, alpha=0.3)

# 4. Correlation Heatmap
ax4 = axes[1, 0]
correlation_matrix = optimizer.cov_matrix / np.outer(optimizer.sigma, optimizer.sigma)
im = ax4.imshow(correlation_matrix, cmap='RdBu_r', aspect='auto', vmin=-1, vmax=1)
ax4.set_xticks(range(len(optimizer.asset_names)))
ax4.set_yticks(range(len(optimizer.asset_names)))
ax4.set_xticklabels(optimizer.asset_names, rotation=45)
ax4.set_yticklabels(optimizer.asset_names)
ax4.set_title('Asset Correlation Matrix')

# Add correlation values to heatmap
for i in range(len(optimizer.asset_names)):
    for j in range(len(optimizer.asset_names)):
        text = ax4.text(j, i, f'{correlation_matrix[i, j]:.2f}',
                       ha="center", va="center", color="black", fontsize=8)

plt.colorbar(im, ax=ax4, shrink=0.8)

# 5. VaR vs Expected Return Trade-off
ax5 = axes[1, 1]
if len(efficient_frontier_data) > 0:
    ax5.plot(efficient_frontier_data['actual_var'] * 100, 
             efficient_frontier_data['expected_return'] * 100,
             'b-o', label='Gross Return', markersize=4)
    ax5.plot(efficient_frontier_data['actual_var'] * 100, 
             efficient_frontier_data['net_return'] * 100,
             'r--s', label='Net Return (after costs)', markersize=4)
    ax5.scatter(portfolio_var * 100, expected_return * 100, 
               c='blue', s=100, marker='*', label='Optimal (Gross)')
    ax5.scatter(portfolio_var * 100, net_return * 100, 
               c='red', s=100, marker='*', label='Optimal (Net)')
ax5.set_xlabel('Portfolio VaR (%)')
ax5.set_ylabel('Expected Return (%)')
ax5.set_title('Return vs Risk Trade-off')
ax5.legend(fontsize=8)
ax5.grid(True, alpha=0.3)

# 6. Portfolio Return Distribution
ax6 = axes[1, 2]
if result.success:
    portfolio_returns = np.dot(optimizer.return_scenarios, optimal_weights)
    ax6.hist(portfolio_returns * 100, bins=50, density=True, alpha=0.7, 
             color='skyblue', edgecolor='black')
    ax6.axvline(-portfolio_var * 100, color='red', linestyle='--', 
                label=f'VaR (5%): {portfolio_var*100:.2f}%')
    ax6.axvline(np.mean(portfolio_returns) * 100, color='green', linestyle='-', 
                label=f'Mean: {np.mean(portfolio_returns)*100:.2f}%')
    ax6.set_xlabel('Daily Return (%)')
    ax6.set_ylabel('Density')
    ax6.set_title('Portfolio Return Distribution')
    ax6.legend()
    ax6.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Summary statistics
print("\n=== Summary Analysis ===")
if len(efficient_frontier_data) > 0:
    print(f"Efficient Frontier Statistics:")
    print(f"  Number of efficient portfolios: {len(efficient_frontier_data)}")
    print(f"  VaR range: {efficient_frontier_data['actual_var'].min()*100:.2f}% - {efficient_frontier_data['actual_var'].max()*100:.2f}%")
    print(f"  Net return range: {efficient_frontier_data['net_return'].min()*100:.2f}% - {efficient_frontier_data['net_return'].max()*100:.2f}%")
    print(f"  Average transaction cost: {efficient_frontier_data['transaction_cost'].mean()*100:.3f}%")

if result.success:
    portfolio_returns = np.dot(optimizer.return_scenarios, optimal_weights)
    print(f"\nOptimal Portfolio Risk Metrics:")
    print(f"  VaR (1%): {-np.percentile(portfolio_returns, 1)*100:.3f}%")
    print(f"  VaR (5%): {-np.percentile(portfolio_returns, 5)*100:.3f}%")
    print(f"  VaR (10%): {-np.percentile(portfolio_returns, 10)*100:.3f}%")
    print(f"  Expected Shortfall (5%): {-np.mean(portfolio_returns[portfolio_returns <= np.percentile(portfolio_returns, 5)])*100:.3f}%")
    print(f"  Portfolio Volatility: {np.std(portfolio_returns)*np.sqrt(252)*100:.2f}% (annualized)")
    print(f"  Sharpe Ratio: {(net_return*252)/(np.std(portfolio_returns)*np.sqrt(252)):.3f}")

def setup_market_data(self):

def generate_nonlinear_returns(self):

def nonlinear_transaction_cost(self, weights):
    return np.sum(self.transaction_costs * np.power(weights, 1.5))

def portfolio_var(self, weights, confidence_level=0.05):

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.integrate import solve_ivp
import warnings
warnings.filterwarnings('ignore')

class TwoLinkRobot:
    """
    2-DOF planar robotic arm for trajectory optimization
    """
    def __init__(self, L1=1.0, L2=1.0, m1=1.0, m2=1.0):
        """
        Initialize robot parameters
        L1, L2: link lengths
        m1, m2: link masses
        """
        self.L1 = L1
        self.L2 = L2
        self.m1 = m1
        self.m2 = m2
        self.g = 9.81  # gravity
        
    def forward_kinematics(self, q):
        """
        Compute end-effector position given joint angles
        q: [q1, q2] joint angles
        Returns: [x, y] end-effector position
        """
        q1, q2 = q[0], q[1]
        x = self.L1 * np.cos(q1) + self.L2 * np.cos(q1 + q2)
        y = self.L1 * np.sin(q1) + self.L2 * np.sin(q1 + q2)
        return np.array([x, y])
    
    def mass_matrix(self, q):
        """
        Compute mass matrix M(q)
        """
        q1, q2 = q[0], q[1]
        
        # Mass matrix elements
        M11 = (self.m1 + self.m2) * self.L1**2 + self.m2 * self.L2**2 + 2 * self.m2 * self.L1 * self.L2 * np.cos(q2)
        M12 = self.m2 * self.L2**2 + self.m2 * self.L1 * self.L2 * np.cos(q2)
        M21 = M12
        M22 = self.m2 * self.L2**2
        
        return np.array([[M11, M12], [M21, M22]])
    
    def coriolis_matrix(self, q, dq):
        """
        Compute Coriolis matrix C(q, dq)
        """
        q1, q2 = q[0], q[1]
        dq1, dq2 = dq[0], dq[1]
        
        # Coriolis terms
        h = -self.m2 * self.L1 * self.L2 * np.sin(q2)
        C11 = h * dq2
        C12 = h * (dq1 + dq2)
        C21 = -h * dq1
        C22 = 0
        
        return np.array([[C11, C12], [C21, C22]])
    
    def gravity_vector(self, q):
        """
        Compute gravity vector G(q)
        """
        q1, q2 = q[0], q[1]
        
        G1 = (self.m1 + self.m2) * self.g * self.L1 * np.cos(q1) + self.m2 * self.g * self.L2 * np.cos(q1 + q2)
        G2 = self.m2 * self.g * self.L2 * np.cos(q1 + q2)
        
        return np.array([G1, G2])

class TrajectoryOptimizer:
    """
    Trajectory optimization for 2-DOF robot
    """
    def __init__(self, robot, T=2.0, N=50):
        """
        Initialize optimizer
        robot: TwoLinkRobot instance
        T: time horizon
        N: number of discretization points
        """
        self.robot = robot
        self.T = T
        self.N = N
        self.dt = T / (N - 1)
        
        # Obstacle parameters (circular obstacle)
        self.obstacle_center = np.array([1.0, 0.5])
        self.obstacle_radius = 0.3
        
    def interpolate_trajectory(self, q_waypoints):
        """
        Create smooth trajectory using polynomial interpolation
        """
        t = np.linspace(0, self.T, self.N)
        
        # Use cubic polynomial for smooth trajectory
        q_traj = np.zeros((self.N, 2))
        dq_traj = np.zeros((self.N, 2))
        ddq_traj = np.zeros((self.N, 2))
        
        for i in range(2):  # for each joint
            # Cubic polynomial coefficients
            # q(t) = a0 + a1*t + a2*t^2 + a3*t^3
            # Boundary conditions: q(0) = q_start, q(T) = q_end, dq(0) = 0, dq(T) = 0
            
            q_start = q_waypoints[0, i]
            q_end = q_waypoints[-1, i]
            
            # Solve for coefficients
            A = np.array([[0, 0, 0, 1],
                         [self.T**3, self.T**2, self.T, 1],
                         [0, 0, 1, 0],
                         [3*self.T**2, 2*self.T, 1, 0]])
            b = np.array([q_start, q_end, 0, 0])
            coeffs = np.linalg.solve(A, b)
            
            # Generate trajectory
            q_traj[:, i] = coeffs[3] + coeffs[2]*t + coeffs[1]*t**2 + coeffs[0]*t**3
            dq_traj[:, i] = coeffs[2] + 2*coeffs[1]*t + 3*coeffs[0]*t**2
            ddq_traj[:, i] = 2*coeffs[1] + 6*coeffs[0]*t
            
        return q_traj, dq_traj, ddq_traj
    
    def compute_torques(self, q_traj, dq_traj, ddq_traj):
        """
        Compute required torques using inverse dynamics
        """
        torques = np.zeros((self.N, 2))
        
        for i in range(self.N):
            q = q_traj[i]
            dq = dq_traj[i]
            ddq = ddq_traj[i]
            
            M = self.robot.mass_matrix(q)
            C = self.robot.coriolis_matrix(q, dq)
            G = self.robot.gravity_vector(q)
            
            # Inverse dynamics: tau = M*ddq + C*dq + G
            torques[i] = M @ ddq + C @ dq + G
            
        return torques
    
    def objective_function(self, q_waypoints_flat):
        """
        Objective function for optimization
        """
        # Reshape waypoints
        q_waypoints = q_waypoints_flat.reshape(-1, 2)
        
        # Generate trajectory
        q_traj, dq_traj, ddq_traj = self.interpolate_trajectory(q_waypoints)
        
        # Compute torques
        torques = self.compute_torques(q_traj, dq_traj, ddq_traj)
        
        # Energy cost (integral of torque squared)
        energy_cost = np.sum(torques**2) * self.dt
        
        # Smoothness cost (integral of acceleration squared)
        smoothness_cost = np.sum(ddq_traj**2) * self.dt
        
        # Total cost
        total_cost = energy_cost + 0.1 * smoothness_cost
        
        return total_cost
    
    def obstacle_constraint(self, q_waypoints_flat):
        """
        Obstacle avoidance constraint
        """
        q_waypoints = q_waypoints_flat.reshape(-1, 2)
        q_traj, _, _ = self.interpolate_trajectory(q_waypoints)
        
        constraints = []
        for i in range(self.N):
            # End-effector position
            ee_pos = self.robot.forward_kinematics(q_traj[i])
            
            # Distance to obstacle center
            dist = np.linalg.norm(ee_pos - self.obstacle_center)
            
            # Constraint: distance should be greater than obstacle radius
            constraints.append(dist - self.obstacle_radius - 0.1)  # 0.1m safety margin
            
        return np.array(constraints)
    
    def optimize_trajectory(self, q_start, q_end, n_waypoints=5):
        """
        Optimize trajectory from start to end configuration
        """
        # Initialize waypoints with linear interpolation
        t_waypoints = np.linspace(0, 1, n_waypoints)
        q_waypoints_init = np.zeros((n_waypoints, 2))
        
        for i in range(2):
            q_waypoints_init[:, i] = q_start[i] + t_waypoints * (q_end[i] - q_start[i])
        
        # Fix start and end points
        q_waypoints_init[0] = q_start
        q_waypoints_init[-1] = q_end
        
        # Flatten for optimization
        x0 = q_waypoints_init.flatten()
        
        # Bounds (joint limits)
        bounds = []
        for i in range(n_waypoints * 2):
            if i < 2 or i >= (n_waypoints - 1) * 2:  # Fix start and end points
                bounds.append((x0[i], x0[i]))
            else:
                bounds.append((-np.pi, np.pi))
        
        # Constraints
        constraints = {'type': 'ineq', 'fun': self.obstacle_constraint}
        
        # Optimize
        result = minimize(self.objective_function, x0, method='SLSQP', 
                         bounds=bounds, constraints=constraints,
                         options={'maxiter': 1000, 'ftol': 1e-6})
        
        # Extract optimized waypoints
        q_waypoints_opt = result.x.reshape(-1, 2)
        
        return q_waypoints_opt, result

# Initialize robot and optimizer
robot = TwoLinkRobot(L1=1.0, L2=1.0, m1=1.0, m2=1.0)
optimizer = TrajectoryOptimizer(robot, T=3.0, N=100)

# Define start and end configurations
q_start = np.array([0.0, 0.0])  # Straight configuration
q_end = np.array([np.pi/2, -np.pi/4])  # Target configuration

print("Starting trajectory optimization...")
print(f"Start configuration: q1={q_start[0]:.3f}, q2={q_start[1]:.3f}")
print(f"End configuration: q1={q_end[0]:.3f}, q2={q_end[1]:.3f}")

# Optimize trajectory
q_waypoints_opt, opt_result = optimizer.optimize_trajectory(q_start, q_end, n_waypoints=8)

print(f"\nOptimization completed!")
print(f"Success: {opt_result.success}")
print(f"Final cost: {opt_result.fun:.6f}")
print(f"Number of iterations: {opt_result.nit}")

# Generate final optimized trajectory
q_traj_opt, dq_traj_opt, ddq_traj_opt = optimizer.interpolate_trajectory(q_waypoints_opt)
torques_opt = optimizer.compute_torques(q_traj_opt, dq_traj_opt, ddq_traj_opt)

# Generate comparison trajectory (straight line in joint space)
q_waypoints_straight = np.array([q_start, q_end])
q_traj_straight, dq_traj_straight, ddq_traj_straight = optimizer.interpolate_trajectory(q_waypoints_straight)
torques_straight = optimizer.compute_torques(q_traj_straight, dq_traj_straight, ddq_traj_straight)

print("\nTrajectory generation completed!")
print("Ready for visualization...")

# Visualization
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Robot Trajectory Optimization Results', fontsize=16, fontweight='bold')

# Time vector
t = np.linspace(0, optimizer.T, optimizer.N)

# Plot 1: Joint trajectories
axes[0, 0].plot(t, q_traj_straight[:, 0], 'r--', label='Straight - Joint 1', linewidth=2)
axes[0, 0].plot(t, q_traj_straight[:, 1], 'b--', label='Straight - Joint 2', linewidth=2)
axes[0, 0].plot(t, q_traj_opt[:, 0], 'r-', label='Optimized - Joint 1', linewidth=2)
axes[0, 0].plot(t, q_traj_opt[:, 1], 'b-', label='Optimized - Joint 2', linewidth=2)
axes[0, 0].set_xlabel('Time [s]')
axes[0, 0].set_ylabel('Joint Angle [rad]')
axes[0, 0].set_title('Joint Trajectories')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Plot 2: Joint velocities
axes[0, 1].plot(t, dq_traj_straight[:, 0], 'r--', label='Straight - Joint 1', linewidth=2)
axes[0, 1].plot(t, dq_traj_straight[:, 1], 'b--', label='Straight - Joint 2', linewidth=2)
axes[0, 1].plot(t, dq_traj_opt[:, 0], 'r-', label='Optimized - Joint 1', linewidth=2)
axes[0, 1].plot(t, dq_traj_opt[:, 1], 'b-', label='Optimized - Joint 2', linewidth=2)
axes[0, 1].set_xlabel('Time [s]')
axes[0, 1].set_ylabel('Joint Velocity [rad/s]')
axes[0, 1].set_title('Joint Velocities')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Plot 3: Required torques
axes[0, 2].plot(t, torques_straight[:, 0], 'r--', label='Straight - Joint 1', linewidth=2)
axes[0, 2].plot(t, torques_straight[:, 1], 'b--', label='Straight - Joint 2', linewidth=2)
axes[0, 2].plot(t, torques_opt[:, 0], 'r-', label='Optimized - Joint 1', linewidth=2)
axes[0, 2].plot(t, torques_opt[:, 1], 'b-', label='Optimized - Joint 2', linewidth=2)
axes[0, 2].set_xlabel('Time [s]')
axes[0, 2].set_ylabel('Torque [Nm]')
axes[0, 2].set_title('Required Torques')
axes[0, 2].legend()
axes[0, 2].grid(True, alpha=0.3)

# Plot 4: End-effector trajectory in Cartesian space
ee_traj_straight = np.array([robot.forward_kinematics(q) for q in q_traj_straight])
ee_traj_opt = np.array([robot.forward_kinematics(q) for q in q_traj_opt])

axes[1, 0].plot(ee_traj_straight[:, 0], ee_traj_straight[:, 1], 'r--', label='Straight', linewidth=3)
axes[1, 0].plot(ee_traj_opt[:, 0], ee_traj_opt[:, 1], 'g-', label='Optimized', linewidth=3)

# Plot obstacle
circle = plt.Circle(optimizer.obstacle_center, optimizer.obstacle_radius, 
                   color='red', alpha=0.3, label='Obstacle')
axes[1, 0].add_patch(circle)

# Mark start and end points
axes[1, 0].plot(ee_traj_straight[0, 0], ee_traj_straight[0, 1], 'ko', markersize=8, label='Start')
axes[1, 0].plot(ee_traj_straight[-1, 0], ee_traj_straight[-1, 1], 'ks', markersize=8, label='End')

axes[1, 0].set_xlabel('X [m]')
axes[1, 0].set_ylabel('Y [m]')
axes[1, 0].set_title('End-Effector Trajectory')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].set_aspect('equal')

# Plot 5: Energy consumption comparison
energy_straight = np.cumsum(np.sum(torques_straight**2, axis=1)) * optimizer.dt
energy_opt = np.cumsum(np.sum(torques_opt**2, axis=1)) * optimizer.dt

axes[1, 1].plot(t, energy_straight, 'r--', label='Straight', linewidth=3)
axes[1, 1].plot(t, energy_opt, 'g-', label='Optimized', linewidth=3)
axes[1, 1].set_xlabel('Time [s]')
axes[1, 1].set_ylabel('Cumulative Energy [J]')
axes[1, 1].set_title('Energy Consumption')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

# Plot 6: Robot arm animation (final positions)
def plot_robot_arm(ax, q, color, label, alpha=1.0):
    # Joint positions
    x0, y0 = 0, 0  # Base
    x1 = robot.L1 * np.cos(q[0])
    y1 = robot.L1 * np.sin(q[0])
    x2 = x1 + robot.L2 * np.cos(q[0] + q[1])
    y2 = y1 + robot.L2 * np.sin(q[0] + q[1])
    
    # Plot links
    ax.plot([x0, x1], [y0, y1], color=color, linewidth=4, alpha=alpha, label=f'{label} Link 1')
    ax.plot([x1, x2], [y1, y2], color=color, linewidth=4, alpha=alpha, label=f'{label} Link 2')
    
    # Plot joints
    ax.plot([x0, x1, x2], [y0, y1, y2], 'o', color=color, markersize=8, alpha=alpha)
    
    return x2, y2

# Plot initial and final configurations
plot_robot_arm(axes[1, 2], q_start, 'blue', 'Initial', alpha=0.5)
plot_robot_arm(axes[1, 2], q_end, 'red', 'Final', alpha=0.8)

# Plot obstacle
circle2 = plt.Circle(optimizer.obstacle_center, optimizer.obstacle_radius, 
                    color='red', alpha=0.3)
axes[1, 2].add_patch(circle2)

axes[1, 2].set_xlabel('X [m]')
axes[1, 2].set_ylabel('Y [m]')
axes[1, 2].set_title('Robot Arm Configurations')
axes[1, 2].grid(True, alpha=0.3)
axes[1, 2].set_aspect('equal')
axes[1, 2].set_xlim(-0.5, 2.5)
axes[1, 2].set_ylim(-0.5, 2.5)

plt.tight_layout()
plt.show()

# Print performance comparison
total_energy_straight = energy_straight[-1]
total_energy_opt = energy_opt[-1]
energy_saving = (total_energy_straight - total_energy_opt) / total_energy_straight * 100

print(f"\n=== Performance Comparison ===")
print(f"Total Energy - Straight trajectory: {total_energy_straight:.4f} J")
print(f"Total Energy - Optimized trajectory: {total_energy_opt:.4f} J")
print(f"Energy saving: {energy_saving:.2f}%")

# Check obstacle avoidance
min_distance_straight = min([np.linalg.norm(robot.forward_kinematics(q) - optimizer.obstacle_center) 
                           for q in q_traj_straight])
min_distance_opt = min([np.linalg.norm(robot.forward_kinematics(q) - optimizer.obstacle_center) 
                       for q in q_traj_opt])

print(f"\nMinimum distance to obstacle:")
print(f"Straight trajectory: {min_distance_straight:.4f} m")
print(f"Optimized trajectory: {min_distance_opt:.4f} m")
print(f"Obstacle radius: {optimizer.obstacle_radius:.4f} m")

if min_distance_straight < optimizer.obstacle_radius:
    print("⚠️  Straight trajectory collides with obstacle!")
else:
    print("✅ Straight trajectory avoids obstacle")

if min_distance_opt < optimizer.obstacle_radius:
    print("⚠️  Optimized trajectory collides with obstacle!")
else:
    print("✅ Optimized trajectory avoids obstacle")

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

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

class EconomicLoadDispatch:
    """
    Economic Load Dispatch solver for power systems
    """
    
    def __init__(self, generators_data, demand, loss_coefficients=None):
        """
        Initialize the ELD problem
        
        Parameters:
        generators_data: list of dictionaries containing generator parameters
        demand: total power demand (MW)
        loss_coefficients: transmission loss coefficients (optional)
        """
        self.generators = generators_data
        self.demand = demand
        self.n_gen = len(generators_data)
        self.loss_coeffs = loss_coefficients
        
    def cost_function(self, power_output):
        """
        Calculate total generation cost
        
        Parameters:
        power_output: array of power outputs for each generator
        
        Returns:
        total_cost: total fuel cost
        """
        total_cost = 0
        for i, gen in enumerate(self.generators):
            a, b, c = gen['a'], gen['b'], gen['c']
            P = power_output[i]
            total_cost += a * P**2 + b * P + c
        return total_cost
    
    def transmission_loss(self, power_output):
        """
        Calculate transmission losses using B-coefficients method
        
        Parameters:
        power_output: array of power outputs
        
        Returns:
        loss: transmission loss (MW)
        """
        if self.loss_coeffs is None:
            return 0
        
        # Simplified loss calculation: P_loss = sum(B_ii * P_i^2)
        loss = 0
        for i, B_ii in enumerate(self.loss_coeffs):
            loss += B_ii * power_output[i]**2
        return loss
    
    def power_balance_constraint(self, power_output):
        """
        Power balance constraint: sum(P_i) = P_demand + P_loss
        
        Parameters:
        power_output: array of power outputs
        
        Returns:
        constraint_violation: difference from power balance
        """
        total_generation = np.sum(power_output)
        transmission_loss = self.transmission_loss(power_output)
        return total_generation - self.demand - transmission_loss
    
    def solve_eld(self):
        """
        Solve the Economic Load Dispatch problem using optimization
        
        Returns:
        result: optimization result containing optimal power dispatch
        """
        # Initial guess: equal distribution
        x0 = np.full(self.n_gen, self.demand / self.n_gen)
        
        # Bounds for each generator
        bounds = [(gen['Pmin'], gen['Pmax']) for gen in self.generators]
        
        # Constraints
        constraints = [
            {'type': 'eq', 'fun': self.power_balance_constraint}
        ]
        
        # Solve optimization problem
        result = minimize(
            fun=self.cost_function,
            x0=x0,
            method='SLSQP',
            bounds=bounds,
            constraints=constraints,
            options={'disp': True, 'maxiter': 1000}
        )
        
        return result
    
    def calculate_lambda(self, power_output):
        """
        Calculate incremental cost (lambda) for each generator
        
        Parameters:
        power_output: optimal power outputs
        
        Returns:
        lambdas: incremental costs for each generator
        """
        lambdas = []
        for i, gen in enumerate(self.generators):
            a, b = gen['a'], gen['b']
            P = power_output[i]
            lambda_i = 2 * a * P + b
            lambdas.append(lambda_i)
        return np.array(lambdas)

# Define generator data for our example problem
generators_data = [
    {'name': 'Generator 1', 'a': 0.0080, 'b': 7.0, 'c': 200, 'Pmin': 50, 'Pmax': 200},
    {'name': 'Generator 2', 'a': 0.0090, 'b': 6.5, 'c': 180, 'Pmin': 30, 'Pmax': 150},
    {'name': 'Generator 3', 'a': 0.0070, 'b': 8.0, 'c': 220, 'Pmin': 40, 'Pmax': 180},
    {'name': 'Generator 4', 'a': 0.0085, 'b': 7.5, 'c': 190, 'Pmin': 20, 'Pmax': 120}
]

# Power demand
total_demand = 400  # MW

# Transmission loss coefficients (B-coefficients)
loss_coefficients = [0.0001, 0.00012, 0.00008, 0.00015]

# Create and solve ELD problem
print("=" * 60)
print("ECONOMIC LOAD DISPATCH OPTIMIZATION")
print("=" * 60)
print(f"Total Power Demand: {total_demand} MW")
print(f"Number of Generators: {len(generators_data)}")
print()

# Initialize ELD solver
eld_solver = EconomicLoadDispatch(generators_data, total_demand, loss_coefficients)

# Solve the optimization problem
print("Solving Economic Load Dispatch...")
result = eld_solver.solve_eld()

# Extract results
optimal_dispatch = result.x
optimal_cost = result.fun
transmission_loss = eld_solver.transmission_loss(optimal_dispatch)
incremental_costs = eld_solver.calculate_lambda(optimal_dispatch)

# Display results
print("\n" + "=" * 60)
print("OPTIMIZATION RESULTS")
print("=" * 60)

print(f"Optimization Status: {'Success' if result.success else 'Failed'}")
print(f"Total Generation Cost: ${optimal_cost:.2f}")
print(f"Transmission Loss: {transmission_loss:.2f} MW")
print(f"Total Generation: {np.sum(optimal_dispatch):.2f} MW")
print(f"Net Power to Load: {np.sum(optimal_dispatch) - transmission_loss:.2f} MW")
print()

# Create detailed results table
results_df = pd.DataFrame({
    'Generator': [gen['name'] for gen in generators_data],
    'Power Output (MW)': optimal_dispatch,
    'Individual Cost ($)': [eld_solver.cost_function([0]*i + [optimal_dispatch[i]] + [0]*(len(generators_data)-i-1)) 
                           - eld_solver.cost_function([0]*len(generators_data)) for i in range(len(generators_data))],
    'Incremental Cost ($/MWh)': incremental_costs,
    'Capacity Utilization (%)': [optimal_dispatch[i]/(generators_data[i]['Pmax']-generators_data[i]['Pmin'])*100 
                                for i in range(len(generators_data))]
})

print("DETAILED DISPATCH RESULTS:")
print("-" * 80)
print(results_df.to_string(index=False, float_format='%.2f'))

# Verification of optimality condition
print(f"\nOptimality Check:")
print(f"Average Incremental Cost: {np.mean(incremental_costs):.3f} $/MWh")
print(f"Standard Deviation: {np.std(incremental_costs):.6f} $/MWh")
print("(All incremental costs should be approximately equal for optimal dispatch)")

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

# 1. Power Dispatch Bar Chart
ax1 = plt.subplot(3, 3, 1)
generators_names = [gen['name'] for gen in generators_data]
colors = plt.cm.Set3(np.linspace(0, 1, len(generators_data)))
bars = plt.bar(generators_names, optimal_dispatch, color=colors, alpha=0.8, edgecolor='black')
plt.title('Optimal Power Dispatch', fontsize=14, fontweight='bold')
plt.ylabel('Power Output (MW)', fontsize=12)
plt.xticks(rotation=45)
plt.grid(axis='y', alpha=0.3)

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

# 2. Cost Functions Visualization
ax2 = plt.subplot(3, 3, 2)
P_range = np.linspace(0, 250, 300)
for i, gen in enumerate(generators_data):
    cost_curve = gen['a'] * P_range**2 + gen['b'] * P_range + gen['c']
    plt.plot(P_range, cost_curve, label=gen['name'], linewidth=2)
    # Mark optimal operating point
    optimal_cost_point = gen['a'] * optimal_dispatch[i]**2 + gen['b'] * optimal_dispatch[i] + gen['c']
    plt.plot(optimal_dispatch[i], optimal_cost_point, 'o', markersize=8, color='red')

plt.title('Generator Cost Functions', fontsize=14, fontweight='bold')
plt.xlabel('Power Output (MW)', fontsize=12)
plt.ylabel('Cost ($)', fontsize=12)
plt.legend()
plt.grid(True, alpha=0.3)

# 3. Incremental Cost Comparison
ax3 = plt.subplot(3, 3, 3)
plt.bar(generators_names, incremental_costs, color=colors, alpha=0.8, edgecolor='black')
plt.title('Incremental Costs at Optimal Dispatch', fontsize=14, fontweight='bold')
plt.ylabel('Incremental Cost ($/MWh)', fontsize=12)
plt.xticks(rotation=45)
plt.grid(axis='y', alpha=0.3)
# Add horizontal line for average
avg_lambda = np.mean(incremental_costs)
plt.axhline(y=avg_lambda, color='red', linestyle='--', linewidth=2, 
            label=f'Average: {avg_lambda:.2f}')
plt.legend()

# 4. Capacity Utilization
ax4 = plt.subplot(3, 3, 4)
max_capacities = [gen['Pmax'] for gen in generators_data]
utilization = (optimal_dispatch / max_capacities) * 100
bars = plt.bar(generators_names, utilization, color=colors, alpha=0.8, edgecolor='black')
plt.title('Generator Capacity Utilization', fontsize=14, fontweight='bold')
plt.ylabel('Utilization (%)', fontsize=12)
plt.xticks(rotation=45)
plt.grid(axis='y', alpha=0.3)
plt.ylim(0, 100)

# Add value labels
for bar, value in zip(bars, utilization):
    plt.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1, 
             f'{value:.1f}%', ha='center', va='bottom', fontweight='bold')

# 5. Cost Breakdown Pie Chart
ax5 = plt.subplot(3, 3, 5)
individual_costs = []
for i in range(len(generators_data)):
    gen = generators_data[i]
    cost = gen['a'] * optimal_dispatch[i]**2 + gen['b'] * optimal_dispatch[i] + gen['c']
    individual_costs.append(cost)

plt.pie(individual_costs, labels=generators_names, autopct='%1.1f%%', 
        colors=colors, startangle=90)
plt.title('Cost Distribution Among Generators', fontsize=14, fontweight='bold')

# 6. Sensitivity Analysis - Demand vs Total Cost
ax6 = plt.subplot(3, 3, 6)
demand_range = np.linspace(300, 500, 20)
total_costs = []

for demand_test in demand_range:
    eld_test = EconomicLoadDispatch(generators_data, demand_test, loss_coefficients)
    result_test = eld_test.solve_eld()
    if result_test.success:
        total_costs.append(result_test.fun)
    else:
        total_costs.append(np.nan)

plt.plot(demand_range, total_costs, 'b-', linewidth=2, marker='o', markersize=4)
plt.plot(total_demand, optimal_cost, 'ro', markersize=10, label='Current Operating Point')
plt.title('Sensitivity Analysis: Demand vs Total Cost', fontsize=14, fontweight='bold')
plt.xlabel('Power Demand (MW)', fontsize=12)
plt.ylabel('Total Cost ($)', fontsize=12)
plt.grid(True, alpha=0.3)
plt.legend()

# 7. Generator Efficiency Comparison
ax7 = plt.subplot(3, 3, 7)
# Calculate average cost per MW for each generator at optimal dispatch
avg_costs = []
for i, gen in enumerate(generators_data):
    if optimal_dispatch[i] > 0:
        cost = gen['a'] * optimal_dispatch[i]**2 + gen['b'] * optimal_dispatch[i] + gen['c']
        avg_cost = cost / optimal_dispatch[i]
        avg_costs.append(avg_cost)
    else:
        avg_costs.append(0)

bars = plt.bar(generators_names, avg_costs, color=colors, alpha=0.8, edgecolor='black')
plt.title('Average Cost per MW at Optimal Dispatch', fontsize=14, fontweight='bold')
plt.ylabel('Average Cost ($/MW)', fontsize=12)
plt.xticks(rotation=45)
plt.grid(axis='y', alpha=0.3)

# 8. Load Factor Analysis
ax8 = plt.subplot(3, 3, 8)
min_outputs = [gen['Pmin'] for gen in generators_data]
load_factors = ((optimal_dispatch - min_outputs) / (max_capacities - np.array(min_outputs))) * 100
bars = plt.bar(generators_names, load_factors, color=colors, alpha=0.8, edgecolor='black')
plt.title('Load Factor (Above Minimum)', fontsize=14, fontweight='bold')
plt.ylabel('Load Factor (%)', fontsize=12)
plt.xticks(rotation=45)
plt.grid(axis='y', alpha=0.3)

# 9. 3D Surface Plot - Cost vs Two Generators
ax9 = plt.subplot(3, 3, 9, projection='3d')
P1_range = np.linspace(generators_data[0]['Pmin'], generators_data[0]['Pmax'], 30)
P2_range = np.linspace(generators_data[1]['Pmin'], generators_data[1]['Pmax'], 30)
P1_mesh, P2_mesh = np.meshgrid(P1_range, P2_range)

# Calculate remaining power for other generators
remaining_power = total_demand - P1_mesh - P2_mesh
# Simplified cost calculation for visualization
cost_surface = (generators_data[0]['a'] * P1_mesh**2 + generators_data[0]['b'] * P1_mesh + generators_data[0]['c'] +
                generators_data[1]['a'] * P2_mesh**2 + generators_data[1]['b'] * P2_mesh + generators_data[1]['c'])

surface = ax9.plot_surface(P1_mesh, P2_mesh, cost_surface, alpha=0.7, cmap='viridis')
ax9.scatter([optimal_dispatch[0]], [optimal_dispatch[1]], 
           [generators_data[0]['a'] * optimal_dispatch[0]**2 + generators_data[0]['b'] * optimal_dispatch[0] + generators_data[0]['c'] +
            generators_data[1]['a'] * optimal_dispatch[1]**2 + generators_data[1]['b'] * optimal_dispatch[1] + generators_data[1]['c']], 
           color='red', s=100, label='Optimal Point')
ax9.set_xlabel(f'{generators_data[0]["name"]} (MW)')
ax9.set_ylabel(f'{generators_data[1]["name"]} (MW)')
ax9.set_zlabel('Combined Cost ($)')
ax9.set_title('Cost Surface for Two Generators', fontsize=12, fontweight='bold')

plt.tight_layout()
plt.show()

# Additional Analysis: Economic Merit Order
print("\n" + "=" * 60)
print("ECONOMIC MERIT ORDER ANALYSIS")
print("=" * 60)

# Calculate full load average costs
merit_order = []
for i, gen in enumerate(generators_data):
    full_load_cost = gen['a'] * gen['Pmax']**2 + gen['b'] * gen['Pmax'] + gen['c']
    avg_cost = full_load_cost / gen['Pmax']
    merit_order.append({
        'Generator': gen['name'],
        'Average Cost ($/MW)': avg_cost,
        'Min Cost ($/MWh)': gen['b'],  # Marginal cost at minimum load
        'Capacity (MW)': gen['Pmax']
    })

merit_df = pd.DataFrame(merit_order)
merit_df = merit_df.sort_values('Average Cost ($/MW)')
print("Merit Order (Sorted by Average Cost):")
print("-" * 60)
print(merit_df.to_string(index=False, float_format='%.3f'))

print(f"\nNote: The economic dispatch optimizes the total system cost,")
print(f"which may differ from simple merit order due to transmission losses")
print(f"and generator constraints.")

def cost_function(self, power_output):

def transmission_loss(self, power_output):

def power_balance_constraint(self, power_output):

def calculate_lambda(self, power_output):

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.stats import linregress
import pandas as pd
from sklearn.metrics import r2_score
import warnings
warnings.filterwarnings('ignore')

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

# Simulate experimental data
def generate_pk_data():
    """
    Generate synthetic pharmacokinetic data for a one-compartment model
    with first-order elimination following IV administration
    """
    # True parameters (what we want to estimate)
    true_C0 = 100.0  # mg/L - initial concentration
    true_ke = 0.1    # h^-1 - elimination rate constant
    dose = 500       # mg - administered dose
    
    # Time points (hours)
    time_points = np.array([0.5, 1, 2, 4, 6, 8, 12, 16, 24, 32, 48])
    
    # Calculate true concentrations
    true_concentrations = true_C0 * np.exp(-true_ke * time_points)
    
    # Add realistic measurement noise (CV = 10%)
    noise_cv = 0.10
    noise = np.random.normal(1, noise_cv, len(time_points))
    observed_concentrations = true_concentrations * noise
    
    # Ensure no negative concentrations
    observed_concentrations = np.maximum(observed_concentrations, 0.1)
    
    return time_points, observed_concentrations, true_C0, true_ke, dose

# Generate data
time, conc, true_C0, true_ke, dose = generate_pk_data()

print("Generated Pharmacokinetic Data:")
print("Time (h)\tConcentration (mg/L)")
for t, c in zip(time, conc):
    print(f"{t:4.1f}\t\t{c:6.2f}")

# Method 1: Linear regression on log-transformed data
def linear_regression_method(time, conc):
    """
    Estimate parameters using linear regression on log-transformed data
    ln(C) = ln(C0) - ke * t
    """
    # Transform concentrations to natural log
    ln_conc = np.log(conc)
    
    # Perform linear regression
    slope, intercept, r_value, p_value, std_err = linregress(time, ln_conc)
    
    # Extract parameters
    ke_est = -slope  # elimination rate constant
    C0_est = np.exp(intercept)  # initial concentration
    r_squared = r_value**2
    
    return ke_est, C0_est, r_squared

# Method 2: Non-linear least squares fitting
def pk_model(t, C0, ke):
    """One-compartment pharmacokinetic model"""
    return C0 * np.exp(-ke * t)

def nonlinear_fitting_method(time, conc):
    """
    Estimate parameters using non-linear least squares fitting
    """
    # Initial parameter guesses
    initial_guess = [conc[0], 0.1]  # C0 = first concentration, ke = 0.1
    
    # Perform curve fitting
    popt, pcov = curve_fit(pk_model, time, conc, p0=initial_guess)
    
    # Extract parameters and calculate confidence intervals
    C0_est, ke_est = popt
    param_errors = np.sqrt(np.diag(pcov))
    C0_error, ke_error = param_errors
    
    # Calculate R-squared
    y_pred = pk_model(time, C0_est, ke_est)
    r_squared = r2_score(conc, y_pred)
    
    return C0_est, ke_est, C0_error, ke_error, r_squared

# Apply both methods
print("\n" + "="*60)
print("PARAMETER ESTIMATION RESULTS")
print("="*60)

# Method 1: Linear regression
ke_lr, C0_lr, r2_lr = linear_regression_method(time, conc)
print(f"\nMethod 1: Linear Regression on Log-Transformed Data")
print(f"Elimination rate constant (ke): {ke_lr:.4f} h⁻¹")
print(f"Initial concentration (C0): {C0_lr:.2f} mg/L")
print(f"R-squared: {r2_lr:.4f}")

# Method 2: Non-linear fitting
C0_nlf, ke_nlf, C0_err, ke_err, r2_nlf = nonlinear_fitting_method(time, conc)
print(f"\nMethod 2: Non-linear Least Squares Fitting")
print(f"Initial concentration (C0): {C0_nlf:.2f} ± {C0_err:.2f} mg/L")
print(f"Elimination rate constant (ke): {ke_nlf:.4f} ± {ke_err:.4f} h⁻¹")
print(f"R-squared: {r2_nlf:.4f}")

# Calculate derived pharmacokinetic parameters
def calculate_pk_parameters(C0, ke, dose):
    """
    Calculate derived pharmacokinetic parameters
    """
    # Half-life
    t_half = np.log(2) / ke
    
    # Volume of distribution
    Vd = dose / C0
    
    # Clearance
    CL = ke * Vd
    
    # Area under the curve (AUC) from 0 to infinity
    AUC_inf = C0 / ke
    
    return t_half, Vd, CL, AUC_inf

# Calculate parameters using non-linear fitting results
t_half, Vd, CL, AUC_inf = calculate_pk_parameters(C0_nlf, ke_nlf, dose)

print(f"\nDerived Pharmacokinetic Parameters:")
print(f"Half-life (t₁/₂): {t_half:.2f} hours")
print(f"Volume of distribution (Vd): {Vd:.2f} L")
print(f"Clearance (CL): {CL:.2f} L/h")
print(f"AUC₀₋∞: {AUC_inf:.2f} mg·h/L")

# Compare with true values
print(f"\nComparison with True Values:")
print(f"True C0: {true_C0:.2f} mg/L, Estimated: {C0_nlf:.2f} mg/L")
print(f"True ke: {true_ke:.4f} h⁻¹, Estimated: {ke_nlf:.4f} h⁻¹")

# Create comprehensive visualization
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))

# Plot 1: Semi-log plot with both fitting methods
time_smooth = np.linspace(0, 50, 1000)
conc_lr = pk_model(time_smooth, C0_lr, ke_lr)
conc_nlf = pk_model(time_smooth, C0_nlf, ke_nlf)
conc_true = pk_model(time_smooth, true_C0, true_ke)

ax1.semilogy(time, conc, 'ro', markersize=8, label='Observed Data', alpha=0.8)
ax1.semilogy(time_smooth, conc_true, 'k--', linewidth=2, label=f'True Model (ke={true_ke:.3f})')
ax1.semilogy(time_smooth, conc_lr, 'b-', linewidth=2, label=f'Linear Regression (ke={ke_lr:.3f})')
ax1.semilogy(time_smooth, conc_nlf, 'g-', linewidth=2, label=f'Non-linear Fitting (ke={ke_nlf:.3f})')
ax1.set_xlabel('Time (hours)')
ax1.set_ylabel('Concentration (mg/L)')
ax1.set_title('Pharmacokinetic Model Fitting\n(Semi-logarithmic Scale)')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2: Linear scale
ax2.plot(time, conc, 'ro', markersize=8, label='Observed Data', alpha=0.8)
ax2.plot(time_smooth, conc_true, 'k--', linewidth=2, label='True Model')
ax2.plot(time_smooth, conc_lr, 'b-', linewidth=2, label='Linear Regression')
ax2.plot(time_smooth, conc_nlf, 'g-', linewidth=2, label='Non-linear Fitting')
ax2.set_xlabel('Time (hours)')
ax2.set_ylabel('Concentration (mg/L)')
ax2.set_title('Pharmacokinetic Model Fitting\n(Linear Scale)')
ax2.legend()
ax2.grid(True, alpha=0.3)

# Plot 3: Residuals analysis
residuals_lr = conc - pk_model(time, C0_lr, ke_lr)
residuals_nlf = conc - pk_model(time, C0_nlf, ke_nlf)

ax3.scatter(time, residuals_lr, color='blue', alpha=0.7, s=60, label='Linear Regression')
ax3.scatter(time, residuals_nlf, color='green', alpha=0.7, s=60, label='Non-linear Fitting')
ax3.axhline(y=0, color='red', linestyle='--', alpha=0.8)
ax3.set_xlabel('Time (hours)')
ax3.set_ylabel('Residuals (mg/L)')
ax3.set_title('Residuals Analysis')
ax3.legend()
ax3.grid(True, alpha=0.3)

# Plot 4: Parameter comparison
methods = ['Linear\nRegression', 'Non-linear\nFitting', 'True\nValue']
ke_values = [ke_lr, ke_nlf, true_ke]
C0_values = [C0_lr, C0_nlf, true_C0]

x_pos = np.arange(len(methods))
width = 0.35

bars1 = ax4.bar(x_pos - width/2, ke_values, width, label='ke (h⁻¹)', alpha=0.8, color='skyblue')
ax4_twin = ax4.twinx()
bars2 = ax4_twin.bar(x_pos + width/2, C0_values, width, label='C₀ (mg/L)', alpha=0.8, color='lightcoral')

ax4.set_xlabel('Estimation Method')
ax4.set_ylabel('Elimination Rate Constant (h⁻¹)', color='blue')
ax4_twin.set_ylabel('Initial Concentration (mg/L)', color='red')
ax4.set_title('Parameter Estimation Comparison')
ax4.set_xticks(x_pos)
ax4.set_xticklabels(methods)
ax4.legend(loc='upper left')
ax4_twin.legend(loc='upper right')

# Add value labels on bars
for i, (ke, C0) in enumerate(zip(ke_values, C0_values)):
    ax4.text(i - width/2, ke + 0.005, f'{ke:.3f}', ha='center', va='bottom', fontsize=9)
    ax4_twin.text(i + width/2, C0 + 2, f'{C0:.1f}', ha='center', va='bottom', fontsize=9)

plt.tight_layout()
plt.show()

# Summary statistics table
print(f"\n" + "="*80)
print("SUMMARY OF PHARMACOKINETIC ANALYSIS")
print("="*80)

summary_data = {
    'Parameter': ['C₀ (mg/L)', 'kₑ (h⁻¹)', 't₁/₂ (h)', 'Vd (L)', 'CL (L/h)', 'AUC₀₋∞ (mg·h/L)', 'R²'],
    'True Value': [true_C0, true_ke, np.log(2)/true_ke, dose/true_C0, 
                   true_ke * (dose/true_C0), true_C0/true_ke, 1.000],
    'Linear Regression': [C0_lr, ke_lr, np.log(2)/ke_lr, dose/C0_lr, 
                         ke_lr * (dose/C0_lr), C0_lr/ke_lr, r2_lr],
    'Non-linear Fitting': [C0_nlf, ke_nlf, np.log(2)/ke_nlf, dose/C0_nlf, 
                          ke_nlf * (dose/C0_nlf), C0_nlf/ke_nlf, r2_nlf]
}

summary_df = pd.DataFrame(summary_data)
print(summary_df.round(3))

print(f"\n" + "="*80)
print("INTERPRETATION AND CONCLUSIONS")
print("="*80)
print(f"1. Both methods provided good estimates of the pharmacokinetic parameters")
print(f"2. Non-linear fitting generally provides more accurate parameter estimates")
print(f"3. The drug shows typical first-order elimination kinetics")
print(f"4. Half-life of ~{t_half:.1f} hours suggests {2-3 if t_half < 8 else 1-2} times daily dosing")
print(f"5. Clearance of {CL:.1f} L/h indicates moderate hepatic/renal elimination")