|
| 1 | +import numpy as np |
| 2 | +import matplotlib.pyplot as plt |
| 3 | +import os |
| 4 | +import sys |
| 5 | + |
| 6 | +# Add the project root to sys.path so we can import src |
| 7 | +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))) |
| 8 | + |
| 9 | +from src.portfolio.moead import MOEADOptimizer |
| 10 | + |
| 11 | + |
| 12 | +def proj_simple_normalization(weights): |
| 13 | + """ |
| 14 | + Project weights by clipping negatives to 0 and normalizing the sum. |
| 15 | + """ |
| 16 | + weights = np.maximum(weights, 0) |
| 17 | + s = np.sum(weights) |
| 18 | + if s == 0: |
| 19 | + return np.ones(len(weights)) / len(weights) |
| 20 | + return weights / s |
| 21 | + |
| 22 | + |
| 23 | +def map_to_2d(weights): |
| 24 | + w2 = weights[:, 1] |
| 25 | + w3 = weights[:, 2] |
| 26 | + |
| 27 | + x = w2 + 0.5 * w3 |
| 28 | + y = (np.sqrt(3.0) / 2.0) * w3 |
| 29 | + return x, y |
| 30 | + |
| 31 | + |
| 32 | +def plot_triangle(ax): |
| 33 | + triangle = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3) / 2], [0, 0]]) |
| 34 | + ax.plot(triangle[:, 0], triangle[:, 1], "k-", lw=1.5) |
| 35 | + ax.set_aspect("equal") |
| 36 | + ax.axis("off") |
| 37 | + |
| 38 | + |
| 39 | +def run_simulation(): |
| 40 | + np.random.seed(42) |
| 41 | + N_points = 2000 |
| 42 | + N_iterations = 200 # Increased to allow operators to reach stationary distribution |
| 43 | + |
| 44 | + # 1. Start with a uniform distribution of points on the 3-simplex |
| 45 | + initial_points = np.random.dirichlet(np.ones(3), N_points) |
| 46 | + |
| 47 | + points_euclidean = initial_points.copy() |
| 48 | + points_norm = initial_points.copy() |
| 49 | + |
| 50 | + optimizer = MOEADOptimizer(problem=None) |
| 51 | + |
| 52 | + # 2. Iterate |
| 53 | + for i in range(N_iterations): |
| 54 | + # Create random pairs for crossover |
| 55 | + indices = np.random.permutation(N_points) |
| 56 | + p1_idx = indices[: N_points // 2] |
| 57 | + p2_idx = indices[N_points // 2 :] |
| 58 | + |
| 59 | + # We process Euclidean population |
| 60 | + offspring_euclidean = np.empty((N_points, 3)) |
| 61 | + for j in range(len(p1_idx)): |
| 62 | + p1 = points_euclidean[p1_idx[j]] |
| 63 | + p2 = points_euclidean[p2_idx[j]] |
| 64 | + # Generate two children per pair to keep population size constant |
| 65 | + offspring_euclidean[2 * j] = optimizer._sbx_crossover(p1, p2) |
| 66 | + # Second child: technically MOEAD generates one randomly, we'll force the other to keep pops equal, |
| 67 | + # but simpler to just call it again |
| 68 | + offspring_euclidean[2 * j + 1] = optimizer._sbx_crossover(p1, p2) |
| 69 | + |
| 70 | + for k in range(N_points): |
| 71 | + offspring_euclidean[k] = optimizer._polynomial_mutation( |
| 72 | + offspring_euclidean[k] |
| 73 | + ) |
| 74 | + # Apply Euclidean projection from MOEAD |
| 75 | + points_euclidean[k] = optimizer._repair(offspring_euclidean[k]) |
| 76 | + |
| 77 | + # We process Normalization population |
| 78 | + offspring_norm = np.empty((N_points, 3)) |
| 79 | + for j in range(len(p1_idx)): |
| 80 | + p1 = points_norm[p1_idx[j]] |
| 81 | + p2 = points_norm[p2_idx[j]] |
| 82 | + offspring_norm[2 * j] = optimizer._sbx_crossover(p1, p2) |
| 83 | + offspring_norm[2 * j + 1] = optimizer._sbx_crossover(p1, p2) |
| 84 | + |
| 85 | + for k in range(N_points): |
| 86 | + offspring_norm[k] = optimizer._polynomial_mutation(offspring_norm[k]) |
| 87 | + points_norm[k] = proj_simple_normalization(offspring_norm[k]) |
| 88 | + |
| 89 | + # 3. Visualization |
| 90 | + fig, axes = plt.subplots(1, 3, figsize=(18, 5)) |
| 91 | + |
| 92 | + # Using hexbin to clearly show the density and fix the "no difference" plotting issue |
| 93 | + |
| 94 | + ax = axes[0] |
| 95 | + plot_triangle(ax) |
| 96 | + x, y = map_to_2d(initial_points) |
| 97 | + hb = ax.hexbin(x, y, gridsize=30, cmap="viridis", mincnt=1) |
| 98 | + fig.colorbar(hb, ax=ax, label="Count") |
| 99 | + ax.set_title("Initial (Uniform)") |
| 100 | + |
| 101 | + ax = axes[1] |
| 102 | + plot_triangle(ax) |
| 103 | + x, y = map_to_2d(points_euclidean) |
| 104 | + hb = ax.hexbin(x, y, gridsize=30, cmap="viridis", mincnt=1) |
| 105 | + fig.colorbar(hb, ax=ax, label="Count") |
| 106 | + ax.set_title( |
| 107 | + f"Euclidean Projection\n(After {N_iterations} iters with SBX & PolyMut)" |
| 108 | + ) |
| 109 | + |
| 110 | + ax = axes[2] |
| 111 | + plot_triangle(ax) |
| 112 | + x, y = map_to_2d(points_norm) |
| 113 | + hb = ax.hexbin(x, y, gridsize=30, cmap="viridis", mincnt=1) |
| 114 | + fig.colorbar(hb, ax=ax, label="Count") |
| 115 | + ax.set_title( |
| 116 | + f"Simple Normalization\n(After {N_iterations} iters with SBX & PolyMut)" |
| 117 | + ) |
| 118 | + |
| 119 | + plt.tight_layout() |
| 120 | + output_path = os.path.join("plots", "repair_bias_simulation.png") |
| 121 | + os.makedirs("plots", exist_ok=True) |
| 122 | + plt.savefig(output_path, dpi=150, bbox_inches="tight") |
| 123 | + print(f"Result saved to {output_path}") |
| 124 | + |
| 125 | + |
| 126 | +if __name__ == "__main__": |
| 127 | + run_simulation() |
0 commit comments