import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from matplotlib.patches import Circle
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

# Problem parameters
START = np.array([0.0, 0.0])
GOAL = np.array([10.0, 10.0])
T = 5.0  # Total time (seconds)
N = 50   # Number of time steps
dt = T / N

# Obstacles: [x, y, radius]
OBSTACLES = [
    [3.0, 3.0, 1.0],
    [7.0, 5.0, 1.2],
    [5.0, 8.0, 0.8]
]

# Physical limits
V_MAX = 5.0   # Maximum velocity
A_MAX = 3.0   # Maximum acceleration

# Cost function weights
ALPHA = 1.0   # Velocity cost weight
BETA = 10.0   # Acceleration cost weight
GAMMA = 100.0 # Obstacle penalty weight
DELTA = 50.0  # Constraint violation penalty

def trajectory_from_params(params):
    """
    Convert optimization parameters to trajectory.
    params: flattened array of [x1, y1, x2, y2, ..., x_{N-1}, y_{N-1}]
    Returns: (N+1) x 2 array including start and goal
    """
    waypoints = params.reshape(-1, 2)
    trajectory = np.vstack([START, waypoints, GOAL])
    return trajectory

def compute_velocity(trajectory):
    """Compute velocity using finite differences"""
    return np.diff(trajectory, axis=0) / dt

def compute_acceleration(velocity):
    """Compute acceleration using finite differences"""
    return np.diff(velocity, axis=0) / dt

def objective_function(params):
    """
    Compute total cost: energy + penalties for constraints
    """
    trajectory = trajectory_from_params(params)
    
    # Compute velocity and acceleration
    velocity = compute_velocity(trajectory)
    acceleration = compute_acceleration(velocity)
    
    # Energy cost (velocity and acceleration terms)
    velocity_cost = ALPHA * np.sum(np.linalg.norm(velocity, axis=1)**2) * dt
    acceleration_cost = BETA * np.sum(np.linalg.norm(acceleration, axis=1)**2) * dt
    
    # Obstacle avoidance penalty (soft constraint)
    obstacle_penalty = 0.0
    for obs in OBSTACLES:
        obs_center = np.array(obs[:2])
        obs_radius = obs[2]
        # Check all waypoints
        for point in trajectory:
            distance = np.linalg.norm(point - obs_center)
            safety_margin = 0.3  # Additional safety buffer
            if distance < obs_radius + safety_margin:
                # Quadratic penalty that increases as we get closer
                violation = (obs_radius + safety_margin - distance)
                obstacle_penalty += GAMMA * violation**2
    
    # Velocity constraint penalty
    velocity_penalty = 0.0
    for v in velocity:
        v_norm = np.linalg.norm(v)
        if v_norm > V_MAX:
            velocity_penalty += DELTA * (v_norm - V_MAX)**2
    
    # Acceleration constraint penalty
    acceleration_penalty = 0.0
    for a in acceleration:
        a_norm = np.linalg.norm(a)
        if a_norm > A_MAX:
            acceleration_penalty += DELTA * (a_norm - A_MAX)**2
    
    total_cost = (velocity_cost + acceleration_cost + 
                  obstacle_penalty + velocity_penalty + acceleration_penalty)
    
    return total_cost

def optimize_trajectory():
    """
    Optimize the trajectory using scipy's minimize
    """
    # Initial guess: linear interpolation between start and goal
    initial_waypoints = np.linspace(START, GOAL, N+1)[1:-1]
    initial_params = initial_waypoints.flatten()
    
    print("Starting optimization...")
    print(f"Initial cost: {objective_function(initial_params):.2f}")
    
    # Optimize using L-BFGS-B method (handles bounds well)
    result = minimize(
        objective_function,
        initial_params,
        method='L-BFGS-B',
        options={'maxiter': 500, 'disp': True}
    )
    
    print(f"\nOptimization completed!")
    print(f"Final cost: {result.fun:.2f}")
    print(f"Success: {result.success}")
    
    return trajectory_from_params(result.x), result

def analyze_trajectory(trajectory):
    """
    Compute and display trajectory statistics
    """
    velocity = compute_velocity(trajectory)
    acceleration = compute_acceleration(velocity)
    
    v_norms = np.linalg.norm(velocity, axis=1)
    a_norms = np.linalg.norm(acceleration, axis=1)
    
    print("\n" + "="*60)
    print("TRAJECTORY ANALYSIS")
    print("="*60)
    print(f"Total path length: {np.sum(v_norms) * dt:.3f} meters")
    print(f"Max velocity: {np.max(v_norms):.3f} m/s (limit: {V_MAX} m/s)")
    print(f"Max acceleration: {np.max(a_norms):.3f} m/s² (limit: {A_MAX} m/s²)")
    print(f"Average velocity: {np.mean(v_norms):.3f} m/s")
    print(f"Average acceleration: {np.mean(a_norms):.3f} m/s²")
    
    # Check obstacle clearance
    print("\nObstacle clearance:")
    for i, obs in enumerate(OBSTACLES):
        obs_center = np.array(obs[:2])
        obs_radius = obs[2]
        min_distance = np.min([np.linalg.norm(point - obs_center) 
                               for point in trajectory])
        clearance = min_distance - obs_radius
        print(f"  Obstacle {i+1}: {clearance:.3f} meters")
    
    return velocity, acceleration, v_norms, a_norms

def plot_results(trajectory, velocity, acceleration, v_norms, a_norms):
    """
    Create comprehensive visualization of the optimized trajectory
    """
    time_points = np.linspace(0, T, len(trajectory))
    time_v = np.linspace(0, T, len(velocity))
    time_a = np.linspace(0, T, len(acceleration))
    
    fig = plt.figure(figsize=(16, 12))
    
    # 1. Trajectory in workspace
    ax1 = plt.subplot(2, 3, 1)
    ax1.plot(trajectory[:, 0], trajectory[:, 1], 'b-', linewidth=2, label='Optimized Path')
    ax1.plot(START[0], START[1], 'go', markersize=12, label='Start')
    ax1.plot(GOAL[0], GOAL[1], 'r*', markersize=15, label='Goal')
    
    # Draw obstacles
    for i, obs in enumerate(OBSTACLES):
        circle = Circle((obs[0], obs[1]), obs[2], color='red', alpha=0.3)
        ax1.add_patch(circle)
        ax1.plot(obs[0], obs[1], 'rx', markersize=10)
    
    ax1.set_xlabel('X position (m)', fontsize=11)
    ax1.set_ylabel('Y position (m)', fontsize=11)
    ax1.set_title('Optimized Trajectory in Workspace', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.axis('equal')
    
    # 2. X and Y positions over time
    ax2 = plt.subplot(2, 3, 2)
    ax2.plot(time_points, trajectory[:, 0], 'b-', linewidth=2, label='X position')
    ax2.plot(time_points, trajectory[:, 1], 'r-', linewidth=2, label='Y position')
    ax2.set_xlabel('Time (s)', fontsize=11)
    ax2.set_ylabel('Position (m)', fontsize=11)
    ax2.set_title('Position vs Time', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 3. Velocity magnitude over time
    ax3 = plt.subplot(2, 3, 3)
    ax3.plot(time_v, v_norms, 'g-', linewidth=2, label='Velocity magnitude')
    ax3.axhline(y=V_MAX, color='r', linestyle='--', linewidth=2, label=f'Limit ({V_MAX} m/s)')
    ax3.set_xlabel('Time (s)', fontsize=11)
    ax3.set_ylabel('Velocity (m/s)', fontsize=11)
    ax3.set_title('Velocity Profile', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 4. Acceleration magnitude over time
    ax4 = plt.subplot(2, 3, 4)
    ax4.plot(time_a, a_norms, 'm-', linewidth=2, label='Acceleration magnitude')
    ax4.axhline(y=A_MAX, color='r', linestyle='--', linewidth=2, label=f'Limit ({A_MAX} m/s²)')
    ax4.set_xlabel('Time (s)', fontsize=11)
    ax4.set_ylabel('Acceleration (m/s²)', fontsize=11)
    ax4.set_title('Acceleration Profile', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 5. Velocity components
    ax5 = plt.subplot(2, 3, 5)
    ax5.plot(time_v, velocity[:, 0], 'b-', linewidth=2, label='X velocity')
    ax5.plot(time_v, velocity[:, 1], 'r-', linewidth=2, label='Y velocity')
    ax5.set_xlabel('Time (s)', fontsize=11)
    ax5.set_ylabel('Velocity (m/s)', fontsize=11)
    ax5.set_title('Velocity Components', fontsize=12, fontweight='bold')
    ax5.legend()
    ax5.grid(True, alpha=0.3)
    
    # 6. Energy consumption over time
    ax6 = plt.subplot(2, 3, 6)
    kinetic_energy = 0.5 * v_norms**2  # Assuming unit mass
    cumulative_energy = np.cumsum(kinetic_energy) * dt
    ax6.plot(time_v, cumulative_energy, 'orange', linewidth=2)
    ax6.set_xlabel('Time (s)', fontsize=11)
    ax6.set_ylabel('Cumulative Energy (J)', fontsize=11)
    ax6.set_title('Energy Consumption', fontsize=12, fontweight='bold')
    ax6.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('trajectory_optimization.png', dpi=150, bbox_inches='tight')
    plt.show()
    
    print("\n" + "="*60)
    print("Visualization complete! Graph saved as 'trajectory_optimization.png'")
    print("="*60)

# Main execution
print("="*60)
print("ROBOT TRAJECTORY OPTIMIZATION")
print("Minimizing Energy While Ensuring Safety")
print("="*60)
print(f"\nProblem Setup:")
print(f"  Start position: {START}")
print(f"  Goal position: {GOAL}")
print(f"  Time horizon: {T} seconds")
print(f"  Number of waypoints: {N-1}")
print(f"  Number of obstacles: {len(OBSTACLES)}")
print(f"  Max velocity: {V_MAX} m/s")
print(f"  Max acceleration: {A_MAX} m/s²")
print("\n")

# Run optimization
optimized_trajectory, result = optimize_trajectory()

# Analyze results
velocity, acceleration, v_norms, a_norms = analyze_trajectory(optimized_trajectory)

# Visualize
plot_results(optimized_trajectory, velocity, acceleration, v_norms, a_norms)

print("\n✓ Optimization complete! Check the graphs above for detailed results.")

START = np.array([0.0, 0.0])
GOAL = np.array([10.0, 10.0])
T = 5.0  # Total time
N = 50   # Number of discretization points

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_val_score, learning_curve
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error, r2_score
import seaborn as sns
from scipy import stats

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

# ============================================================================
# 1. Generate Synthetic Gene Expression Data
# ============================================================================

def generate_gene_expression_data(n_samples=200, n_genes=15, noise_level=0.3):
    """
    Generate synthetic gene expression data with realistic characteristics.
    
    Parameters:
    -----------
    n_samples : int
        Number of samples (e.g., different conditions or time points)
    n_genes : int
        Number of regulator genes
    noise_level : float
        Standard deviation of Gaussian noise added to measurements
    
    Returns:
    --------
    X : ndarray, shape (n_samples, n_genes)
        Gene expression levels for regulator genes
    y : ndarray, shape (n_samples,)
        Target gene expression levels
    true_coefficients : ndarray
        True coefficients used to generate the target
    """
    
    # Generate regulator gene expression levels
    # Use correlated features to mimic real biological networks
    mean = np.zeros(n_genes)
    
    # Create correlation structure
    correlation = 0.3
    cov = np.eye(n_genes) * (1 - correlation) + np.ones((n_genes, n_genes)) * correlation
    
    X = np.random.multivariate_normal(mean, cov, n_samples)
    
    # Create true coefficients with sparse structure (only some genes matter)
    true_coefficients = np.zeros(n_genes)
    # First 5 genes have strong effects
    true_coefficients[0:5] = [2.5, -1.8, 3.2, -2.1, 1.5]
    # Next 3 genes have moderate effects
    true_coefficients[5:8] = [0.8, -0.6, 0.9]
    # Remaining genes have no effect (noise)
    
    # Generate target gene expression
    y_true = X @ true_coefficients
    
    # Add measurement noise
    noise = np.random.normal(0, noise_level, n_samples)
    y = y_true + noise
    
    return X, y, true_coefficients

# Generate data
X, y, true_coefs = generate_gene_expression_data(n_samples=200, n_genes=15, noise_level=0.5)

# Split into training and test sets
split_idx = 150
X_train, X_test = X[:split_idx], X[split_idx:]
y_train, y_test = y[:split_idx], y[split_idx:]

# Standardize features
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

print("Dataset Information:")
print(f"Training samples: {X_train.shape[0]}")
print(f"Test samples: {X_test.shape[0]}")
print(f"Number of genes: {X_train.shape[1]}")
print(f"Target expression range: [{y.min():.2f}, {y.max():.2f}]")
print("\n" + "="*70 + "\n")

# ============================================================================
# 2. Hyperparameter Optimization via Cross-Validation
# ============================================================================

# Define range of alpha values to test (logarithmic scale)
alphas = np.logspace(-3, 3, 50)

# Store cross-validation scores
cv_scores_mean = []
cv_scores_std = []

print("Performing hyperparameter optimization...")
print("Testing {} alpha values from {:.4f} to {:.2f}".format(len(alphas), alphas[0], alphas[-1]))

for alpha in alphas:
    model = Ridge(alpha=alpha)
    # 5-fold cross-validation
    scores = cross_val_score(model, X_train_scaled, y_train, 
                            cv=5, scoring='neg_mean_squared_error')
    # Convert to RMSE
    rmse_scores = np.sqrt(-scores)
    cv_scores_mean.append(rmse_scores.mean())
    cv_scores_std.append(rmse_scores.std())

cv_scores_mean = np.array(cv_scores_mean)
cv_scores_std = np.array(cv_scores_std)

# Find optimal alpha
optimal_idx = np.argmin(cv_scores_mean)
optimal_alpha = alphas[optimal_idx]
optimal_cv_score = cv_scores_mean[optimal_idx]

print(f"\nOptimal alpha: {optimal_alpha:.4f}")
print(f"CV RMSE at optimal alpha: {optimal_cv_score:.4f}")
print("\n" + "="*70 + "\n")

# ============================================================================
# 3. Train Models with Different Hyperparameters
# ============================================================================

# Compare three models: under-regularized, optimal, over-regularized
alpha_underreg = alphas[max(0, optimal_idx - 15)]
alpha_optimal = optimal_alpha
alpha_overreg = alphas[min(len(alphas)-1, optimal_idx + 15)]

models = {
    'Under-regularized': Ridge(alpha=alpha_underreg),
    'Optimal': Ridge(alpha=alpha_optimal),
    'Over-regularized': Ridge(alpha=alpha_overreg)
}

results = {}
for name, model in models.items():
    model.fit(X_train_scaled, y_train)
    
    # Predictions
    y_train_pred = model.predict(X_train_scaled)
    y_test_pred = model.predict(X_test_scaled)
    
    # Metrics
    train_rmse = np.sqrt(mean_squared_error(y_train, y_train_pred))
    test_rmse = np.sqrt(mean_squared_error(y_test, y_test_pred))
    train_r2 = r2_score(y_train, y_train_pred)
    test_r2 = r2_score(y_test, y_test_pred)
    
    results[name] = {
        'model': model,
        'alpha': model.alpha,
        'train_rmse': train_rmse,
        'test_rmse': test_rmse,
        'train_r2': train_r2,
        'test_r2': test_r2,
        'coefficients': model.coef_,
        'y_test_pred': y_test_pred
    }
    
    print(f"{name} (alpha={model.alpha:.4f}):")
    print(f"  Training RMSE: {train_rmse:.4f}")
    print(f"  Test RMSE:     {test_rmse:.4f}")
    print(f"  Training R²:   {train_r2:.4f}")
    print(f"  Test R²:       {test_r2:.4f}")
    print(f"  Overfitting gap: {abs(train_rmse - test_rmse):.4f}")
    print()

# ============================================================================
# 4. Comprehensive Visualization
# ============================================================================

fig = plt.figure(figsize=(18, 12))
gs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3)

# Plot 1: Cross-validation curve
ax1 = fig.add_subplot(gs[0, :])
ax1.semilogx(alphas, cv_scores_mean, 'b-', linewidth=2, label='Mean CV RMSE')
ax1.fill_between(alphas, 
                 cv_scores_mean - cv_scores_std,
                 cv_scores_mean + cv_scores_std,
                 alpha=0.2, color='b', label='±1 std dev')
ax1.axvline(optimal_alpha, color='r', linestyle='--', linewidth=2, 
           label=f'Optimal α = {optimal_alpha:.4f}')
ax1.axvline(alpha_underreg, color='orange', linestyle=':', linewidth=1.5, alpha=0.7)
ax1.axvline(alpha_overreg, color='purple', linestyle=':', linewidth=1.5, alpha=0.7)
ax1.set_xlabel('Regularization Parameter (α)', fontsize=12, fontweight='bold')
ax1.set_ylabel('Cross-Validation RMSE', fontsize=12, fontweight='bold')
ax1.set_title('Hyperparameter Optimization: Cross-Validation Curve', 
             fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# Plot 2: Model comparison - Training vs Test RMSE
ax2 = fig.add_subplot(gs[1, 0])
model_names = list(results.keys())
train_rmses = [results[name]['train_rmse'] for name in model_names]
test_rmses = [results[name]['test_rmse'] for name in model_names]

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

bars1 = ax2.bar(x_pos - width/2, train_rmses, width, label='Training', 
               color='steelblue', alpha=0.8)
bars2 = ax2.bar(x_pos + width/2, test_rmses, width, label='Test', 
               color='coral', alpha=0.8)

ax2.set_xlabel('Model Configuration', fontsize=11, fontweight='bold')
ax2.set_ylabel('RMSE', fontsize=11, fontweight='bold')
ax2.set_title('Training vs Test Error', fontsize=12, fontweight='bold')
ax2.set_xticks(x_pos)
ax2.set_xticklabels(model_names, rotation=15, ha='right')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3, axis='y')

# Add value labels on bars
for bars in [bars1, bars2]:
    for bar in bars:
        height = bar.get_height()
        ax2.text(bar.get_x() + bar.get_width()/2., height,
                f'{height:.3f}', ha='center', va='bottom', fontsize=8)

# Plot 3: R² scores comparison
ax3 = fig.add_subplot(gs[1, 1])
train_r2s = [results[name]['train_r2'] for name in model_names]
test_r2s = [results[name]['test_r2'] for name in model_names]

bars1 = ax3.bar(x_pos - width/2, train_r2s, width, label='Training', 
               color='green', alpha=0.7)
bars2 = ax3.bar(x_pos + width/2, test_r2s, width, label='Test', 
               color='red', alpha=0.7)

ax3.set_xlabel('Model Configuration', fontsize=11, fontweight='bold')
ax3.set_ylabel('R² Score', fontsize=11, fontweight='bold')
ax3.set_title('Coefficient of Determination (R²)', fontsize=12, fontweight='bold')
ax3.set_xticks(x_pos)
ax3.set_xticklabels(model_names, rotation=15, ha='right')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3, axis='y')
ax3.set_ylim([0, 1.0])

# Add value labels
for bars in [bars1, bars2]:
    for bar in bars:
        height = bar.get_height()
        ax3.text(bar.get_x() + bar.get_width()/2., height,
                f'{height:.3f}', ha='center', va='bottom', fontsize=8)

# Plot 4: Overfitting analysis
ax4 = fig.add_subplot(gs[1, 2])
overfitting_gaps = [abs(results[name]['train_rmse'] - results[name]['test_rmse']) 
                    for name in model_names]
colors = ['orange', 'green', 'purple']
bars = ax4.bar(model_names, overfitting_gaps, color=colors, alpha=0.7)

ax4.set_xlabel('Model Configuration', fontsize=11, fontweight='bold')
ax4.set_ylabel('|Train RMSE - Test RMSE|', fontsize=11, fontweight='bold')
ax4.set_title('Overfitting Gap Analysis', fontsize=12, fontweight='bold')
ax4.set_xticklabels(model_names, rotation=15, ha='right')
ax4.grid(True, alpha=0.3, axis='y')

for bar in bars:
    height = bar.get_height()
    ax4.text(bar.get_x() + bar.get_width()/2., height,
            f'{height:.3f}', ha='center', va='bottom', fontsize=9)

# Plot 5: Coefficient comparison
ax5 = fig.add_subplot(gs[2, :])
gene_indices = np.arange(len(true_coefs))
width = 0.2

bars1 = ax5.bar(gene_indices - width*1.5, true_coefs, width, 
               label='True Coefficients', color='black', alpha=0.8)
bars2 = ax5.bar(gene_indices - width*0.5, results['Under-regularized']['coefficients'], 
               width, label='Under-regularized', color='orange', alpha=0.7)
bars3 = ax5.bar(gene_indices + width*0.5, results['Optimal']['coefficients'], 
               width, label='Optimal', color='green', alpha=0.7)
bars4 = ax5.bar(gene_indices + width*1.5, results['Over-regularized']['coefficients'], 
               width, label='Over-regularized', color='purple', alpha=0.7)

ax5.set_xlabel('Gene Index', fontsize=11, fontweight='bold')
ax5.set_ylabel('Coefficient Value', fontsize=11, fontweight='bold')
ax5.set_title('Learned Coefficients vs True Coefficients', fontsize=12, fontweight='bold')
ax5.set_xticks(gene_indices)
ax5.legend(fontsize=9, loc='upper right')
ax5.grid(True, alpha=0.3, axis='y')
ax5.axhline(y=0, color='k', linestyle='-', linewidth=0.5)

plt.suptitle('Gene Expression Model: Hyperparameter Optimization Analysis', 
            fontsize=16, fontweight='bold', y=0.995)

plt.tight_layout()
plt.show()

# ============================================================================
# 5. Prediction Quality Visualization
# ============================================================================

fig2, axes = plt.subplots(1, 3, figsize=(18, 5))

for idx, (name, result) in enumerate(results.items()):
    ax = axes[idx]
    
    # Scatter plot of predicted vs actual
    ax.scatter(y_test, result['y_test_pred'], alpha=0.6, s=80, 
              color=colors[idx], edgecolors='black', linewidth=0.5)
    
    # Perfect prediction line
    min_val = min(y_test.min(), result['y_test_pred'].min())
    max_val = max(y_test.max(), result['y_test_pred'].max())
    ax.plot([min_val, max_val], [min_val, max_val], 'r--', 
           linewidth=2, label='Perfect Prediction')
    
    # Regression line
    slope, intercept, r_value, p_value, std_err = stats.linregress(y_test, result['y_test_pred'])
    line_x = np.array([min_val, max_val])
    line_y = slope * line_x + intercept
    ax.plot(line_x, line_y, 'b-', linewidth=2, alpha=0.7, label='Fitted Line')
    
    ax.set_xlabel('Actual Expression', fontsize=11, fontweight='bold')
    ax.set_ylabel('Predicted Expression', fontsize=11, fontweight='bold')
    ax.set_title(f'{name}\nR² = {result["test_r2"]:.3f}, RMSE = {result["test_rmse"]:.3f}', 
                fontsize=12, fontweight='bold')
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)
    
    # Add diagonal reference
    ax.set_aspect('equal', adjustable='box')

plt.suptitle('Prediction Quality: Actual vs Predicted Gene Expression', 
            fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# ============================================================================
# 6. Learning Curves
# ============================================================================

fig3, axes = plt.subplots(1, 3, figsize=(18, 5))

for idx, (name, result) in enumerate(results.items()):
    ax = axes[idx]
    
    train_sizes, train_scores, val_scores = learning_curve(
        Ridge(alpha=result['alpha']), X_train_scaled, y_train,
        cv=5, scoring='neg_mean_squared_error',
        train_sizes=np.linspace(0.1, 1.0, 10), n_jobs=-1
    )
    
    train_rmse_mean = np.sqrt(-train_scores.mean(axis=1))
    train_rmse_std = np.sqrt(-train_scores).std(axis=1)
    val_rmse_mean = np.sqrt(-val_scores.mean(axis=1))
    val_rmse_std = np.sqrt(-val_scores).std(axis=1)
    
    ax.plot(train_sizes, train_rmse_mean, 'o-', color='blue', 
           linewidth=2, markersize=8, label='Training RMSE')
    ax.fill_between(train_sizes, 
                    train_rmse_mean - train_rmse_std,
                    train_rmse_mean + train_rmse_std,
                    alpha=0.2, color='blue')
    
    ax.plot(train_sizes, val_rmse_mean, 'o-', color='red', 
           linewidth=2, markersize=8, label='Validation RMSE')
    ax.fill_between(train_sizes, 
                    val_rmse_mean - val_rmse_std,
                    val_rmse_mean + val_rmse_std,
                    alpha=0.2, color='red')
    
    ax.set_xlabel('Training Set Size', fontsize=11, fontweight='bold')
    ax.set_ylabel('RMSE', fontsize=11, fontweight='bold')
    ax.set_title(f'{name} (α = {result["alpha"]:.4f})', 
                fontsize=12, fontweight='bold')
    ax.legend(fontsize=9, loc='best')
    ax.grid(True, alpha=0.3)

plt.suptitle('Learning Curves: Model Performance vs Training Set Size', 
            fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

print("="*70)
print("Analysis Complete!")
print("="*70)

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

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

# Define objective functions
def binding_affinity(x):
    """
    Target binding affinity (lower is better - stronger binding)
    x[0]: Molecular weight proxy (0-10)
    x[1]: Hydrophobicity index (0-10)
    x[2]: H-bond donors/acceptors (0-10)
    """
    return -10 * np.exp(-0.1 * (x[0] - 5)**2) - 8 * np.exp(-0.05 * (x[1] - 3)**2) + 0.5 * x[2]

def off_target_binding(x):
    """
    Off-target binding (lower is better - fewer side effects)
    """
    return 5 * np.exp(-0.08 * (x[0] - 7)**2) + 3 * np.sqrt(x[1] + 0.1)

def solubility_penalty(x):
    """
    Solubility represented by LogP
    Ideal range: 2-3 (for oral bioavailability)
    Penalty increases as we deviate from ideal range
    """
    logP = 0.5 * x[0] + 0.3 * x[1] - 0.2 * x[2]
    target_logP = 2.5
    return np.abs(logP - target_logP)

def combined_objective(x, weights):
    """
    Weighted sum of all objectives
    weights: [w_affinity, w_offtarget, w_solubility]
    """
    affinity = binding_affinity(x)
    offtarget = off_target_binding(x)
    solubility = solubility_penalty(x)
    
    # Normalize objectives (approximate ranges)
    affinity_norm = (affinity + 18) / 10  # Range roughly -18 to -8
    offtarget_norm = offtarget / 10  # Range roughly 0 to 10
    solubility_norm = solubility / 3  # Range roughly 0 to 3
    
    return weights[0] * affinity_norm + weights[1] * offtarget_norm + weights[2] * solubility_norm

# Optimization using Differential Evolution
bounds = [(0, 10), (0, 10), (0, 10)]

# Different weight scenarios
scenarios = {
    'Balanced': [1.0, 1.0, 1.0],
    'Efficacy-focused': [2.0, 0.5, 0.5],
    'Safety-focused': [0.5, 2.0, 1.0],
    'Solubility-focused': [0.5, 0.5, 2.0]
}

results = {}
print("="*70)
print("DRUG MOLECULAR DESIGN OPTIMIZATION RESULTS")
print("="*70)

for scenario_name, weights in scenarios.items():
    result = differential_evolution(
        lambda x: combined_objective(x, weights),
        bounds,
        maxiter=1000,
        popsize=15,
        tol=1e-7,
        seed=42
    )
    
    x_opt = result.x
    results[scenario_name] = {
        'x': x_opt,
        'affinity': binding_affinity(x_opt),
        'offtarget': off_target_binding(x_opt),
        'solubility': solubility_penalty(x_opt),
        'logP': 0.5 * x_opt[0] + 0.3 * x_opt[1] - 0.2 * x_opt[2]
    }
    
    print(f"\n{scenario_name} Approach:")
    print(f"  Molecular descriptors: x1={x_opt[0]:.3f}, x2={x_opt[1]:.3f}, x3={x_opt[2]:.3f}")
    print(f"  Binding affinity: {results[scenario_name]['affinity']:.3f} (more negative = stronger)")
    print(f"  Off-target binding: {results[scenario_name]['offtarget']:.3f} (lower = fewer side effects)")
    print(f"  Solubility penalty: {results[scenario_name]['solubility']:.3f} (lower = better)")
    print(f"  LogP value: {results[scenario_name]['logP']:.3f} (ideal: 2-3)")

print("\n" + "="*70)

# Visualization
fig = plt.figure(figsize=(16, 12))

# 1. Radar chart comparing scenarios
ax1 = plt.subplot(2, 3, 1, projection='polar')
categories = ['Binding\nAffinity', 'Low\nOff-target', 'Solubility']
N = len(categories)
angles = [n / float(N) * 2 * np.pi for n in range(N)]
angles += angles[:1]

for scenario_name, data in results.items():
    # Normalize for radar chart (0-1, higher is better)
    values = [
        1 - (data['affinity'] + 18) / 10,  # Convert to 0-1, higher better
        1 - data['offtarget'] / 10,  # Lower is better
        1 - data['solubility'] / 3   # Lower is better
    ]
    values += values[:1]
    ax1.plot(angles, values, 'o-', linewidth=2, label=scenario_name)
    ax1.fill(angles, values, alpha=0.15)

ax1.set_xticks(angles[:-1])
ax1.set_xticklabels(categories, size=9)
ax1.set_ylim(0, 1)
ax1.set_title('Multi-Objective Performance Comparison', size=11, weight='bold', pad=20)
ax1.legend(loc='upper right', bbox_to_anchor=(1.3, 1.1), fontsize=8)
ax1.grid(True)

# 2. Bar chart - Binding affinity
ax2 = plt.subplot(2, 3, 2)
scenario_names = list(results.keys())
affinities = [results[s]['affinity'] for s in scenario_names]
colors = sns.color_palette("husl", len(scenario_names))
bars = ax2.bar(scenario_names, affinities, color=colors, alpha=0.7, edgecolor='black')
ax2.set_ylabel('Binding Energy', fontweight='bold')
ax2.set_title('Target Binding Affinity\n(More negative = Stronger binding)', fontweight='bold')
ax2.axhline(y=-15, color='green', linestyle='--', alpha=0.5, label='Strong binding threshold')
ax2.legend(fontsize=8)
ax2.grid(axis='y', alpha=0.3)
plt.setp(ax2.xaxis.get_majorticklabels(), rotation=45, ha='right')

# 3. Bar chart - Off-target binding
ax3 = plt.subplot(2, 3, 3)
offtargets = [results[s]['offtarget'] for s in scenario_names]
bars = ax3.bar(scenario_names, offtargets, color=colors, alpha=0.7, edgecolor='black')
ax3.set_ylabel('Off-target Binding', fontweight='bold')
ax3.set_title('Off-Target Binding\n(Lower = Fewer side effects)', fontweight='bold')
ax3.axhline(y=5, color='orange', linestyle='--', alpha=0.5, label='Acceptable threshold')
ax3.legend(fontsize=8)
ax3.grid(axis='y', alpha=0.3)
plt.setp(ax3.xaxis.get_majorticklabels(), rotation=45, ha='right')

# 4. LogP visualization with ideal range
ax4 = plt.subplot(2, 3, 4)
logP_values = [results[s]['logP'] for s in scenario_names]
bars = ax4.bar(scenario_names, logP_values, color=colors, alpha=0.7, edgecolor='black')
ax4.axhspan(2, 3, alpha=0.2, color='green', label='Ideal range (2-3)')
ax4.set_ylabel('LogP Value', fontweight='bold')
ax4.set_title('Lipophilicity (LogP)\n(Affects oral bioavailability)', fontweight='bold')
ax4.legend(fontsize=8)
ax4.grid(axis='y', alpha=0.3)
plt.setp(ax4.xaxis.get_majorticklabels(), rotation=45, ha='right')

# 5. Molecular descriptors comparison
ax5 = plt.subplot(2, 3, 5)
x_labels = ['MW proxy', 'Hydrophobicity', 'H-bond D/A']
x_pos = np.arange(len(x_labels))
width = 0.2

for i, scenario_name in enumerate(scenario_names):
    x_vals = results[scenario_name]['x']
    ax5.bar(x_pos + i*width, x_vals, width, label=scenario_name, alpha=0.7)

ax5.set_xlabel('Molecular Descriptors', fontweight='bold')
ax5.set_ylabel('Descriptor Value', fontweight='bold')
ax5.set_title('Optimized Molecular Descriptors', fontweight='bold')
ax5.set_xticks(x_pos + width * 1.5)
ax5.set_xticklabels(x_labels)
ax5.legend(fontsize=8)
ax5.grid(axis='y', alpha=0.3)

# 6. 3D scatter plot of molecular descriptor space
ax6 = fig.add_subplot(2, 3, 6, projection='3d')
for i, scenario_name in enumerate(scenario_names):
    x_vals = results[scenario_name]['x']
    ax6.scatter(x_vals[0], x_vals[1], x_vals[2], 
               c=[colors[i]], s=200, marker='o', 
               edgecolors='black', linewidth=2,
               label=scenario_name, alpha=0.8)

ax6.set_xlabel('MW proxy (x1)', fontweight='bold')
ax6.set_ylabel('Hydrophobicity (x2)', fontweight='bold')
ax6.set_zlabel('H-bond D/A (x3)', fontweight='bold')
ax6.set_title('Molecular Descriptor Space', fontweight='bold', pad=20)
ax6.legend(fontsize=8)
ax6.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('drug_optimization_results.png', dpi=300, bbox_inches='tight')
plt.show()

# Generate Pareto front exploration (Affinity vs Off-target tradeoff)
print("\n" + "="*70)
print("PARETO FRONT ANALYSIS: Efficacy vs Safety Tradeoff")
print("="*70)

n_points = 20
weight_range = np.linspace(0, 1, n_points)
pareto_results = []

for w_affinity in weight_range:
    w_offtarget = 1 - w_affinity
    weights = [w_affinity, w_offtarget, 0.5]  # Fixed solubility weight
    
    result = differential_evolution(
        lambda x: combined_objective(x, weights),
        bounds,
        maxiter=500,
        popsize=10,
        seed=42
    )
    
    x_opt = result.x
    pareto_results.append({
        'affinity': binding_affinity(x_opt),
        'offtarget': off_target_binding(x_opt),
        'weight_affinity': w_affinity
    })

# Pareto front visualization
fig2, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Pareto front
affinities = [r['affinity'] for r in pareto_results]
offtargets = [r['offtarget'] for r in pareto_results]
weights_aff = [r['weight_affinity'] for r in pareto_results]

scatter = ax1.scatter(affinities, offtargets, c=weights_aff, 
                     cmap='RdYlGn_r', s=100, edgecolors='black', linewidth=1.5)
ax1.plot(affinities, offtargets, 'b--', alpha=0.3, linewidth=1)
ax1.set_xlabel('Binding Affinity (more negative = stronger)', fontweight='bold')
ax1.set_ylabel('Off-target Binding (lower = safer)', fontweight='bold')
ax1.set_title('Pareto Front: Efficacy vs Safety Tradeoff', fontweight='bold', fontsize=12)
ax1.grid(True, alpha=0.3)
cbar = plt.colorbar(scatter, ax=ax1)
cbar.set_label('Efficacy Weight', fontweight='bold')

# Highlight optimal scenarios
for scenario_name in ['Efficacy-focused', 'Safety-focused', 'Balanced']:
    ax1.scatter(results[scenario_name]['affinity'], 
               results[scenario_name]['offtarget'],
               s=300, marker='*', edgecolors='red', linewidth=2,
               label=scenario_name, zorder=5)
ax1.legend(fontsize=9)

# Weight sensitivity analysis
ax2.plot(weights_aff, affinities, 'o-', label='Binding Affinity', linewidth=2, markersize=6)
ax2_twin = ax2.twinx()
ax2_twin.plot(weights_aff, offtargets, 's-', color='orange', 
             label='Off-target Binding', linewidth=2, markersize=6)

ax2.set_xlabel('Efficacy Weight in Objective', fontweight='bold')
ax2.set_ylabel('Binding Affinity', fontweight='bold', color='blue')
ax2_twin.set_ylabel('Off-target Binding', fontweight='bold', color='orange')
ax2.set_title('Sensitivity to Objective Weighting', fontweight='bold', fontsize=12)
ax2.grid(True, alpha=0.3)
ax2.tick_params(axis='y', labelcolor='blue')
ax2_twin.tick_params(axis='y', labelcolor='orange')

lines1, labels1 = ax2.get_legend_handles_labels()
lines2, labels2 = ax2_twin.get_legend_handles_labels()
ax2.legend(lines1 + lines2, labels1 + labels2, loc='center right', fontsize=9)

plt.tight_layout()
plt.savefig('pareto_front_analysis.png', dpi=300, bbox_inches='tight')
plt.show()

print("\nOptimization complete! Graphs saved and displayed.")
print("="*70)

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize, fsolve

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

# Physical constants
MU = 1.327e20  # Sun's gravitational parameter (m^3/s^2)
AU = 1.496e11  # Astronomical Unit (m)

# True orbital parameters (these are what we're trying to recover)
TRUE_PARAMS = {
    'a': 2.5 * AU,           # Semi-major axis
    'e': 0.15,               # Eccentricity
    'omega': np.pi/4,        # Argument of periapsis
    'M0': 0.0,               # Mean anomaly at epoch
    't0': 0.0                # Epoch time
}

def solve_kepler(M, e, tol=1e-10):
    """
    Solve Kepler's equation: M = E - e*sin(E)
    Using Newton-Raphson method
    """
    E = M  # Initial guess
    for _ in range(50):
        f = E - e * np.sin(E) - M
        fp = 1 - e * np.cos(E)
        E_new = E - f / fp
        if abs(E_new - E) < tol:
            return E_new
        E = E_new
    return E

def orbital_position(t, a, e, omega, M0, t0):
    """
    Calculate the position (x, y) at time t given orbital parameters
    """
    # Mean motion
    n = np.sqrt(MU / a**3)
    
    # Mean anomaly at time t
    M = M0 + n * (t - t0)
    
    # Solve Kepler's equation for eccentric anomaly
    E = solve_kepler(M, e)
    
    # Position in orbital plane
    x_orb = a * (np.cos(E) - e)
    y_orb = a * np.sqrt(1 - e**2) * np.sin(E)
    
    # Rotate by argument of periapsis
    x = x_orb * np.cos(omega) - y_orb * np.sin(omega)
    y = x_orb * np.sin(omega) + y_orb * np.cos(omega)
    
    return np.array([x, y])

# Generate synthetic observations
N_OBS = 20
time_span = 365.25 * 24 * 3600  # 1 year in seconds
obs_times = np.linspace(0, time_span, N_OBS)

# True positions
true_positions = np.array([orbital_position(t, TRUE_PARAMS['a'], TRUE_PARAMS['e'], 
                                            TRUE_PARAMS['omega'], TRUE_PARAMS['M0'], 
                                            TRUE_PARAMS['t0']) 
                          for t in obs_times])

# Add observation noise
NOISE_LEVEL = 0.05 * AU  # 0.05 AU noise
observed_positions = true_positions + np.random.normal(0, NOISE_LEVEL, true_positions.shape)

def residuals(params, times, observations):
    """
    Calculate residuals between observed and predicted positions
    """
    a, e, omega, M0 = params
    
    # Constraint: eccentricity must be in [0, 1)
    if e < 0 or e >= 1:
        return np.inf
    
    # Constraint: semi-major axis must be positive
    if a <= 0:
        return np.inf
    
    predicted = np.array([orbital_position(t, a, e, omega, M0, TRUE_PARAMS['t0']) 
                         for t in times])
    
    diff = observations - predicted
    return np.sum(diff**2)  # Sum of squared residuals

def negative_log_likelihood(params, times, observations, sigma):
    """
    Negative log-likelihood function (assuming Gaussian noise)
    This is equivalent to weighted sum of squared residuals
    """
    ssr = residuals(params, times, observations)
    n = len(observations)
    return n * np.log(2 * np.pi * sigma**2) + ssr / sigma**2

# Initial guess (deliberately different from true values)
initial_guess = [2.8 * AU, 0.2, np.pi/3, 0.5]

# Perform optimization (Maximum Likelihood Estimation)
print("=" * 60)
print("ORBITAL DETERMINATION VIA MAXIMUM LIKELIHOOD ESTIMATION")
print("=" * 60)
print("\nTrue Orbital Parameters:")
print(f"  Semi-major axis (a): {TRUE_PARAMS['a']/AU:.4f} AU")
print(f"  Eccentricity (e):    {TRUE_PARAMS['e']:.4f}")
print(f"  Arg. of periapsis:   {TRUE_PARAMS['omega']:.4f} rad ({np.degrees(TRUE_PARAMS['omega']):.2f}°)")
print(f"  Mean anomaly (M0):   {TRUE_PARAMS['M0']:.4f} rad")

print("\nInitial Guess:")
print(f"  Semi-major axis (a): {initial_guess[0]/AU:.4f} AU")
print(f"  Eccentricity (e):    {initial_guess[1]:.4f}")
print(f"  Arg. of periapsis:   {initial_guess[2]:.4f} rad ({np.degrees(initial_guess[2]):.2f}°)")
print(f"  Mean anomaly (M0):   {initial_guess[3]:.4f} rad")

print("\nOptimizing...")
result = minimize(lambda p: negative_log_likelihood(p, obs_times, observed_positions, NOISE_LEVEL),
                 initial_guess,
                 method='Nelder-Mead',
                 options={'maxiter': 10000, 'xatol': 1e-8, 'fatol': 1e-8})

estimated_params = result.x

print("\nEstimated Orbital Parameters:")
print(f"  Semi-major axis (a): {estimated_params[0]/AU:.4f} AU")
print(f"  Eccentricity (e):    {estimated_params[1]:.4f}")
print(f"  Arg. of periapsis:   {estimated_params[2]:.4f} rad ({np.degrees(estimated_params[2]):.2f}°)")
print(f"  Mean anomaly (M0):   {estimated_params[3]:.4f} rad")

print("\nEstimation Errors:")
print(f"  Δa: {(estimated_params[0] - TRUE_PARAMS['a'])/AU:.6f} AU ({100*(estimated_params[0] - TRUE_PARAMS['a'])/TRUE_PARAMS['a']:.2f}%)")
print(f"  Δe: {estimated_params[1] - TRUE_PARAMS['e']:.6f} ({100*(estimated_params[1] - TRUE_PARAMS['e'])/TRUE_PARAMS['e']:.2f}%)")
print(f"  Δω: {estimated_params[2] - TRUE_PARAMS['omega']:.6f} rad ({np.degrees(estimated_params[2] - TRUE_PARAMS['omega']):.2f}°)")
print(f"  ΔM0: {estimated_params[3] - TRUE_PARAMS['M0']:.6f} rad")

print(f"\nOptimization converged: {result.success}")
print(f"Number of iterations: {result.nit}")
print(f"Final negative log-likelihood: {result.fun:.2f}")

# Generate predicted orbit with estimated parameters
n_points = 200
orbit_times = np.linspace(0, time_span, n_points)

true_orbit = np.array([orbital_position(t, TRUE_PARAMS['a'], TRUE_PARAMS['e'], 
                                        TRUE_PARAMS['omega'], TRUE_PARAMS['M0'], 
                                        TRUE_PARAMS['t0']) 
                      for t in orbit_times])

estimated_orbit = np.array([orbital_position(t, estimated_params[0], estimated_params[1], 
                                             estimated_params[2], estimated_params[3], 
                                             TRUE_PARAMS['t0']) 
                           for t in orbit_times])

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

# Plot 1: Orbits comparison
ax1 = plt.subplot(2, 2, 1)
ax1.plot(true_orbit[:, 0]/AU, true_orbit[:, 1]/AU, 'b-', linewidth=2, label='True Orbit', alpha=0.7)
ax1.plot(estimated_orbit[:, 0]/AU, estimated_orbit[:, 1]/AU, 'r--', linewidth=2, label='Estimated Orbit', alpha=0.7)
ax1.scatter(observed_positions[:, 0]/AU, observed_positions[:, 1]/AU, c='green', s=100, 
           marker='o', edgecolors='black', linewidths=1.5, label='Observations', zorder=5, alpha=0.8)
ax1.scatter(0, 0, c='orange', s=300, marker='*', edgecolors='black', linewidths=2, label='Sun', zorder=10)
ax1.set_xlabel('X Position (AU)', fontsize=12, fontweight='bold')
ax1.set_ylabel('Y Position (AU)', fontsize=12, fontweight='bold')
ax1.set_title('Orbital Determination: True vs Estimated Orbit', fontsize=14, fontweight='bold')
ax1.legend(loc='upper right', fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.axis('equal')
ax1.set_xlim(-3.5, 3.5)
ax1.set_ylim(-3.5, 3.5)

# Plot 2: Residuals over time
ax2 = plt.subplot(2, 2, 2)
predicted_at_obs = np.array([orbital_position(t, estimated_params[0], estimated_params[1], 
                                               estimated_params[2], estimated_params[3], 
                                               TRUE_PARAMS['t0']) 
                            for t in obs_times])
residual_distances = np.linalg.norm(observed_positions - predicted_at_obs, axis=1) / AU
ax2.plot(obs_times/(24*3600), residual_distances, 'ro-', linewidth=2, markersize=8, alpha=0.7)
ax2.axhline(y=NOISE_LEVEL/AU, color='gray', linestyle='--', linewidth=2, label=f'Noise Level ({NOISE_LEVEL/AU:.3f} AU)')
ax2.set_xlabel('Time (days)', fontsize=12, fontweight='bold')
ax2.set_ylabel('Residual Distance (AU)', fontsize=12, fontweight='bold')
ax2.set_title('Observation Residuals', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# Plot 3: Parameter comparison
ax3 = plt.subplot(2, 2, 3)
params_names = ['a (AU)', 'e', 'ω (rad)', 'M₀ (rad)']
true_vals = [TRUE_PARAMS['a']/AU, TRUE_PARAMS['e'], TRUE_PARAMS['omega'], TRUE_PARAMS['M0']]
est_vals = [estimated_params[0]/AU, estimated_params[1], estimated_params[2], estimated_params[3]]
x_pos = np.arange(len(params_names))
width = 0.35
ax3.bar(x_pos - width/2, true_vals, width, label='True', color='blue', alpha=0.7, edgecolor='black')
ax3.bar(x_pos + width/2, est_vals, width, label='Estimated', color='red', alpha=0.7, edgecolor='black')
ax3.set_ylabel('Parameter Value', fontsize=12, fontweight='bold')
ax3.set_title('Orbital Parameters: True vs Estimated', fontsize=14, fontweight='bold')
ax3.set_xticks(x_pos)
ax3.set_xticklabels(params_names, fontsize=11)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3, axis='y')

# Plot 4: Zoomed view of observation region
ax4 = plt.subplot(2, 2, 4)
# Find the range of observations
x_min, x_max = observed_positions[:, 0].min()/AU, observed_positions[:, 0].max()/AU
y_min, y_max = observed_positions[:, 1].min()/AU, observed_positions[:, 1].max()/AU
margin = 0.3
ax4.plot(true_orbit[:, 0]/AU, true_orbit[:, 1]/AU, 'b-', linewidth=2, label='True Orbit', alpha=0.7)
ax4.plot(estimated_orbit[:, 0]/AU, estimated_orbit[:, 1]/AU, 'r--', linewidth=2, label='Estimated Orbit', alpha=0.7)
ax4.scatter(observed_positions[:, 0]/AU, observed_positions[:, 1]/AU, c='green', s=120, 
           marker='o', edgecolors='black', linewidths=1.5, label='Observations', zorder=5, alpha=0.8)
# Draw error bars
for i, obs in enumerate(observed_positions):
    pred = predicted_at_obs[i]
    ax4.plot([obs[0]/AU, pred[0]/AU], [obs[1]/AU, pred[1]/AU], 'k-', alpha=0.3, linewidth=1)
ax4.set_xlabel('X Position (AU)', fontsize=12, fontweight='bold')
ax4.set_ylabel('Y Position (AU)', fontsize=12, fontweight='bold')
ax4.set_title('Zoomed View: Observations with Fit', fontsize=14, fontweight='bold')
ax4.legend(loc='best', fontsize=10)
ax4.grid(True, alpha=0.3)
ax4.set_xlim(x_min - margin, x_max + margin)
ax4.set_ylim(y_min - margin, y_max + margin)

plt.tight_layout()
plt.savefig('orbital_determination.png', dpi=150, bbox_inches='tight')
plt.show()

print("\n" + "=" * 60)
print("Visualization saved as 'orbital_determination.png'")
print("=" * 60)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
import seaborn as sns
from matplotlib import cm
from itertools import product

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

# Define factor ranges
factors = {
    'Temperature (°C)': [50, 90],
    'Catalyst (%)': [0.5, 2.5],
    'Time (min)': [30, 90]
}

# Create Central Composite Design manually
def create_ccd(n_factors, alpha=1.682, center_points=4):
    """
    Create a Central Composite Design
    
    Parameters:
    - n_factors: number of factors
    - alpha: distance of axial points from center (1.682 for orthogonal with 3 factors)
    - center_points: number of center point replicates
    """
    # Factorial points (2^k corners)
    factorial = np.array(list(product([-1, 1], repeat=n_factors)))
    
    # Axial points (2*k star points)
    axial = np.zeros((2*n_factors, n_factors))
    for i in range(n_factors):
        axial[2*i, i] = -alpha
        axial[2*i+1, i] = alpha
    
    # Center points
    center = np.zeros((center_points, n_factors))
    
    # Combine all points
    design = np.vstack([factorial, axial, center])
    
    return design

n_factors = 3
design_matrix = create_ccd(n_factors, alpha=1.682, center_points=4)

print("=" * 70)
print("CENTRAL COMPOSITE DESIGN (CCD) MATRIX")
print("=" * 70)
print(f"Number of factors: {n_factors}")
print(f"Number of experimental runs: {len(design_matrix)}")
print(f"\nDesign breakdown:")
print(f"  - Factorial points (corners): {2**n_factors}")
print(f"  - Axial points (star): {2*n_factors}")
print(f"  - Center points: 4")
print(f"  - Total: {len(design_matrix)}")
print(f"\nDesign Matrix (coded values -1.682 to +1.682):")
print(design_matrix)
print("=" * 70)

# Convert coded values to actual experimental conditions
def decode_values(coded_matrix, factor_ranges):
    decoded = np.zeros_like(coded_matrix)
    for i, (factor_name, (low, high)) in enumerate(factor_ranges.items()):
        center = (high + low) / 2
        half_range = (high - low) / 2
        decoded[:, i] = center + coded_matrix[:, i] * half_range
    return decoded

actual_conditions = decode_values(design_matrix, factors)

# Create DataFrame for better visualization
df_design = pd.DataFrame(
    actual_conditions,
    columns=list(factors.keys())
)
df_design.index = [f'Run {i+1}' for i in range(len(df_design))]

print("\n" + "=" * 70)
print("EXPERIMENTAL CONDITIONS (Actual Values)")
print("=" * 70)
print(df_design.round(2))
print("=" * 70)

# Simulate experimental results
# True underlying model: Yield = 60 + 8*X1 + 5*X2 + 3*X3 - 2*X1*X2 - 3*X1^2 - 2*X2^2 + noise
def true_yield_function(X):
    """
    Simulated true yield function with quadratic terms and interactions
    X: coded values (-1.682 to +1.682)
    """
    X1, X2, X3 = X[:, 0], X[:, 1], X[:, 2]
    yield_value = (60 + 8*X1 + 5*X2 + 3*X3 
                   - 2*X1*X2  # Interaction term
                   - 3*X1**2  # Quadratic term
                   - 2*X2**2  # Quadratic term
                   + np.random.normal(0, 1.5, len(X)))  # Experimental noise
    return yield_value

# Generate experimental responses
yields = true_yield_function(design_matrix)
df_design['Yield (%)'] = yields.round(2)

print("\n" + "=" * 70)
print("EXPERIMENTAL RESULTS")
print("=" * 70)
print(df_design.round(2))
print("=" * 70)

# Calculate summary statistics
print("\n" + "=" * 70)
print("SUMMARY STATISTICS")
print("=" * 70)
print(f"Mean Yield: {yields.mean():.2f}%")
print(f"Std Dev: {yields.std():.2f}%")
print(f"Min Yield: {yields.min():.2f}%")
print(f"Max Yield: {yields.max():.2f}%")
print("=" * 70)

# Build response surface model
poly = PolynomialFeatures(degree=2, include_bias=True)
X_poly = poly.fit_transform(design_matrix)

model = LinearRegression()
model.fit(X_poly, yields)

# Get feature names
feature_names = poly.get_feature_names_out(['X1', 'X2', 'X3'])

# Display model coefficients
coefficients = pd.DataFrame({
    'Term': feature_names,
    'Coefficient': np.concatenate([[model.intercept_], model.coef_[1:]])
})

print("\n" + "=" * 70)
print("RESPONSE SURFACE MODEL COEFFICIENTS")
print("=" * 70)
print("Model Equation:")
print("Yield = β₀ + β₁X₁ + β₂X₂ + β₃X₃ + β₁₂X₁X₂ + β₁₃X₁X₃ + β₂₃X₂X₃")
print("        + β₁₁X₁² + β₂₂X₂² + β₃₃X₃²")
print("\nWhere:")
print("  X₁ = Temperature (coded)")
print("  X₂ = Catalyst Concentration (coded)")
print("  X₃ = Reaction Time (coded)")
print("\n" + coefficients.to_string(index=False))
print("=" * 70)

# Calculate R-squared
y_pred = model.predict(X_poly)
r_squared = model.score(X_poly, yields)
print(f"\nModel R² Score: {r_squared:.4f}")
print(f"Adjusted R² Score: {1 - (1-r_squared)*(len(yields)-1)/(len(yields)-X_poly.shape[1]):.4f}")
print("=" * 70)

# Find optimal conditions
# Create a fine grid for optimization
grid_size = 20
x1_grid = np.linspace(-1.682, 1.682, grid_size)
x2_grid = np.linspace(-1.682, 1.682, grid_size)
x3_grid = np.linspace(-1.682, 1.682, grid_size)

# Search for maximum yield
max_yield = -np.inf
optimal_coded = None

for x1 in x1_grid:
    for x2 in x2_grid:
        for x3 in x3_grid:
            X_test = np.array([[x1, x2, x3]])
            X_test_poly = poly.transform(X_test)
            predicted_yield = model.predict(X_test_poly)[0]
            if predicted_yield > max_yield:
                max_yield = predicted_yield
                optimal_coded = X_test[0]

optimal_actual = decode_values(optimal_coded.reshape(1, -1), factors)[0]

print("\n" + "=" * 70)
print("OPTIMIZATION RESULTS")
print("=" * 70)
print("Optimal Conditions (Coded Values):")
print(f"  X₁ (Temperature): {optimal_coded[0]:.3f}")
print(f"  X₂ (Catalyst): {optimal_coded[1]:.3f}")
print(f"  X₃ (Time): {optimal_coded[2]:.3f}")
print("\nOptimal Conditions (Actual Values):")
print(f"  Temperature: {optimal_actual[0]:.2f}°C")
print(f"  Catalyst: {optimal_actual[1]:.2f}%")
print(f"  Time: {optimal_actual[2]:.2f} min")
print(f"\nPredicted Maximum Yield: {max_yield:.2f}%")
print("=" * 70)

# Visualization
fig = plt.figure(figsize=(18, 12))

# 1. Main Effects Plot
ax1 = fig.add_subplot(2, 3, 1)
for i, (factor_name, factor_col) in enumerate(zip(factors.keys(), df_design.columns[:3])):
    sorted_idx = df_design[factor_col].argsort()
    ax1.plot(df_design[factor_col].iloc[sorted_idx], 
             df_design['Yield (%)'].iloc[sorted_idx], 
             'o-', label=factor_name, linewidth=2, markersize=6, alpha=0.7)
ax1.set_xlabel('Factor Value', fontsize=11, fontweight='bold')
ax1.set_ylabel('Yield (%)', fontsize=11, fontweight='bold')
ax1.set_title('Main Effects on Yield', fontsize=13, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 2. Residuals Plot
ax2 = fig.add_subplot(2, 3, 2)
residuals = yields - y_pred
ax2.scatter(y_pred, residuals, alpha=0.6, s=80, edgecolors='k', linewidth=0.5)
ax2.axhline(y=0, color='r', linestyle='--', linewidth=2)
ax2.set_xlabel('Predicted Yield (%)', fontsize=11, fontweight='bold')
ax2.set_ylabel('Residuals', fontsize=11, fontweight='bold')
ax2.set_title('Residual Plot (Should be Random)', fontsize=13, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 3. Predicted vs Actual
ax3 = fig.add_subplot(2, 3, 3)
ax3.scatter(yields, y_pred, alpha=0.6, s=80, edgecolors='k', linewidth=0.5)
min_val, max_val = yields.min(), yields.max()
ax3.plot([min_val, max_val], [min_val, max_val], 'r--', linewidth=2, label='Perfect Fit')
ax3.set_xlabel('Actual Yield (%)', fontsize=11, fontweight='bold')
ax3.set_ylabel('Predicted Yield (%)', fontsize=11, fontweight='bold')
ax3.set_title(f'Predicted vs Actual (R²={r_squared:.4f})', fontsize=13, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

# 4. 3D Response Surface (Temperature vs Catalyst)
ax4 = fig.add_subplot(2, 3, 4, projection='3d')
x1_surf = np.linspace(-1.682, 1.682, 30)
x2_surf = np.linspace(-1.682, 1.682, 30)
X1_mesh, X2_mesh = np.meshgrid(x1_surf, x2_surf)
X3_fixed = 0  # Fix time at center point

X_surf = np.column_stack([X1_mesh.ravel(), X2_mesh.ravel(), 
                          np.full(X1_mesh.ravel().shape, X3_fixed)])
X_surf_poly = poly.transform(X_surf)
Z = model.predict(X_surf_poly).reshape(X1_mesh.shape)

surf = ax4.plot_surface(X1_mesh, X2_mesh, Z, cmap='viridis', alpha=0.8, 
                        edgecolor='none', antialiased=True)
ax4.scatter(design_matrix[:, 0], design_matrix[:, 1], yields, 
           c='red', s=50, alpha=0.8, edgecolors='k', linewidth=0.5)
ax4.set_xlabel('Temperature (coded)', fontsize=10, fontweight='bold')
ax4.set_ylabel('Catalyst (coded)', fontsize=10, fontweight='bold')
ax4.set_zlabel('Yield (%)', fontsize=10, fontweight='bold')
ax4.set_title('Response Surface\n(Time at center)', fontsize=12, fontweight='bold')
fig.colorbar(surf, ax=ax4, shrink=0.5)

# 5. Contour Plot (Temperature vs Catalyst)
ax5 = fig.add_subplot(2, 3, 5)
contour = ax5.contourf(X1_mesh, X2_mesh, Z, levels=20, cmap='viridis', alpha=0.8)
ax5.scatter(design_matrix[:, 0], design_matrix[:, 1], 
           c='red', s=60, alpha=0.9, edgecolors='white', linewidth=1.5, 
           label='Exp. Points')
ax5.scatter(optimal_coded[0], optimal_coded[1], 
           c='yellow', s=200, marker='*', edgecolors='black', linewidth=2,
           label='Optimum', zorder=5)
ax5.set_xlabel('Temperature (coded)', fontsize=11, fontweight='bold')
ax5.set_ylabel('Catalyst (coded)', fontsize=11, fontweight='bold')
ax5.set_title('Contour Plot (Time at center)', fontsize=13, fontweight='bold')
ax5.legend()
fig.colorbar(contour, ax=ax5)

# 6. Coefficient Importance
ax6 = fig.add_subplot(2, 3, 6)
coef_importance = coefficients[coefficients['Term'] != '1'].copy()
coef_importance['Abs_Coef'] = np.abs(coef_importance['Coefficient'])
coef_importance = coef_importance.sort_values('Abs_Coef', ascending=True)

colors = ['green' if x > 0 else 'red' for x in coef_importance['Coefficient']]
ax6.barh(coef_importance['Term'], coef_importance['Coefficient'], color=colors, alpha=0.7, edgecolor='black')
ax6.set_xlabel('Coefficient Value', fontsize=11, fontweight='bold')
ax6.set_ylabel('Model Term', fontsize=11, fontweight='bold')
ax6.set_title('Coefficient Magnitudes\n(Green=Positive, Red=Negative)', fontsize=13, fontweight='bold')
ax6.axvline(x=0, color='black', linestyle='-', linewidth=0.8)
ax6.grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.savefig('doe_analysis.png', dpi=300, bbox_inches='tight')
plt.show()

print("\n" + "=" * 70)
print("INTERPRETATION GUIDE")
print("=" * 70)
print("1. Main Effects Plot: Shows how each factor affects yield individually")
print("2. Residual Plot: Random scatter indicates good model fit")
print("3. Predicted vs Actual: Points near diagonal line = accurate predictions")
print("4. 3D Response Surface: Visualizes yield landscape")
print("5. Contour Plot: Top view showing yield levels, star = optimal point")
print("6. Coefficient Bar Chart: Larger bars = stronger effects")
print("=" * 70)
print("\nAll visualizations generated successfully!")
print("=" * 70)