"""Feasibility-First Search Framework (P3-M2). Pure Python implementation (no numpy/scipy/sklearn dependency). Implements: - LHS (Latin Hypercube Sampling) for initial space coverage - Active learning batch selection (distance + uncertainty heuristic) - Constrained feasibility prioritization - Local trust region refinement - Search state management with checkpoint support Per third-party review: default mode is feasibility-first (constrained Bayesian / active learning), not fixed full-factorial DoE. """ import math import random import hashlib import json from dataclasses import dataclass, field from typing import Dict, List, Optional, Tuple, Any, Callable from datetime import datetime from .l0_prescreening import L0PreScreeningEngine, FeasibilityReport @dataclass class ParameterRange: """Definition of a scan parameter range.""" name: str min_value: float max_value: float step: Optional[float] = None # if set, values are quantized unit: str = "" description: str = "" def normalize(self, value: float) -> float: """Normalize value to [0, 1] range.""" if self.max_value == self.min_value: return 0.5 return (value - self.min_value) / (self.max_value - self.min_value) def denormalize(self, norm: float) -> float: """Denormalize from [0, 1] to actual range.""" value = self.min_value + norm * (self.max_value - self.min_value) if self.step is not None and self.step > 0: value = round(value / self.step) * self.step value = max(self.min_value, min(self.max_value, value)) return value def random_value(self, rng: random.Random) -> float: """Generate a random value within range.""" return self.denormalize(rng.random()) @dataclass class SearchPoint: """A single point in the search space.""" id: int params: Dict[str, float] status: str = "pending" # pending / running / ok / failed / infeasible metrics: Dict[str, float] = field(default_factory=dict) feasibility_report: Optional[Dict[str, Any]] = None surrogate_prediction: Optional[float] = None surrogate_uncertainty: Optional[float] = None batch_id: int = 0 created_at: str = "" def param_hash(self) -> str: """Generate a hash of parameter values for caching/dedup.""" sorted_params = sorted(self.params.items()) param_str = json.dumps(sorted_params, sort_keys=True) return hashlib.md5(param_str.encode()).hexdigest()[:12] @dataclass class SearchState: """Complete state of an adaptive search run.""" run_id: str parameters: List[ParameterRange] points: List[SearchPoint] = field(default_factory=list) current_batch: int = 0 total_budget: int = 80 used_budget: int = 0 batch_size: int = 4 search_method: str = "active_learning" trust_region_active: bool = False trust_region_center: Optional[Dict[str, float]] = None trust_region_radius: float = 0.3 # fraction of normalized range best_feasible_point: Optional[SearchPoint] = None best_objective_value: float = float('inf') objective_metric: str = "tavg_nm" # metric to optimize objective_direction: str = "maximize" # maximize / minimize convergence_status: str = "searching" # searching / converged / stalled / budget_exhausted history: List[Dict[str, Any]] = field(default_factory=list) created_at: str = "" updated_at: str = "" def get_param_names(self) -> List[str]: return [p.name for p in self.parameters] def get_param_by_name(self, name: str) -> Optional[ParameterRange]: for p in self.parameters: if p.name == name: return p return None def get_pending_points(self) -> List[SearchPoint]: return [p for p in self.points if p.status == "pending"] def get_completed_points(self) -> List[SearchPoint]: return [p for p in self.points if p.status in ("ok", "failed", "infeasible")] def get_feasible_points(self) -> List[SearchPoint]: return [p for p in self.points if p.status == "ok" and p.feasibility_report and p.feasibility_report.get("feasible", False)] def remaining_budget(self) -> int: return self.total_budget - self.used_budget class FeasibilityFirstSearch: """Feasibility-first adaptive search engine. Implements the recommended default path from third-party review: L0 pre-screening + initial samples + active learning batch selection + local trust region refinement. """ def __init__( self, parameters: List[ParameterRange], l0_engine: Optional[L0PreScreeningEngine] = None, total_budget: int = 80, batch_size: int = 4, initial_samples: int = 16, objective_metric: str = "tavg_nm", objective_direction: str = "maximize", seed: int = 42, ): self.rng = random.Random(seed) self.parameters = parameters self.l0_engine = l0_engine or L0PreScreeningEngine() self.total_budget = total_budget self.batch_size = batch_size self.initial_samples = initial_samples self.objective_metric = objective_metric self.objective_direction = objective_direction self.state = SearchState( run_id=f"search_{datetime.now().strftime('%Y%m%d_%H%M%S')}", parameters=parameters, total_budget=total_budget, batch_size=batch_size, objective_metric=objective_metric, objective_direction=objective_direction, created_at=datetime.now().isoformat(), ) def latin_hypercube_sample(self, n_samples: int) -> List[Dict[str, float]]: """Generate Latin Hypercube Samples in normalized space. Pure Python implementation (no numpy dependency). """ n_params = len(self.parameters) # Create n_samples intervals per dimension intervals = [] for _ in range(n_params): # Shuffle interval assignments assignments = list(range(n_samples)) self.rng.shuffle(assignments) intervals.append(assignments) samples = [] for i in range(n_samples): sample = {} for j, param in enumerate(self.parameters): # Random point within assigned interval interval_idx = intervals[j][i] norm = (interval_idx + self.rng.random()) / n_samples sample[param.name] = param.denormalize(norm) samples.append(sample) return samples def _min_distance_to_existing(self, params: Dict[str, float], existing: List[Dict[str, float]]) -> float: """Calculate minimum normalized Euclidean distance to existing points.""" if not existing: return float('inf') min_dist = float('inf') for ex in existing: dist_sq = 0.0 for param in self.parameters: v1 = param.normalize(params.get(param.name, param.min_value)) v2 = param.normalize(ex.get(param.name, param.min_value)) dist_sq += (v1 - v2) ** 2 dist = math.sqrt(dist_sq) min_dist = min(min_dist, dist) return min_dist def _is_duplicate(self, params: Dict[str, float], existing_points: List[SearchPoint]) -> bool: """Check if parameter set is a duplicate of an existing point.""" param_hash = hashlib.md5( json.dumps(sorted(params.items()), sort_keys=True).encode() ).hexdigest()[:12] return any(p.param_hash() == param_hash for p in existing_points) def generate_initial_batch(self) -> List[SearchPoint]: """Generate initial LHS samples, filtered by L0 feasibility. Prioritizes feasible points but includes some infeasible for boundary learning (active learning needs both classes). """ # Generate more samples than needed, then select diverse subset n_candidates = max(self.initial_samples * 3, 20) candidates = self.latin_hypercube_sample(n_candidates) # L0 pre-screen all candidates scored = [] for params in candidates: report = self.l0_engine.evaluate(params) scored.append((params, report)) # Separate feasible and infeasible feasible = [(p, r) for p, r in scored if r.feasible] infeasible = [(p, r) for p, r in scored if not r.feasible] # Select diverse subset: mostly feasible, some infeasible for boundary n_feasible = min(len(feasible), int(self.initial_samples * 0.75)) n_infeasible = min(len(infeasible), self.initial_samples - n_feasible) selected = [] selected_params = [] # Greedy max-distance selection for feasible points feasible.sort(key=lambda x: x[1].pass_rate, reverse=True) for params, report in feasible: if len(selected) >= n_feasible: break # Check duplicate against already selected params is_dup = any( hashlib.md5(json.dumps(sorted(p.items()), sort_keys=True).encode()).hexdigest()[:12] == hashlib.md5(json.dumps(sorted(params.items()), sort_keys=True).encode()).hexdigest()[:12] for p, _ in selected ) if not is_dup: dist = self._min_distance_to_existing(params, selected_params) if dist > 0.1 or len(selected) < 3: selected.append((params, report)) selected_params.append(params) # Add infeasible points for boundary learning for params, report in infeasible: if len(selected) >= self.initial_samples: break is_dup = any( hashlib.md5(json.dumps(sorted(p.items()), sort_keys=True).encode()).hexdigest()[:12] == hashlib.md5(json.dumps(sorted(params.items()), sort_keys=True).encode()).hexdigest()[:12] for p, _ in selected ) if not is_dup: selected.append((params, report)) selected_params.append(params) # Create SearchPoint objects points = [] for i, (params, report) in enumerate(selected): point = SearchPoint( id=len(self.state.points) + i, params=params, status="infeasible" if not report.feasible else "pending", feasibility_report=report.to_dict(), batch_id=0, created_at=datetime.now().isoformat(), ) if not report.feasible: # Mark as completed (L0 rejected) without using solver budget point.status = "infeasible" points.append(point) self.state.points.extend(points) self.state.current_batch = 1 self.state.used_budget += sum(1 for p in points if p.status == "pending") self.state.updated_at = datetime.now().isoformat() return [p for p in points if p.status == "pending"] def select_next_batch(self) -> List[SearchPoint]: """Select next batch of points using active learning. Strategy: 1. If feasible points found, activate trust region around best 2. Select points balancing: - Exploration: high distance from existing points - Exploitation: near best feasible point (trust region) - Uncertainty: points near feasibility boundary 3. L0 pre-screen all candidates, reject obviously infeasible """ if self.state.remaining_budget() <= 0: self.state.convergence_status = "budget_exhausted" return [] feasible_points = self.state.get_feasible_points() completed_params = [p.params for p in self.state.get_completed_points()] # Activate trust region if we have enough feasible points if len(feasible_points) >= 3 and not self.state.trust_region_active: # C4 fix: select anchor by optimization direction if self.objective_direction == "minimize": best = min(feasible_points, key=lambda p: p.metrics.get(self.objective_metric, float("inf"))) else: best = max(feasible_points, key=lambda p: p.metrics.get(self.objective_metric, 0)) self.state.trust_region_active = True self.state.trust_region_center = best.params self.state.best_feasible_point = best self.state.history.append({ "event": "trust_region_activated", "center": best.params, "batch": self.state.current_batch, }) # Generate candidates n_candidates = 50 candidates = [] if self.state.trust_region_active and self.state.trust_region_center: # 70% from trust region, 30% from global exploration n_trust = int(n_candidates * 0.7) n_global = n_candidates - n_trust # Trust region samples (Gaussian-like around center) for _ in range(n_trust): params = {} for param in self.parameters: center_norm = param.normalize(self.state.trust_region_center[param.name]) # Sample with decreasing radius as search progresses radius = self.state.trust_region_radius * max(0.3, 1.0 - self.state.used_budget / self.total_budget) sample_norm = center_norm + self.rng.gauss(0, radius * 0.3) sample_norm = max(0.0, min(1.0, sample_norm)) params[param.name] = param.denormalize(sample_norm) candidates.append(params) # Global exploration samples candidates.extend(self.latin_hypercube_sample(n_global)) else: # No trust region yet: full LHS exploration candidates = self.latin_hypercube_sample(n_candidates) # Score and select candidates scored_candidates = [] for params in candidates: if self._is_duplicate(params, self.state.points): continue # L0 pre-screen report = self.l0_engine.evaluate(params) if not report.feasible: continue # Skip obviously infeasible # Distance score (exploration) dist = self._min_distance_to_existing(params, completed_params) dist_score = min(dist / math.sqrt(len(self.parameters)), 1.0) # Trust region proximity score (exploitation) trust_score = 0.0 if self.state.trust_region_center: dist_to_center_sq = 0.0 for param in self.parameters: v1 = param.normalize(params[param.name]) v2 = param.normalize(self.state.trust_region_center[param.name]) dist_to_center_sq += (v1 - v2) ** 2 dist_to_center = math.sqrt(dist_to_center_sq) trust_score = max(0, 1.0 - dist_to_center / self.state.trust_region_radius) # Combined score (balance exploration and exploitation) if self.state.trust_region_active: score = 0.4 * dist_score + 0.6 * trust_score else: score = dist_score scored_candidates.append((params, report, score)) # Sort by score and select top batch_size scored_candidates.sort(key=lambda x: x[2], reverse=True) batch_size = min(self.batch_size, self.state.remaining_budget(), len(scored_candidates)) selected = [] for i in range(batch_size): params, report, score = scored_candidates[i] point = SearchPoint( id=len(self.state.points) + i, params=params, status="pending", feasibility_report=report.to_dict(), batch_id=self.state.current_batch, created_at=datetime.now().isoformat(), ) selected.append(point) self.state.points.extend(selected) self.state.used_budget += len(selected) self.state.current_batch += 1 self.state.updated_at = datetime.now().isoformat() return selected def report_result(self, point_id: int, metrics: Dict[str, float], status: str = "ok") -> None: """Report simulation result for a point. Args: point_id: ID of the search point metrics: Dictionary of metric values status: ok / failed """ for point in self.state.points: if point.id == point_id: point.metrics = metrics point.status = status # Update best feasible point if status == "ok" and self.objective_metric in metrics: value = metrics[self.objective_metric] is_better = ( (self.objective_direction == "maximize" and value > self.state.best_objective_value) or (self.objective_direction == "minimize" and value < self.state.best_objective_value) ) if is_better or self.state.best_feasible_point is None: self.state.best_objective_value = value self.state.best_feasible_point = point # C5 fix: migrate trust region center to new best point if self.state.trust_region_active: self.state.trust_region_center = point.params # Check convergence self._check_convergence() break self.state.updated_at = datetime.now().isoformat() def _check_convergence(self) -> None: """Check if search has converged or stalled.""" feasible = self.state.get_feasible_points() if len(feasible) < 5: return # Check if objective has improved in last N points recent = feasible[-10:] if len(feasible) >= 10 else feasible if len(recent) >= 5: values = [p.metrics.get(self.objective_metric, 0) for p in recent] if self.objective_direction == "maximize": improvement = max(values) - max(values[:-3]) if len(values) > 3 else 0 else: improvement = min(values[:-3]) - min(values) if len(values) > 3 else 0 if improvement < 0.001 and self.state.trust_region_active: # Shrink trust region self.state.trust_region_radius *= 0.7 if self.state.trust_region_radius < 0.05: self.state.convergence_status = "converged" self.state.history.append({ "event": "converged", "reason": "trust_region_shrunk_below_threshold", "batch": self.state.current_batch, }) def get_state_summary(self) -> Dict[str, Any]: """Get a summary of current search state.""" return { "run_id": self.state.run_id, "search_method": self.state.search_method, "convergence_status": self.state.convergence_status, "total_budget": self.state.total_budget, "used_budget": self.state.used_budget, "remaining_budget": self.state.remaining_budget(), "current_batch": self.state.current_batch, "total_points": len(self.state.points), "pending_points": len(self.state.get_pending_points()), "completed_points": len(self.state.get_completed_points()), "feasible_points": len(self.state.get_feasible_points()), "trust_region_active": self.state.trust_region_active, "trust_region_radius": round(self.state.trust_region_radius, 4), "best_objective_value": self.state.best_objective_value if self.state.best_objective_value != float('inf') else None, "best_point_params": self.state.best_feasible_point.params if self.state.best_feasible_point else None, } def export_state(self) -> Dict[str, Any]: """Export full search state for checkpointing.""" return { "state": { "run_id": self.state.run_id, "current_batch": self.state.current_batch, "total_budget": self.state.total_budget, "used_budget": self.state.used_budget, "convergence_status": self.state.convergence_status, "trust_region_active": self.state.trust_region_active, "trust_region_center": self.state.trust_region_center, "trust_region_radius": self.state.trust_region_radius, "best_objective_value": self.state.best_objective_value, }, "parameters": [ {"name": p.name, "min": p.min_value, "max": p.max_value, "step": p.step, "unit": p.unit} for p in self.parameters ], "points": [ { "id": p.id, "params": p.params, "status": p.status, "metrics": p.metrics, "batch_id": p.batch_id, "feasible": p.feasibility_report.get("feasible") if p.feasibility_report else None, } for p in self.state.points ], }