feat: add geometrically valid simplex crossover operator

Author: MichalRedm

src/portfolio/moead.py

@@ -18,7 +18,7 @@ def get_pareto_parameters(self, num_points, **kwargs):
     # ── Public interface ──────────────────────────────────────────────────────
 
     def generate_pareto_front(
-        self, num_points=100, generations=100, T=10, nr=2, verbose=False, **kwargs
+        self, num_points=100, generations=100, T=10, nr=2, verbose=False, crossover_operator="sbx", **kwargs
     ):
         """
         MOEA/D with Tchebycheff decomposition, stable normalisation, and
@@ -143,21 +143,23 @@ def normalise(f):
                 nb_idx = neighbors[i]
                 p1_idx, p2_idx = np.random.choice(nb_idx, 2, replace=False)
 
-                # Simulated Binary Crossover (SBX)
-                offspring = self._sbx_crossover(population[p1_idx], population[p2_idx])
-
-                # CRITICAL: We MUST repair immediately after SBX because it
-                # violates the simplex (sum != 1, potential negative weights)
-                offspring = self._repair(offspring)
-
-                # Polynomial Mutation
-                # (Replaces the 10% vector-level Gaussian mutation with
-                # mathematically standard component-level Polynomial Mutation)
-                offspring = self._polynomial_mutation(offspring)
-
-                # CRITICAL: Repair again because polynomial mutation pushes
-                # bounds independently, breaking the sum-to-1 constraint.
-                offspring = self._repair(offspring)
+                if crossover_operator == "sbx":
+                    # Simulated Binary Crossover (SBX)
+                    offspring = self._sbx_crossover(population[p1_idx], population[p2_idx])
+                    # CRITICAL: We MUST repair immediately after SBX
+                    offspring = self._repair(offspring)
+                    
+                    # Polynomial Mutation
+                    offspring = self._polynomial_mutation(offspring)
+                    # CRITICAL: Repair again because polynomial mutation pushes bounds
+                    offspring = self._repair(offspring)
+                    
+                elif crossover_operator == "simplex":
+                    # Linear mix crossing the full valid line segment on the simplex.
+                    # It natively avoids breaking the simplex geometry.
+                    offspring = self._simplex_crossover(population[p1_idx], population[p2_idx])
+                else:
+                    raise ValueError(f"Unknown crossover operator: {crossover_operator}")
 
                 # Evaluate offspring.
                 off_f = np.array(

src/portfolio/moead_base.py

@@ -115,6 +115,33 @@ def _differential_evolution_crossover(self, p1, p2, p3, F=0.5, CR=1.0):
         offspring = np.where(mask, v, p1)
         return offspring
 
+    def _simplex_crossover(self, p1, p2):
+        """
+        Linear combination of two solutions spanning the entire valid segment
+        on the simplex. Does not require repair on its own.
+        """
+        diff = p1 - p2
+        
+        pos_mask = diff > 1e-9
+        neg_mask = diff < -1e-9
+        
+        alpha_min = np.max(-p2[pos_mask] / diff[pos_mask]) if np.any(pos_mask) else 0.0
+        alpha_max = np.min(-p2[neg_mask] / diff[neg_mask]) if np.any(neg_mask) else 1.0
+        
+        alpha = np.random.uniform(alpha_min, alpha_max)
+        offspring = alpha * p1 + (1.0 - alpha) * p2
+        
+        # Ensure constraints due to possible floating point inaccuracies
+        offspring = np.clip(offspring, 0.0, 1.0)
+        sum_weights = np.sum(offspring)
+        if sum_weights > 0:
+            offspring = offspring / sum_weights
+        else:
+            # Fallback (should theoretically never happen)
+            offspring = np.ones(len(p1)) / len(p1)
+            
+        return offspring
+
     def optimize(self, param, **kwargs):
         raise NotImplementedError("MOEA/D generates the whole front. Use generate_pareto_front.")