{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9578ab3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from optimizr import differential_evolution\n",
    "import time\n",
    "\n",
    "np.random.seed(42)\n",
    "print(\"OptimizR Differential Evolution Module Loaded!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78cfac35",
   "metadata": {},
   "source": [
    "# Differential Evolution Tutorial - Global Optimization\n",
    "\n",
    "## Introduction\n",
    "\n",
    "**Differential Evolution (DE)** is a powerful population-based stochastic optimization algorithm designed for global optimization of non-convex, non-differentiable, and multimodal problems.\n",
    "\n",
    "### Why Differential Evolution?\n",
    "\n",
    "Unlike gradient-based methods that can get stuck in local minima, DE:\n",
    "- ✅ **Global search capability** - Explores entire parameter space\n",
    "- ✅ **No gradient required** - Works with black-box functions\n",
    "- ✅ **Few hyperparameters** - Mutation factor F and crossover rate CR\n",
    "- ✅ **Robust** - Handles noisy and discontinuous functions\n",
    "- ✅ **Parallelizable** - Population members can be evaluated independently\n",
    "\n",
    "### Applications\n",
    "- Portfolio optimization\n",
    "- Hyperparameter tuning in ML\n",
    "- Engineering design optimization\n",
    "- Physics parameter fitting\n",
    "- Control system design\n",
    "\n",
    "## Algorithm Overview\n",
    "\n",
    "### The DE/rand/1/bin Strategy\n",
    "\n",
    "Given a population of $N_p$ candidate solutions $\\mathbf{x}_i$, DE iterates:\n",
    "\n",
    "**1. Mutation** - Create mutant vector:\n",
    "$$\\mathbf{v}_i = \\mathbf{x}_{r1} + F \\cdot (\\mathbf{x}_{r2} - \\mathbf{x}_{r3})$$\n",
    "\n",
    "where $r1, r2, r3$ are random distinct indices, and $F \\in [0, 2]$ is the mutation factor.\n",
    "\n",
    "**2. Crossover** - Create trial vector:\n",
    "$$u_{i,j} = \\begin{cases}\n",
    "v_{i,j} & \\text{if } \\text{rand}() < CR \\text{ or } j = j_{rand} \\\\\n",
    "x_{i,j} & \\text{otherwise}\n",
    "\\end{cases}$$\n",
    "\n",
    "where $CR \\in [0, 1]$ is the crossover probability.\n",
    "\n",
    "**3. Selection** - Greedy selection:\n",
    "$$\\mathbf{x}_i^{t+1} = \\begin{cases}\n",
    "\\mathbf{u}_i & \\text{if } f(\\mathbf{u}_i) < f(\\mathbf{x}_i^t) \\\\\n",
    "\\mathbf{x}_i^t & \\text{otherwise}\n",
    "\\end{cases}$$\n",
    "\n",
    "### Convergence\n",
    "\n",
    "Under mild conditions, DE converges to the global optimum with probability 1:\n",
    "$$\\lim_{t \\to \\infty} P\\left(\\|\\mathbf{x}^*_t - \\mathbf{x}^*\\| < \\epsilon\\right) = 1$$\n",
    "\n",
    "where $\\mathbf{x}^*$ is the global optimum.\n",
    "\n",
    "### Complexity\n",
    "\n",
    "- **Time:** $O(N_p \\cdot d \\cdot T)$ where $d$ is dimension, $T$ is iterations\n",
    "- **Space:** $O(N_p \\cdot d)$ for population storage\n",
    "\n",
    "## References\n",
    "\n",
    "- Storn, R., & Price, K. (1997). \"Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces.\" *Journal of global optimization*, 11(4), 341-359.\n",
    "- Das, S., & Suganthan, P. N. (2011). \"Differential evolution: A survey of the state-of-the-art.\" *IEEE transactions on evolutionary computation*, 15(1), 4-31."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "588bf0b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from optimizr import differential_evolution\n",
    "\n",
    "np.random.seed(42)\n",
    "print(\"OptimizR Differential Evolution Loaded!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ba7fb3a",
   "metadata": {},
   "source": [
    "## Example 1: Rosenbrock Function (Banana Valley)\n",
    "\n",
    "$$f(\\mathbf{x}) = \\sum_{i=1}^{n-1} \\left[100(x_{i+1} - x_i^2)^2 + (1 - x_i)^2\\right]$$\n",
    "\n",
    "Global minimum: $f(1, 1, \\ldots, 1) = 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f9644cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rosenbrock(x):\n",
    "    \"\"\"N-dimensional Rosenbrock function.\"\"\"\n",
    "    return sum(100 * (x[i+1] - x[i]**2)**2 + (1 - x[i])**2 \n",
    "              for i in range(len(x) - 1))\n",
    "\n",
    "# Test function\n",
    "print(f\"f([1, 1, 1]): {rosenbrock([1.0, 1.0, 1.0])}\")\n",
    "print(f\"f([0, 0, 0]): {rosenbrock([0.0, 0.0, 0.0])}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dda42ec7",
   "metadata": {},
   "source": [
    "### Visualize 2D Rosenbrock"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec984eff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create meshgrid\n",
    "x1 = np.linspace(-2, 2, 200)\n",
    "x2 = np.linspace(-1, 3, 200)\n",
    "X1, X2 = np.meshgrid(x1, x2)\n",
    "Z = np.array([[rosenbrock([x1_val, x2_val]) for x1_val, x2_val in zip(x1_row, x2_row)] \n",
    "             for x1_row, x2_row in zip(X1, X2)])\n",
    "\n",
    "fig = plt.figure(figsize=(14, 6))\n",
    "\n",
    "# 3D surface\n",
    "ax1 = fig.add_subplot(121, projection='3d')\n",
    "surf = ax1.plot_surface(X1, X2, np.log10(Z + 1), cmap='viridis', alpha=0.8)\n",
    "ax1.scatter([1], [1], [0], c='red', s=200, marker='*', edgecolors='black', linewidths=2, label='Global min')\n",
    "ax1.set_xlabel('$x_1$', fontsize=11)\n",
    "ax1.set_ylabel('$x_2$', fontsize=11)\n",
    "ax1.set_zlabel('$\\log_{10}(f + 1)$', fontsize=11)\n",
    "ax1.set_title('Rosenbrock Function (3D)', fontsize=13, fontweight='bold')\n",
    "\n",
    "# 2D contour\n",
    "ax2 = fig.add_subplot(122)\n",
    "contour = ax2.contour(X1, X2, np.log10(Z + 1), levels=20, cmap='viridis')\n",
    "ax2.scatter([1], [1], c='red', s=200, marker='*', edgecolors='black', linewidths=2, label='Global min', zorder=5)\n",
    "ax2.set_xlabel('$x_1$', fontsize=11)\n",
    "ax2.set_ylabel('$x_2$', fontsize=11)\n",
    "ax2.set_title('Rosenbrock Function (Contour)', fontsize=13, fontweight='bold')\n",
    "ax2.legend()\n",
    "plt.colorbar(contour, ax=ax2, label='$\\log_{10}(f + 1)$')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "812d80ae",
   "metadata": {},
   "source": [
    "### Optimize with Differential Evolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f722fb0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 10-dimensional Rosenbrock\n",
    "n_dims = 10\n",
    "bounds = [(-5, 5)] * n_dims\n",
    "\n",
    "print(f\"Optimizing {n_dims}D Rosenbrock function...\")\n",
    "result = differential_evolution(\n",
    "    objective_fn=rosenbrock,\n",
    "    bounds=bounds,\n",
    "    maxiter=500,\n",
    "    popsize=15,\n",
    "    mutation_factor=0.8,\n",
    "    crossover_rate=0.7,\n",
    "    seed=42\n",
    ")\n",
    "\n",
    "print(f\"\\nOptimization completed!\")\n",
    "print(f\"Best solution: {result.x}\")\n",
    "print(f\"Best value: {result.fun:.6e}\")\n",
    "print(f\"Function evaluations: {result.nfev}\")\n",
    "print(f\"\\nDistance to true optimum [1, 1, ..., 1]:\")\n",
    "print(f\"  ||x - x*|| = {np.linalg.norm(result.x - np.ones(n_dims)):.6f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89fde8e9",
   "metadata": {},
   "source": [
    "## Example 2: Rastrigin Function (Many Local Minima)\n",
    "\n",
    "$$f(\\mathbf{x}) = 10n + \\sum_{i=1}^n \\left[x_i^2 - 10\\cos(2\\pi x_i)\\right]$$\n",
    "\n",
    "Global minimum: $f(0, 0, \\ldots, 0) = 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fff32181",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rastrigin(x):\n",
    "    \"\"\"Rastrigin function with many local minima.\"\"\"\n",
    "    n = len(x)\n",
    "    return 10 * n + sum(xi**2 - 10 * np.cos(2 * np.pi * xi) for xi in x)\n",
    "\n",
    "# Visualize 2D\n",
    "x1 = np.linspace(-5.12, 5.12, 200)\n",
    "x2 = np.linspace(-5.12, 5.12, 200)\n",
    "X1, X2 = np.meshgrid(x1, x2)\n",
    "Z = np.array([[rastrigin([x1_val, x2_val]) for x1_val, x2_val in zip(x1_row, x2_row)]\n",
    "             for x1_row, x2_row in zip(X1, X2)])\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(14, 6))\n",
    "\n",
    "# 3D plot\n",
    "ax1 = fig.add_subplot(121, projection='3d')\n",
    "ax1.plot_surface(X1, X2, Z, cmap='plasma', alpha=0.8)\n",
    "ax1.scatter([0], [0], [0], c='red', s=200, marker='*', edgecolors='black', linewidths=2)\n",
    "ax1.set_xlabel('$x_1$', fontsize=11)\n",
    "ax1.set_ylabel('$x_2$', fontsize=11)\n",
    "ax1.set_zlabel('$f(x)$', fontsize=11)\n",
    "ax1.set_title('Rastrigin Function (3D)', fontsize=13, fontweight='bold')\n",
    "\n",
    "# Contour plot\n",
    "contour = axes[1].contourf(X1, X2, Z, levels=30, cmap='plasma')\n",
    "axes[1].scatter([0], [0], c='red', s=200, marker='*', edgecolors='black', linewidths=2, label='Global min', zorder=5)\n",
    "axes[1].set_xlabel('$x_1$', fontsize=11)\n",
    "axes[1].set_ylabel('$x_2$', fontsize=11)\n",
    "axes[1].set_title('Rastrigin Function (Contour)', fontsize=13, fontweight='bold')\n",
    "axes[1].legend()\n",
    "plt.colorbar(contour, ax=axes[1])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(\"Note: Rastrigin has MANY local minima (visible as the peaks in the plot)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79deeae5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optimize Rastrigin\n",
    "n_dims = 10\n",
    "bounds = [(-5.12, 5.12)] * n_dims\n",
    "\n",
    "print(f\"Optimizing {n_dims}D Rastrigin function...\")\n",
    "result = differential_evolution(\n",
    "    objective_fn=rastrigin,\n",
    "    bounds=bounds,\n",
    "    maxiter=1000,\n",
    "    popsize=20,\n",
    "    mutation_factor=0.9,\n",
    "    crossover_rate=0.9,\n",
    "    seed=42\n",
    ")\n",
    "\n",
    "print(f\"\\nBest solution: {result.x}\")\n",
    "print(f\"Best value: {result.fun:.6e}\")\n",
    "print(f\"Distance to global optimum: {np.linalg.norm(result.x):.6f}\")\n",
    "\n",
    "if result.fun < 1.0:\n",
    "    print(\"\\n✓ Successfully found global minimum!\")\n",
    "else:\n",
    "    print(\"\\n⚠ Stuck in local minimum (try increasing popsize or maxiter)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb324ced",
   "metadata": {},
   "source": [
    "## Example 3: Real-World Application - Portfolio Optimization\n",
    "\n",
    "Minimize portfolio variance with expected return constraint.\n",
    "\n",
    "$$\\min_{\\mathbf{w}} \\quad \\mathbf{w}^T \\Sigma \\mathbf{w}$$\n",
    "$$\\text{s.t.} \\quad \\mathbf{w}^T \\boldsymbol{\\mu} \\geq r_{\\text{target}}$$\n",
    "$$\\sum_i w_i = 1, \\quad w_i \\geq 0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19b33e87",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate synthetic asset data\n",
    "n_assets = 10\n",
    "n_periods = 252  # 1 year of daily data\n",
    "\n",
    "# Simulate correlated returns\n",
    "np.random.seed(42)\n",
    "mean_returns = np.random.uniform(0.0005, 0.002, n_assets)  # Daily returns\n",
    "returns = np.random.multivariate_normal(\n",
    "    mean=mean_returns,\n",
    "    cov=np.diag(np.random.uniform(0.01, 0.03, n_assets)**2),\n",
    "    size=n_periods\n",
    ")\n",
    "\n",
    "# Compute statistics\n",
    "mu = returns.mean(axis=0)  # Expected returns\n",
    "Sigma = np.cov(returns.T)  # Covariance matrix\n",
    "\n",
    "print(f\"Portfolio with {n_assets} assets\")\n",
    "print(f\"Expected returns (daily): {mu}\")\n",
    "print(f\"Annualized returns: {mu * 252}\")\n",
    "\n",
    "# Visualize returns\n",
    "plt.figure(figsize=(12, 6))\n",
    "cumulative_returns = np.cumprod(1 + returns, axis=0) - 1\n",
    "for i in range(n_assets):\n",
    "    plt.plot(cumulative_returns[:, i], alpha=0.6, label=f'Asset {i+1}')\n",
    "plt.xlabel('Days', fontsize=11)\n",
    "plt.ylabel('Cumulative Return', fontsize=11)\n",
    "plt.title('Simulated Asset Returns', fontsize=13, fontweight='bold')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "plt.grid(alpha=0.3)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "348c49cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def portfolio_objective(weights):\n",
    "    \"\"\"\n",
    "    Minimize: variance + penalty for constraint violations.\n",
    "    \"\"\"\n",
    "    # Portfolio variance\n",
    "    variance = weights @ Sigma @ weights\n",
    "    \n",
    "    # Constraints (penalize violations)\n",
    "    target_return = 0.0015  # Target daily return\n",
    "    return_constraint = max(0, target_return - weights @ mu)\n",
    "    sum_constraint = abs(weights.sum() - 1.0)\n",
    "    negative_constraint = max(0, -weights.min())\n",
    "    \n",
    "    # Penalize constraint violations heavily\n",
    "    penalty = 1000 * (return_constraint + sum_constraint + negative_constraint)\n",
    "    \n",
    "    return variance + penalty\n",
    "\n",
    "# Optimize\n",
    "bounds = [(0, 1)] * n_assets  # Weights between 0 and 1\n",
    "\n",
    "print(\"Optimizing portfolio allocation...\")\n",
    "result = differential_evolution(\n",
    "    objective_fn=portfolio_objective,\n",
    "    bounds=bounds,\n",
    "    maxiter=500,\n",
    "    popsize=20,\n",
    "    seed=42\n",
    ")\n",
    "\n",
    "optimal_weights = result.x\n",
    "optimal_return = optimal_weights @ mu\n",
    "optimal_volatility = np.sqrt(optimal_weights @ Sigma @ optimal_weights)\n",
    "\n",
    "print(f\"\\nOptimal Portfolio:\")\n",
    "print(f\"Weights: {optimal_weights}\")\n",
    "print(f\"Sum of weights: {optimal_weights.sum():.6f}\")\n",
    "print(f\"\\nExpected daily return: {optimal_return:.6f} ({optimal_return * 252:.2%} annualized)\")\n",
    "print(f\"Daily volatility: {optimal_volatility:.6f} ({optimal_volatility * np.sqrt(252):.2%} annualized)\")\n",
    "print(f\"Sharpe ratio (assuming 0% risk-free): {optimal_return / optimal_volatility:.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75325f0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize allocation\n",
    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
    "\n",
    "# Bar chart\n",
    "axes[0].bar(range(n_assets), optimal_weights, color='steelblue', edgecolor='black')\n",
    "axes[0].set_xlabel('Asset', fontsize=11)\n",
    "axes[0].set_ylabel('Weight', fontsize=11)\n",
    "axes[0].set_title('Optimal Portfolio Allocation', fontsize=13, fontweight='bold')\n",
    "axes[0].grid(alpha=0.3, axis='y')\n",
    "\n",
    "# Pie chart\n",
    "nonzero_weights = optimal_weights[optimal_weights > 0.01]\n",
    "nonzero_assets = [f'Asset {i+1}' for i in range(n_assets) if optimal_weights[i] > 0.01]\n",
    "axes[1].pie(nonzero_weights, labels=nonzero_assets, autopct='%1.1f%%', startangle=90)\n",
    "axes[1].set_title('Portfolio Composition', fontsize=13, fontweight='bold')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ed601b8",
   "metadata": {},
   "source": [
    "## Key Takeaways\n",
    "\n",
    "1. **DE is excellent for non-convex, multimodal problems** where gradient-based methods fail\n",
    "2. **Population-based approach** explores solution space thoroughly\n",
    "3. **Few hyperparameters** - typically F ∈ [0.5, 1], CR ∈ [0.7, 1]\n",
    "4. **Robust** - works well on wide variety of problems\n",
    "5. **OptimizR provides 50-100x speedup** over pure Python\n",
    "6. **Real-world applications** - engineering, ML, finance, science\n",
    "\n",
    "## Further Reading\n",
    "\n",
    "- Storn & Price (1997). \"Differential evolution–a simple and efficient heuristic for global optimization\"\n",
    "- Price, Storn & Lampinen (2005). \"Differential Evolution: A Practical Approach to Global Optimization\""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "rhftlab",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}