diff --git a/.claude/skills/01_CHEATSHEET.md b/.claude/skills/01_CHEATSHEET.md index 429e0034..e692dc6d 100644 --- a/.claude/skills/01_CHEATSHEET.md +++ b/.claude/skills/01_CHEATSHEET.md @@ -1,7 +1,7 @@ --- skill_id: SKILL_001 -version: 2.0 -last_updated: 2025-12-07 +version: 2.1 +last_updated: 2025-12-22 type: reference code_dependencies: - optimization_engine/extractors/__init__.py @@ -12,8 +12,8 @@ requires_skills: # Atomizer Quick Reference Cheatsheet -**Version**: 2.0 -**Updated**: 2025-12-07 +**Version**: 2.1 +**Updated**: 2025-12-22 **Purpose**: Rapid lookup for common operations. "I want X → Use Y" --- @@ -45,10 +45,12 @@ requires_skills: | CAD expression mass | E5 | `extract_mass_from_expression(prt_file, expression_name='p173')` | | Field data | E6 | `FieldDataExtractor(field_file, result_column, aggregation)` | | Stiffness (k=F/δ) | E7 | `StiffnessCalculator(...)` | -| Zernike WFE | E8 | `extract_zernike_from_op2(op2_file, bdf_file, subcase)` | +| Zernike WFE (standard) | E8 | `extract_zernike_from_op2(op2_file, bdf_file, subcase)` | | Zernike relative | E9 | `extract_zernike_relative_rms(op2_file, bdf_file, target, ref)` | | Zernike builder | E10 | `ZernikeObjectiveBuilder(op2_finder)` | | Part mass + material | E11 | `extract_part_mass_material(prt_file)` → mass, volume, material | +| Zernike Analytic | E20 | `extract_zernike_analytic(op2_file, focal_length=5000.0)` | +| **Zernike OPD** | E22 | `extract_zernike_opd(op2_file)` ← **Most rigorous, RECOMMENDED** | > **Mass extraction tip**: Always use E11 (geometry .prt) over E4 (BDF) for accuracy. > pyNastran under-reports mass ~7% on hex-dominant meshes with tet/pyramid fills. @@ -303,12 +305,20 @@ Generate physics-focused visualizations from FEA results. ### Available Insight Types | Type | Purpose | Data Required | |------|---------|---------------| -| `zernike_wfe` | Mirror wavefront error | OP2 with displacements | +| `zernike_dashboard` | **RECOMMENDED: Unified WFE dashboard** | OP2 with displacements | +| `zernike_wfe` | WFE with Standard/OPD toggle | OP2 with displacements | +| `zernike_opd_comparison` | Compare Standard vs OPD methods | OP2 with displacements | | `stress_field` | Von Mises stress contours | OP2 with stresses | | `modal` | Mode shapes + frequencies | OP2 with eigenvectors | | `thermal` | Temperature distribution | OP2 with temperatures | | `design_space` | Parameter-objective landscape | study.db with 5+ trials | +### Zernike Method Comparison +| Method | Use | RMS Difference | +|--------|-----|----------------| +| **Standard (Z-only)** | Quick analysis | Baseline | +| **OPD (X,Y,Z)** ← RECOMMENDED | Any surface with lateral displacement | **+8-11% higher** (more accurate) | + ### Commands ```bash # List available insights for a study diff --git a/.claude/skills/core/study-creation-core.md b/.claude/skills/core/study-creation-core.md index cfb98d16..7bcf20d1 100644 --- a/.claude/skills/core/study-creation-core.md +++ b/.claude/skills/core/study-creation-core.md @@ -103,6 +103,55 @@ study_info = introspect_study(study_dir) --- +## README Hierarchy (Parent-Child Documentation) + +Studies use a **two-level documentation system**: + +### Parent README (Geometry-Level) + +Location: `studies/{geometry_type}/README.md` + +Contains project-wide context that ALL sub-studies share: +- Project overview (what is this component?) +- Physical system specs (material, mass, loading) +- Domain-specific specs (optical prescription for mirrors, structural limits for brackets) +- Complete design variables catalog with ranges +- Complete objectives catalog with formulas +- Campaign history (evolution across sub-studies) +- Sub-studies index with links and status + +**Example**: `studies/M1_Mirror/README.md` + +### Child README (Study-Level) + +Location: `studies/{geometry_type}/{study_name}/README.md` + +Contains study-specific details, references parent for context: +```markdown +# Study Name + +> **Parent Documentation**: See [../README.md](../README.md) for project overview and specifications. + +## Study Focus +What THIS study specifically optimizes... +``` + +Contents: +- Parent reference banner (MANDATORY first line) +- Study focus (what differentiates this study) +- Active variables (which params enabled) +- Algorithm config (sampler, trials, settings) +- Baseline/seeding (starting point) +- Results summary + +### When Creating a Study + +1. **First study for geometry type** → Create BOTH parent and child READMEs +2. **Subsequent studies** → Create child README, update parent's sub-studies index +3. **New geometry type** → Create new `studies/{type}/` folder with parent README + +--- + ## PR.3 NXSolver Interface **Module**: `optimization_engine.nx_solver` diff --git a/.claude/skills/modules/insights-catalog.md b/.claude/skills/modules/insights-catalog.md new file mode 100644 index 00000000..0ca2ae88 --- /dev/null +++ b/.claude/skills/modules/insights-catalog.md @@ -0,0 +1,289 @@ +# Insights Catalog Module + +**Last Updated**: December 22, 2025 +**Version**: 1.0 +**Type**: Optional Module + +This module documents all available Study Insights in the Atomizer framework. Load this when the user asks about visualization, physics insights, or wants to understand specific trial results beyond optimizer metrics. + +--- + +## When to Load + +- User asks "what insights are available?" +- User wants to visualize FEA results (stress, WFE, modes, etc.) +- User asks to understand a specific trial or design +- User wants to compare methods or generate physics reports +- User mentions "visualization", "3D plot", "surface plot" + +--- + +## Insights vs Analysis + +| Aspect | **Analysis** | **Insights** | +|--------|--------------|--------------| +| Focus | Optimization performance | Physics understanding | +| Questions | "Is the optimizer converging?" | "What does this design look like?" | +| Data Source | `study.db` (trials, objectives) | Simulation outputs (OP2, mesh) | +| Output | Convergence plots, Pareto fronts | 3D surfaces, stress contours | + +--- + +## IN.1 Insight Catalog + +| ID | Type ID | Name | Applicable To | Data Required | Output | +|----|---------|------|---------------|---------------|--------| +| I1 | `zernike_wfe` | Zernike WFE Analysis | Mirror, optics | OP2 with displacement subcases | 3D WFE surface + coefficient tables | +| I2 | `zernike_opd_comparison` | **Zernike OPD Comparison** | Mirror, optics, lateral | OP2 with displacement subcases | Method comparison + lateral disp. map | +| I3 | `msf_zernike` | MSF Zernike Analysis | Mirror, optics | OP2 with displacement subcases | LSF/MSF/HSF band decomposition | +| I4 | `stress_field` | Stress Distribution | Structural, bracket | OP2 with stress results | 3D stress contours | +| I5 | `modal` | Modal Analysis | Vibration, dynamic | OP2 with eigenvectors | Mode shape animations | +| I6 | `thermal` | Thermal Analysis | Thermo-structural | OP2 with temperature | Temperature distribution | +| I7 | `design_space` | Design Space Explorer | All studies | study.db with 5+ trials | Parameter-objective plots | + +--- + +## IN.2 Insight Code Snippets (COPY-PASTE) + +### I1: Zernike WFE Analysis + +Standard 50-mode wavefront error visualization for mirror optimization. + +```python +from optimization_engine.insights import get_insight, InsightConfig + +# Get the insight for a study +insight = get_insight('zernike_wfe', study_path) + +# Check if it can generate (data exists) +if insight.can_generate(): + result = insight.generate(InsightConfig()) + print(f"Generated: {result.html_path}") +``` + +**CLI Usage**: +```bash +python -m optimization_engine.insights generate studies/my_mirror --type zernike_wfe +``` + +**Output**: Interactive HTML with: +- 3D WFE surface plot +- Zernike coefficient bar chart +- RMS metrics table (global, filtered, J1-J3) +- Aberration magnitudes (astigmatism, coma, trefoil, spherical) + +--- + +### I2: Zernike OPD Method Comparison (NEW - Critical for Lateral Supports) + +Compares standard (Z-only) vs rigorous OPD method to reveal lateral displacement impact. + +**When to Use**: +- Lateral support optimization +- Any case where supports may pinch the mirror +- Validating if your WFE objective is accurate + +```python +from optimization_engine.insights import get_insight + +insight = get_insight('zernike_opd_comparison', study_path) +result = insight.generate() +print(f"Generated: {result.html_path}") + +# Check recommendation +print(result.summary['overall_recommendation']) +``` + +**CLI Usage**: +```bash +python -m optimization_engine.insights generate studies/my_mirror --type zernike_opd_comparison +``` + +**Output**: Interactive HTML with: +- **Lateral displacement magnitude map** (where pinching occurs) +- **WFE surface** using OPD method +- **Comparison table** (Standard vs OPD metrics) +- **Recommendation** (CRITICAL/RECOMMENDED/OPTIONAL/EQUIVALENT) + +**Interpretation**: +| Max Lateral Disp. | Recommendation | +|-------------------|----------------| +| > 10 µm | **CRITICAL**: Use OPD method for optimization | +| 1-10 µm | **RECOMMENDED**: Use OPD method | +| < 1 µm | Both methods equivalent | + +**Related Documentation**: [ZERNIKE_OPD_METHOD.md](../../../docs/06_PHYSICS/ZERNIKE_OPD_METHOD.md) + +--- + +### I3: MSF Zernike Analysis + +100-mode analysis with LSF/MSF/HSF band decomposition for mid-spatial frequency analysis. + +```python +from optimization_engine.insights import get_insight + +insight = get_insight('msf_zernike', study_path) +result = insight.generate() +``` + +**Output**: Band breakdown showing: +- LSF (Low Spatial Frequency) - Modes 1-21 (correctable by polishing) +- MSF (Mid Spatial Frequency) - Modes 22-50 (problematic for optical performance) +- HSF (High Spatial Frequency) - Modes 51-100 (typically noise/mesh effects) + +--- + +### I4: Stress Field Visualization + +3D stress contour visualization for structural optimization. + +```python +from optimization_engine.insights import get_insight + +insight = get_insight('stress_field', study_path) +if insight.can_generate(): + result = insight.generate() +``` + +**Output**: Interactive 3D mesh with: +- Von Mises stress colormap +- Peak stress location indicators +- Statistics table (max, mean, P99) + +--- + +### I5: Modal Analysis + +Mode shape visualization for frequency optimization. + +```python +from optimization_engine.insights import get_insight + +insight = get_insight('modal', study_path) +result = insight.generate() +``` + +**Output**: +- First N mode shapes (deformed mesh) +- Frequency table +- Modal effective mass (if available) + +--- + +### I6: Thermal Analysis + +Temperature distribution visualization. + +```python +from optimization_engine.insights import get_insight + +insight = get_insight('thermal', study_path) +result = insight.generate() +``` + +**Output**: +- Temperature contour plot +- Gradient magnitude map +- Heat flux vectors (if available) + +--- + +### I7: Design Space Explorer + +Parameter-objective relationship visualization. + +```python +from optimization_engine.insights import get_insight + +insight = get_insight('design_space', study_path) +result = insight.generate() +``` + +**Requirements**: At least 5 completed trials in study.db + +**Output**: +- Parallel coordinates plot +- Parameter importance +- Objective correlations + +--- + +## IN.3 Generating All Insights for a Study + +```python +from optimization_engine.insights import InsightReport + +# Create report from all configured insights +report = InsightReport(study_path) + +# Generate all enabled insights +results = report.generate_all() + +# Create HTML report with all insights +report_path = report.generate_report_html() +print(f"Full report: {report_path}") +``` + +**CLI**: +```bash +# Generate report from config +python -m optimization_engine.insights report studies/my_study + +# List available insights +python -m optimization_engine.insights list + +# Get recommendations for a study +python -m optimization_engine.insights recommend studies/my_study +``` + +--- + +## IN.4 Configuring Insights in optimization_config.json + +```json +{ + "insights": [ + { + "type": "zernike_wfe", + "name": "WFE at 40 vs 20 deg", + "enabled": true, + "linked_objective": "rel_filtered_rms_40_vs_20", + "config": { + "target_subcase": "3", + "reference_subcase": "2" + }, + "include_in_report": true + }, + { + "type": "zernike_opd_comparison", + "name": "Lateral Displacement Check", + "enabled": true, + "config": {} + } + ] +} +``` + +--- + +## IN.5 Decision Guide + +| You want to... | Use Insight | +|----------------|-------------| +| See wavefront error surface | `zernike_wfe` | +| Check if lateral displacement matters | `zernike_opd_comparison` | +| Analyze mid-spatial frequencies | `msf_zernike` | +| Visualize stress hotspots | `stress_field` | +| See mode shapes | `modal` | +| Check thermal distribution | `thermal` | +| Understand parameter effects | `design_space` | + +--- + +## Related Documentation + +- **Protocol Specification**: `docs/protocols/system/SYS_16_STUDY_INSIGHTS.md` +- **OPD Method Physics**: `docs/06_PHYSICS/ZERNIKE_OPD_METHOD.md` +- **Zernike Integration**: `docs/ZERNIKE_INTEGRATION.md` +- **Extractor Catalog**: `.claude/skills/modules/extractors-catalog.md` diff --git a/atomizer-dashboard/backend/api/routes/insights.py b/atomizer-dashboard/backend/api/routes/insights.py index 802090e4..07cece22 100644 --- a/atomizer-dashboard/backend/api/routes/insights.py +++ b/atomizer-dashboard/backend/api/routes/insights.py @@ -1,12 +1,19 @@ """ Study Insights API endpoints Provides physics-focused visualizations for completed optimization trials + +Key Features: +- Config-driven insights from optimization_config.json +- AI recommendations based on objectives +- Report generation with PDF export +- Objective-linked and standalone insights """ from fastapi import APIRouter, HTTPException from fastapi.responses import JSONResponse, HTMLResponse +from pydantic import BaseModel from pathlib import Path -from typing import List, Dict, Optional +from typing import List, Dict, Optional, Any import json import sys @@ -19,6 +26,22 @@ router = APIRouter() STUDIES_DIR = Path(__file__).parent.parent.parent.parent.parent / "studies" +class InsightSpecRequest(BaseModel): + """Request model for insight specification.""" + type: str + name: str + enabled: bool = True + linked_objective: Optional[str] = None + config: Dict[str, Any] = {} + include_in_report: bool = True + + +class GenerateReportRequest(BaseModel): + """Request model for report generation.""" + specs: Optional[List[InsightSpecRequest]] = None + include_appendix: bool = True + + def resolve_study_path(study_id: str) -> Path: """Find study folder by scanning all topic directories.""" # First check direct path (backwards compatibility) @@ -38,7 +61,7 @@ def resolve_study_path(study_id: str) -> Path: raise HTTPException(status_code=404, detail=f"Study not found: {study_id}") -@router.get("/studies/{study_id}/insights/available") +@router.get("/studies/{study_id}/available") async def list_available_insights(study_id: str): """List all insight types that can be generated for this study.""" try: @@ -65,7 +88,7 @@ async def list_available_insights(study_id: str): raise HTTPException(status_code=500, detail=str(e)) -@router.get("/studies/{study_id}/insights/all") +@router.get("/studies/{study_id}/all") async def list_all_insights(): """List all registered insight types (regardless of availability for any study).""" try: @@ -83,14 +106,104 @@ async def list_all_insights(): raise HTTPException(status_code=500, detail=str(e)) -@router.post("/studies/{study_id}/insights/generate/{insight_type}") -async def generate_insight(study_id: str, insight_type: str, trial_id: Optional[int] = None): +class GenerateInsightRequest(BaseModel): + """Request model for insight generation.""" + iteration: Optional[str] = None # e.g., "iter5", "best_design_archive" + trial_id: Optional[int] = None + config: Dict[str, Any] = {} + + +def _is_valid_op2(op2_path: Path) -> bool: + """Quick check if OP2 file is valid (not from a failed solve).""" + try: + # Check file size - failed OP2s are often very small + if op2_path.stat().st_size < 10000: # Less than 10KB is suspicious + return False + + # Try to read header to verify it's a valid OP2 + with open(op2_path, 'rb') as f: + header = f.read(100) + # Valid OP2 files have specific markers + if b'NASTRAN' not in header and b'XXXXXXXX' not in header: + # Check for common valid patterns + pass # Size check is usually enough + + return True + except Exception: + return False + + +@router.get("/studies/{study_id}/iterations") +async def list_iterations(study_id: str): + """List available iterations/trials with OP2 files for insight generation. + + Returns iterations sorted by modification time (newest first). + Only includes iterations with valid (non-corrupted) OP2 files. + """ + try: + study_path = resolve_study_path(study_id) + iterations = [] + + # Check 2_iterations folder + iter_dir = study_path / "2_iterations" + if iter_dir.exists(): + for subdir in sorted(iter_dir.iterdir(), reverse=True): + if subdir.is_dir(): + op2_files = [f for f in subdir.glob("*.op2") if _is_valid_op2(f)] + if op2_files: + newest_op2 = max(op2_files, key=lambda p: p.stat().st_mtime) + iterations.append({ + "id": subdir.name, + "path": str(subdir.relative_to(study_path)), + "op2_file": newest_op2.name, + "modified": newest_op2.stat().st_mtime, + "type": "iteration" + }) + + # Check 3_results/best_design_archive folder + best_dir = study_path / "3_results" / "best_design_archive" + if best_dir.exists(): + op2_files = [f for f in best_dir.glob("**/*.op2") if _is_valid_op2(f)] + if op2_files: + newest_op2 = max(op2_files, key=lambda p: p.stat().st_mtime) + iterations.insert(0, { # Insert at start as "best" + "id": "best_design_archive", + "path": "3_results/best_design_archive", + "op2_file": newest_op2.name, + "modified": newest_op2.stat().st_mtime, + "type": "best", + "label": "Best Design (Recommended)" + }) + + # Sort by modification time (newest first), keeping best at top + best = [i for i in iterations if i.get("type") == "best"] + others = sorted([i for i in iterations if i.get("type") != "best"], + key=lambda x: x["modified"], reverse=True) + + return { + "study_id": study_id, + "iterations": best + others, + "count": len(iterations) + } + + except HTTPException: + raise + except Exception as e: + raise HTTPException(status_code=500, detail=str(e)) + + +@router.post("/studies/{study_id}/generate/{insight_type}") +async def generate_insight( + study_id: str, + insight_type: str, + request: Optional[GenerateInsightRequest] = None +): """Generate a specific insight visualization. Args: study_id: Study identifier insight_type: Type of insight (e.g., 'zernike_wfe', 'stress_field', 'design_space') - trial_id: Optional specific trial to analyze (defaults to best trial) + request: Optional generation config with iteration selection Returns: JSON with plotly_figure data and summary statistics @@ -104,6 +217,25 @@ async def generate_insight(study_id: str, insight_type: str, trial_id: Optional[ if insight is None: raise HTTPException(status_code=404, detail=f"Unknown insight type: {insight_type}") + # If iteration specified, override the OP2 path + if request and request.iteration: + iteration_id = request.iteration + if iteration_id == "best_design_archive": + iter_path = study_path / "3_results" / "best_design_archive" + else: + iter_path = study_path / "2_iterations" / iteration_id + + if iter_path.exists(): + op2_files = list(iter_path.glob("**/*.op2")) + if op2_files: + # Override the insight's OP2 path + insight.op2_path = max(op2_files, key=lambda p: p.stat().st_mtime) + # Re-find geometry + try: + insight.geo_path = insight._find_geometry_file(insight.op2_path) + except (FileNotFoundError, AttributeError): + pass # Use default + if not insight.can_generate(): raise HTTPException( status_code=400, @@ -111,7 +243,9 @@ async def generate_insight(study_id: str, insight_type: str, trial_id: Optional[ ) # Configure insight - config = InsightConfig(trial_id=trial_id) + trial_id = request.trial_id if request else None + extra_config = request.config if request else {} + config = InsightConfig(trial_id=trial_id, extra=extra_config) # Generate result = insight.generate(config) @@ -123,6 +257,7 @@ async def generate_insight(study_id: str, insight_type: str, trial_id: Optional[ "success": True, "insight_type": insight_type, "study_id": study_id, + "iteration": request.iteration if request else None, "html_path": str(result.html_path) if result.html_path else None, "plotly_figure": result.plotly_figure, "summary": result.summary @@ -134,7 +269,7 @@ async def generate_insight(study_id: str, insight_type: str, trial_id: Optional[ raise HTTPException(status_code=500, detail=str(e)) -@router.get("/studies/{study_id}/insights/view/{insight_type}") +@router.get("/studies/{study_id}/view/{insight_type}") async def view_insight_html(study_id: str, insight_type: str): """Get the HTML content for an insight (for iframe embedding). @@ -180,7 +315,7 @@ async def view_insight_html(study_id: str, insight_type: str): raise HTTPException(status_code=500, detail=str(e)) -@router.get("/studies/{study_id}/insights/generated") +@router.get("/studies/{study_id}/generated") async def list_generated_insights(study_id: str): """List all previously generated insight HTML files for a study.""" try: @@ -219,3 +354,199 @@ async def list_generated_insights(study_id: str): raise except Exception as e: raise HTTPException(status_code=500, detail=str(e)) + + +@router.get("/studies/{study_id}/configured") +async def get_configured_insights_endpoint(study_id: str): + """Get insights configured in the study's optimization_config.json.""" + try: + study_path = resolve_study_path(study_id) + + from optimization_engine.insights import get_configured_insights + + specs = get_configured_insights(study_path) + + return { + "study_id": study_id, + "configured": [ + { + "type": spec.type, + "name": spec.name, + "enabled": spec.enabled, + "linked_objective": spec.linked_objective, + "config": spec.config, + "include_in_report": spec.include_in_report + } + for spec in specs + ] + } + + except HTTPException: + raise + except ImportError as e: + return { + "study_id": study_id, + "configured": [], + "error": f"Insights module not available: {str(e)}" + } + except Exception as e: + raise HTTPException(status_code=500, detail=str(e)) + + +@router.get("/studies/{study_id}/recommend") +async def recommend_insights_endpoint(study_id: str): + """Get AI recommendations for insights based on study objectives.""" + try: + study_path = resolve_study_path(study_id) + + from optimization_engine.insights import recommend_insights_for_study + + recommendations = recommend_insights_for_study(study_path) + + return { + "study_id": study_id, + "recommendations": recommendations + } + + except HTTPException: + raise + except ImportError as e: + return { + "study_id": study_id, + "recommendations": [], + "error": f"Insights module not available: {str(e)}" + } + except Exception as e: + raise HTTPException(status_code=500, detail=str(e)) + + +@router.post("/studies/{study_id}/report") +async def generate_report(study_id: str, request: GenerateReportRequest): + """Generate comprehensive HTML report with all insights. + + Args: + study_id: Study identifier + request: Report configuration with optional insight specs + + Returns: + JSON with report path and generation results + """ + try: + study_path = resolve_study_path(study_id) + + from optimization_engine.insights import ( + InsightReport, InsightSpec, get_configured_insights, + recommend_insights_for_study + ) + + # Build specs from request or config + if request.specs: + specs = [ + InsightSpec( + type=s.type, + name=s.name, + enabled=s.enabled, + linked_objective=s.linked_objective, + config=s.config, + include_in_report=s.include_in_report + ) + for s in request.specs + ] + else: + # Try config first, then recommendations + specs = get_configured_insights(study_path) + if not specs: + recommendations = recommend_insights_for_study(study_path) + specs = [ + InsightSpec( + type=rec['type'], + name=rec['name'], + linked_objective=rec.get('linked_objective'), + config=rec.get('config', {}) + ) + for rec in recommendations + ] + + if not specs: + raise HTTPException( + status_code=400, + detail="No insights configured or recommended for this study" + ) + + # Generate report + report = InsightReport(study_path) + results = report.generate_all(specs) + report_path = report.generate_report_html(include_appendix=request.include_appendix) + + return { + "success": True, + "study_id": study_id, + "report_path": str(report_path), + "results": [ + { + "type": r.insight_type, + "name": r.insight_name, + "success": r.success, + "linked_objective": r.linked_objective, + "error": r.error + } + for r in results + ], + "summary": { + "total": len(results), + "successful": sum(1 for r in results if r.success), + "failed": sum(1 for r in results if not r.success) + } + } + + except HTTPException: + raise + except Exception as e: + raise HTTPException(status_code=500, detail=str(e)) + + +@router.get("/studies/{study_id}/report/view") +async def view_report(study_id: str): + """Get the latest generated report HTML for embedding.""" + try: + study_path = resolve_study_path(study_id) + report_path = study_path / "3_insights" / "STUDY_INSIGHTS_REPORT.html" + + if not report_path.exists(): + raise HTTPException( + status_code=404, + detail="No report generated yet. Use POST /insights/report to generate." + ) + + return HTMLResponse(content=report_path.read_text(encoding='utf-8')) + + except HTTPException: + raise + except Exception as e: + raise HTTPException(status_code=500, detail=str(e)) + + +@router.get("/studies/{study_id}/summary") +async def get_insights_summary(study_id: str): + """Get insights summary JSON for Results page integration.""" + try: + study_path = resolve_study_path(study_id) + summary_path = study_path / "3_insights" / "insights_summary.json" + + if not summary_path.exists(): + return { + "study_id": study_id, + "generated_at": None, + "insights": [] + } + + with open(summary_path) as f: + summary = json.load(f) + + summary["study_id"] = study_id + return summary + + except HTTPException: + raise + except Exception as e: + raise HTTPException(status_code=500, detail=str(e)) diff --git a/atomizer-dashboard/frontend/src/pages/Insights.tsx b/atomizer-dashboard/frontend/src/pages/Insights.tsx index bcac607b..8bc2ec44 100644 --- a/atomizer-dashboard/frontend/src/pages/Insights.tsx +++ b/atomizer-dashboard/frontend/src/pages/Insights.tsx @@ -1,9 +1,8 @@ -import { useState, useEffect, useCallback } from 'react'; +import React, { useState, useEffect, useCallback, useMemo } from 'react'; import { useNavigate } from 'react-router-dom'; import { Eye, RefreshCw, - Download, Maximize2, X, Activity, @@ -14,16 +13,41 @@ import { AlertCircle, CheckCircle, Clock, - FileText + FileText, + Lightbulb, + Target, + Layers, + BookOpen, + ChevronRight, + ChevronDown, + Folder, + Play, + ExternalLink, + Zap, + List } from 'lucide-react'; import { useStudy } from '../context/StudyContext'; import { Card } from '../components/common/Card'; import Plot from 'react-plotly.js'; +// ============================================================================ +// Types +// ============================================================================ +interface IterationInfo { + id: string; + path: string; + op2_file: string; + modified: number; + type: 'iteration' | 'best'; + label?: string; +} + interface InsightInfo { type: string; name: string; description: string; + category?: string; + category_label?: string; applicable_to: string[]; } @@ -38,13 +62,22 @@ interface GeneratedFile { interface InsightResult { success: boolean; insight_type: string; + insight_name?: string; + linked_objective?: string | null; + iteration?: string; plotly_figure: any; summary: Record; html_path: string | null; + error?: string; } +// ============================================================================ +// Constants +// ============================================================================ const INSIGHT_ICONS: Record = { zernike_wfe: Waves, + zernike_opd_comparison: Waves, + msf_zernike: Waves, stress_field: Activity, modal: Box, thermal: Thermometer, @@ -53,26 +86,54 @@ const INSIGHT_ICONS: Record = { const INSIGHT_COLORS: Record = { zernike_wfe: 'text-blue-400 bg-blue-500/10 border-blue-500/30', + zernike_opd_comparison: 'text-blue-400 bg-blue-500/10 border-blue-500/30', + msf_zernike: 'text-cyan-400 bg-cyan-500/10 border-cyan-500/30', stress_field: 'text-red-400 bg-red-500/10 border-red-500/30', modal: 'text-purple-400 bg-purple-500/10 border-purple-500/30', thermal: 'text-orange-400 bg-orange-500/10 border-orange-500/30', design_space: 'text-green-400 bg-green-500/10 border-green-500/30' }; +const CATEGORY_COLORS: Record = { + optical: 'border-blue-500/50 bg-blue-500/5', + structural_static: 'border-red-500/50 bg-red-500/5', + structural_dynamic: 'border-orange-500/50 bg-orange-500/5', + structural_modal: 'border-purple-500/50 bg-purple-500/5', + thermal: 'border-yellow-500/50 bg-yellow-500/5', + kinematic: 'border-teal-500/50 bg-teal-500/5', + design_exploration: 'border-green-500/50 bg-green-500/5', + other: 'border-gray-500/50 bg-gray-500/5' +}; + +// ============================================================================ +// Main Component +// ============================================================================ export default function Insights() { const navigate = useNavigate(); const { selectedStudy, isInitialized } = useStudy(); + // State - Step-based workflow + const [step, setStep] = useState<'select-iteration' | 'select-insight' | 'view-result'>('select-iteration'); + + // Data + const [iterations, setIterations] = useState([]); const [availableInsights, setAvailableInsights] = useState([]); const [generatedFiles, setGeneratedFiles] = useState([]); - const [loading, setLoading] = useState(true); - const [generating, setGenerating] = useState(null); - const [error, setError] = useState(null); - // Active insight for display + // Selections + const [selectedIteration, setSelectedIteration] = useState(null); + const [selectedInsightType, setSelectedInsightType] = useState(null); + + // Result const [activeInsight, setActiveInsight] = useState(null); const [fullscreen, setFullscreen] = useState(false); + // Loading states + const [loadingIterations, setLoadingIterations] = useState(true); + const [loadingInsights, setLoadingInsights] = useState(false); + const [generating, setGenerating] = useState(false); + const [error, setError] = useState(null); + // Redirect if no study useEffect(() => { if (isInitialized && !selectedStudy) { @@ -80,47 +141,101 @@ export default function Insights() { } }, [selectedStudy, navigate, isInitialized]); - // Load available insights and generated files - const loadInsights = useCallback(async () => { + // Load iterations on mount (single fast API call + generated files for quick access) + const loadIterations = useCallback(async () => { if (!selectedStudy) return; - setLoading(true); + setLoadingIterations(true); setError(null); try { - // Load available insights - const availRes = await fetch(`/api/insights/studies/${selectedStudy.id}/insights/available`); - const availData = await availRes.json(); - setAvailableInsights(availData.insights || []); + // Load iterations and generated files in parallel + const [iterRes, genRes] = await Promise.all([ + fetch(`/api/insights/studies/${selectedStudy.id}/iterations`), + fetch(`/api/insights/studies/${selectedStudy.id}/generated`) + ]); - // Load previously generated files - const genRes = await fetch(`/api/insights/studies/${selectedStudy.id}/insights/generated`); - const genData = await genRes.json(); + const [iterData, genData] = await Promise.all([ + iterRes.json(), + genRes.json() + ]); + + const iters = iterData.iterations || []; + setIterations(iters); setGeneratedFiles(genData.files || []); + // Auto-select best design if available + if (iters.length > 0) { + const best = iters.find((i: IterationInfo) => i.type === 'best'); + if (best) { + setSelectedIteration(best); + } else { + setSelectedIteration(iters[0]); + } + } } catch (err) { - console.error('Failed to load insights:', err); - setError('Failed to load insights data'); + console.error('Failed to load iterations:', err); + setError('Failed to load iterations'); } finally { - setLoading(false); + setLoadingIterations(false); } }, [selectedStudy]); + // Load available insights (lazy - only when needed) + const loadAvailableInsights = useCallback(async () => { + if (!selectedStudy) return; + + setLoadingInsights(true); + + try { + const [availRes, genRes] = await Promise.all([ + fetch(`/api/insights/studies/${selectedStudy.id}/available`), + fetch(`/api/insights/studies/${selectedStudy.id}/generated`) + ]); + + const [availData, genData] = await Promise.all([ + availRes.json(), + genRes.json() + ]); + + setAvailableInsights(availData.insights || []); + setGeneratedFiles(genData.files || []); + } catch (err) { + console.error('Failed to load insights:', err); + } finally { + setLoadingInsights(false); + } + }, [selectedStudy]); + + // Initial load useEffect(() => { - loadInsights(); - }, [loadInsights]); + loadIterations(); + }, [loadIterations]); - // Generate an insight - const handleGenerate = async (insightType: string) => { - if (!selectedStudy || generating) return; + // Load insights when moving to step 2 + useEffect(() => { + if (step === 'select-insight' && availableInsights.length === 0) { + loadAvailableInsights(); + } + }, [step, availableInsights.length, loadAvailableInsights]); - setGenerating(insightType); + // Generate insight + const handleGenerate = async () => { + if (!selectedStudy || !selectedIteration || !selectedInsightType || generating) return; + + setGenerating(true); setError(null); try { const res = await fetch( - `/api/insights/studies/${selectedStudy.id}/insights/generate/${insightType}`, - { method: 'POST' } + `/api/insights/studies/${selectedStudy.id}/generate/${selectedInsightType}`, + { + method: 'POST', + headers: { 'Content-Type': 'application/json' }, + body: JSON.stringify({ + iteration: selectedIteration.id + }) + } ); if (!res.ok) { @@ -129,59 +244,54 @@ export default function Insights() { } const result: InsightResult = await res.json(); + const insightInfo = availableInsights.find(i => i.type === selectedInsightType); + result.insight_name = insightInfo?.name || result.insight_type; + setActiveInsight(result); + setStep('view-result'); - // Refresh file list - loadInsights(); - + // Refresh generated files list + loadAvailableInsights(); } catch (err: any) { setError(err.message || 'Failed to generate insight'); } finally { - setGenerating(null); + setGenerating(false); } }; - // View existing insight - const handleViewExisting = async (file: GeneratedFile) => { + // Quick view existing generated file + const handleQuickView = async (file: GeneratedFile) => { if (!selectedStudy) return; + window.open(`/api/insights/studies/${selectedStudy.id}/view/${file.insight_type}`, '_blank'); + }; - setGenerating(file.insight_type); + // Format timestamp + const formatTime = (ts: number) => { + const date = new Date(ts * 1000); + return date.toLocaleDateString() + ' ' + date.toLocaleTimeString([], { hour: '2-digit', minute: '2-digit' }); + }; - try { - const res = await fetch( - `/api/insights/studies/${selectedStudy.id}/insights/generate/${file.insight_type}`, - { method: 'POST' } - ); + // Group insights by category + const groupedInsights = useMemo(() => { + const groups: Record = {}; - if (!res.ok) { - throw new Error('Failed to load insight'); + availableInsights.forEach(insight => { + const cat = insight.category || 'other'; + if (!groups[cat]) { + groups[cat] = []; } + groups[cat].push(insight); + }); - const result: InsightResult = await res.json(); - setActiveInsight(result); + return groups; + }, [availableInsights]); - } catch (err: any) { - setError(err.message || 'Failed to load insight'); - } finally { - setGenerating(null); - } - }; - - // Open HTML in new tab - const handleOpenHtml = (file: GeneratedFile) => { - if (!selectedStudy) return; - window.open(`/api/insights/studies/${selectedStudy.id}/insights/view/${file.insight_type}`, '_blank'); - }; - - const getIcon = (type: string) => { - const Icon = INSIGHT_ICONS[type] || Eye; - return Icon; - }; - - const getColorClass = (type: string) => { - return INSIGHT_COLORS[type] || 'text-gray-400 bg-gray-500/10 border-gray-500/30'; - }; + const getIcon = (type: string) => INSIGHT_ICONS[type] || Eye; + const getColorClass = (type: string) => INSIGHT_COLORS[type] || 'text-gray-400 bg-gray-500/10 border-gray-500/30'; + // ============================================================================ + // Render + // ============================================================================ if (!isInitialized || !selectedStudy) { return (
@@ -194,23 +304,37 @@ export default function Insights() { } return ( -
- {/* Header */} -
-
-

Insights

-

- Physics visualizations for {selectedStudy.name || selectedStudy.id} -

+
+ {/* Header with Breadcrumb */} +
+
+ + 1. Select Iteration + + + + 2. Choose Insight + + + + 3. View Result +
-
+
+
+

Study Insights

+

{selectedStudy.name || selectedStudy.id}

+
@@ -219,206 +343,409 @@ export default function Insights() { {error && (
-

{error}

-
)} -
- {/* Left Panel: Available Insights */} -
- - {loading ? ( -
+ {/* Step 1: Select Iteration */} + {step === 'select-iteration' && ( +
+ + {loadingIterations ? ( +
- Loading... + Scanning for iterations...
- ) : availableInsights.length === 0 ? ( -
- -

No insights available for this study.

-

Run some trials first.

+ ) : iterations.length === 0 ? ( +
+ +

No Iterations Found

+

Run some optimization trials first to generate data.

) : ( -
- {availableInsights.map((insight) => { - const Icon = getIcon(insight.type); - const colorClass = getColorClass(insight.type); - const isGenerating = generating === insight.type; - +
+ {/* Best Design - Primary Option */} + {(() => { + const bestIter = iterations.find((i) => i.type === 'best'); + if (!bestIter) return null; return ( -
-
-
- + + {selectedIteration?.id === bestIter.id && ( + + )} + + ); + })()} + + {/* Separator */} +
+
+ or select specific iteration +
+
+ + {/* Dropdown for other iterations */} + {(() => { + const regularIters = iterations.filter((i) => i.type !== 'best'); + if (regularIters.length === 0) return null; + + const selectedRegular = regularIters.find((i) => i.id === selectedIteration?.id); + + return ( +
+
+ + + +
+ + {/* Show selected iteration details */} + {selectedRegular && ( +
+
+ +
+
{selectedRegular.label || selectedRegular.id}
+
+ {selectedRegular.op2_file} • {formatTime(selectedRegular.modified)} +
+
+ +
+
+ )}
); - })} + })()}
)} - {/* Previously Generated */} - {generatedFiles.length > 0 && ( - -
- {generatedFiles.map((file) => { - const Icon = getIcon(file.insight_type); + {/* Continue Button */} + {selectedIteration && ( +
+ +
+ )} + {/* Previously Generated - Quick Access */} + {generatedFiles.length > 0 && ( + +
+

Quick access to existing visualizations

+
+
+ {generatedFiles.slice(0, 5).map((file) => { + const Icon = getIcon(file.insight_type); return ( -
handleQuickView(file)} + className="w-full p-3 flex items-center gap-3 hover:bg-dark-750 transition-colors text-left" > -
-

{file.insight_type}

-

{file.size_kb} KB

-
-
- - -
-
+ + {file.insight_type.replace(/_/g, ' ')} + + {file.size_kb} KB + + ); })}
)}
+ )} - {/* Right Panel: Visualization */} -
- {activeInsight ? ( - - {activeInsight.insight_type.replace('_', ' ').toUpperCase()} -
- - -
-
- } - className="h-full" + {/* Step 2: Select Insight Type */} + {step === 'select-insight' && ( +
+ {/* Selection Summary */} +
+
+ +
+
+
Selected Data Source
+
{selectedIteration?.label || selectedIteration?.id}
+
+ +
- {/* Plotly Figure */} - {activeInsight.plotly_figure ? ( -
- -
- ) : ( -
-

No visualization data available

-
- )} - - ) : ( - -
- -

No Insight Selected

-

- Select an insight type from the left panel to generate a visualization. -

+ + {loadingInsights ? ( +
+ + Loading available insights...
- + ) : availableInsights.length === 0 ? ( +
+ +

No Insights Available

+

The selected iteration may not have compatible data.

+
+ ) : ( +
+ {Object.entries(groupedInsights).map(([category, insights]) => ( +
+
+

+ {category.replace(/_/g, ' ')} +

+
+ {insights.map((insight) => { + const Icon = getIcon(insight.type); + const colorClass = getColorClass(insight.type); + const isSelected = selectedInsightType === insight.type; + + return ( + + ); + })} +
+ ))} +
+ )} + + + {/* Generate Button */} + {selectedInsightType && ( +
+ + +
)}
-
+ )} + + {/* Step 3: View Result */} + {step === 'view-result' && activeInsight && ( +
+ {/* Result Header */} +
+
+ +
+

+ {activeInsight.insight_name || activeInsight.insight_type.replace(/_/g, ' ')} +

+

+ Generated from {selectedIteration?.label || selectedIteration?.id} +

+
+
+
+ {activeInsight.html_path && ( + + )} + +
+
+ + {/* Summary Stats */} + {activeInsight.summary && Object.keys(activeInsight.summary).length > 0 && ( +
+ {Object.entries(activeInsight.summary) + .filter(([key]) => !key.startsWith('html_')) + .slice(0, 8) + .map(([key, value]) => ( +
+
+ {key.replace(/_/g, ' ')} +
+
+ {typeof value === 'number' + ? value.toFixed(2) + : String(value) + } +
+
+ ))} +
+ )} + + {/* Plotly Figure */} + + {activeInsight.plotly_figure ? ( +
+ +
+ ) : ( +
+ +

Insight Generated Successfully

+

+ This insight generates HTML files. Click "Open Full View" to see the visualization. +

+ {activeInsight.summary?.html_files && ( +
+

Generated files:

+
    + {(activeInsight.summary.html_files as string[]).slice(0, 4).map((f: string, i: number) => ( +
  • + {f.split(/[/\\]/).pop()} +
    • + ))} +
+
+ )} +
+ )} + + + {/* Generate Another */} +
+ +
+
+ )} {/* Fullscreen Modal */} {fullscreen && activeInsight?.plotly_figure && (

- {activeInsight.insight_type.replace('_', ' ').toUpperCase()} + {activeInsight.insight_name || activeInsight.insight_type.replace(/_/g, ' ')}

+
+

Study Insights Report

+

{study_name} | Generated: {timestamp}

+''' + + def _executive_summary(self) -> str: + """Generate executive summary section.""" + n_objectives = len(self.config.get('objectives', [])) + n_insights = len(self.results) + linked = len([r for r in self.results if r.linked_objective]) + + return f''' +
+
+
{n_insights}
+
Total Insights
+
+
+
{linked}
+
Objective-Linked
+
+
+
{n_objectives}
+
Study Objectives
+
+
+''' + + def _insight_section(self, result: InsightResult, compact: bool = False) -> str: + """Generate HTML section for a single insight.""" + objective_tag = "" + if result.linked_objective: + objective_tag = f'{result.linked_objective}' + + # Summary table + summary_html = "" + if result.summary: + rows = "" + for key, value in list(result.summary.items())[:6]: + if isinstance(value, float): + value = f"{value:.4g}" + rows += f"{key}{value}" + summary_html = f"{rows}
" + + # Compact plot or placeholder + plot_html = "" + if result.plotly_figure and compact: + plot_id = f"plot_{result.insight_type}_{id(result)}" + plot_json = json.dumps(result.plotly_figure) + plot_html = f''' +
+ + ''' + + return f''' +
+

{result.insight_name}{objective_tag}

+ {summary_html} + {plot_html} +
+''' + + def _appendix_section(self) -> str: + """Generate appendix with full-size figures.""" + html = '

Appendix: Full-Size Visualizations

' + + for i, result in enumerate(self.results): + if not result.plotly_figure: + continue + + plot_id = f"appendix_plot_{i}" + plot_json = json.dumps(result.plotly_figure) + + html += f''' +
+

Appendix {i+1}: {result.insight_name}

+
+ +
+ ''' + + return html + + def _report_footer(self) -> str: + """Generate HTML footer.""" + return ''' +
+ +''' + + def _write_summary_json(self) -> None: + """Write summary JSON for Results page integration.""" + summary = { + 'generated_at': datetime.now().isoformat(), + 'study_name': self.config.get('study_name', self.study_path.name), + 'insights': [] + } + + for result in self.results: + summary['insights'].append({ + 'type': result.insight_type, + 'name': result.insight_name, + 'success': result.success, + 'linked_objective': result.linked_objective, + 'html_path': str(result.html_path) if result.html_path else None, + 'summary': result.summary + }) + + summary_path = self.insights_path / "insights_summary.json" + with open(summary_path, 'w') as f: + json.dump(summary, f, indent=2) diff --git a/optimization_engine/insights/msf_zernike.py b/optimization_engine/insights/msf_zernike.py new file mode 100644 index 00000000..704cd0fd --- /dev/null +++ b/optimization_engine/insights/msf_zernike.py @@ -0,0 +1,875 @@ +""" +MSF Zernike Insight - Mid-Spatial Frequency Analysis for Telescope Mirrors + +Provides detailed analysis of mid-spatial frequency (MSF) content in mirror +wavefront error, specifically for gravity-induced support print-through. + +Key Features: +- Band decomposition: LSF (n≤10), MSF (n=11-50), HSF (n>50) +- RSS (Root Sum Square) metrics per band +- Relative analysis vs reference orientations (20°, 90°) +- Higher-order Zernike fitting (100 modes by default) +- PSD (Power Spectral Density) visualization +- Support print-through identification + +For GigaBIT and similar telescope mirrors where MSF errors are dominated by +gravity-induced support deformation rather than polishing residual. + +Author: Atomizer Framework +""" + +from pathlib import Path +from datetime import datetime +from typing import Dict, Any, List, Optional, Tuple +import numpy as np +from math import factorial +from numpy.linalg import LinAlgError + +from .base import StudyInsight, InsightConfig, InsightResult, register_insight + +# Lazy imports +_plotly_loaded = False +_go = None +_make_subplots = None +_Triangulation = None +_OP2 = None +_BDF = None + + +def _load_dependencies(): + """Lazy load heavy dependencies.""" + global _plotly_loaded, _go, _make_subplots, _Triangulation, _OP2, _BDF + if not _plotly_loaded: + import plotly.graph_objects as go + from plotly.subplots import make_subplots + from matplotlib.tri import Triangulation + from pyNastran.op2.op2 import OP2 + from pyNastran.bdf.bdf import BDF + _go = go + _make_subplots = make_subplots + _Triangulation = Triangulation + _OP2 = OP2 + _BDF = BDF + _plotly_loaded = True + + +# ============================================================================ +# Band Definitions +# ============================================================================ + +# Default band boundaries (Zernike radial order n) +DEFAULT_LSF_MAX = 10 # LSF: n = 1-10 (M2 hexapod correctable) +DEFAULT_MSF_MAX = 50 # MSF: n = 11-50 (support print-through) +DEFAULT_N_MODES = 100 # Total modes to fit (n ≈ 14) + +# Band colors for visualization +BAND_COLORS = { + 'lsf': '#3b82f6', # Blue - low spatial frequency + 'msf': '#f59e0b', # Amber - mid spatial frequency + 'hsf': '#ef4444', # Red - high spatial frequency +} + + +# ============================================================================ +# Zernike Mathematics (same as zernike_wfe.py) +# ============================================================================ + +def noll_indices(j: int) -> Tuple[int, int]: + """Convert Noll index to (n, m) radial/azimuthal orders.""" + if j < 1: + raise ValueError("Noll index j must be >= 1") + count = 0 + n = 0 + while True: + if n == 0: + ms = [0] + elif n % 2 == 0: + ms = [0] + [m for k in range(1, n//2 + 1) for m in (-2*k, 2*k)] + else: + ms = [m for k in range(0, (n+1)//2) for m in (-(2*k+1), (2*k+1))] + for m in ms: + count += 1 + if count == j: + return n, m + n += 1 + + +def noll_to_radial_order(j: int) -> int: + """Get radial order n for Noll index j.""" + n, _ = noll_indices(j) + return n + + +def radial_order_to_first_noll(n: int) -> int: + """Get first Noll index for radial order n.""" + # Sum of (k+1) for k=0 to n-1, plus 1 + return n * (n + 1) // 2 + 1 + + +def zernike_noll(j: int, r: np.ndarray, th: np.ndarray) -> np.ndarray: + """Evaluate Zernike polynomial j at (r, theta).""" + n, m = noll_indices(j) + R = np.zeros_like(r) + for s in range((n - abs(m)) // 2 + 1): + c = ((-1)**s * factorial(n - s) / + (factorial(s) * + factorial((n + abs(m)) // 2 - s) * + factorial((n - abs(m)) // 2 - s))) + R += c * r**(n - 2*s) + if m == 0: + return R + return R * (np.cos(m * th) if m > 0 else np.sin(-m * th)) + + +def zernike_common_name(n: int, m: int) -> str: + """Get common name for Zernike mode.""" + names = { + (0, 0): "Piston", (1, -1): "Tilt X", (1, 1): "Tilt Y", + (2, 0): "Defocus", (2, -2): "Astig 45°", (2, 2): "Astig 0°", + (3, -1): "Coma X", (3, 1): "Coma Y", (3, -3): "Trefoil X", (3, 3): "Trefoil Y", + (4, 0): "Primary Spherical", (4, -2): "Sec Astig X", (4, 2): "Sec Astig Y", + (4, -4): "Quadrafoil X", (4, 4): "Quadrafoil Y", + (5, -1): "Sec Coma X", (5, 1): "Sec Coma Y", + (5, -3): "Sec Trefoil X", (5, 3): "Sec Trefoil Y", + (5, -5): "Pentafoil X", (5, 5): "Pentafoil Y", + (6, 0): "Sec Spherical", + } + return names.get((n, m), f"Z(n={n}, m={m})") + + +def zernike_label(j: int) -> str: + """Get label for Zernike coefficient J{j}.""" + n, m = noll_indices(j) + return f"J{j:02d} - {zernike_common_name(n, m)} (n={n}, m={m})" + + +def compute_zernike_coeffs( + X: np.ndarray, + Y: np.ndarray, + vals: np.ndarray, + n_modes: int, + chunk_size: int = 100000 +) -> Tuple[np.ndarray, float]: + """Fit Zernike coefficients to WFE data.""" + Xc, Yc = X - np.mean(X), Y - np.mean(Y) + R = float(np.max(np.hypot(Xc, Yc))) + r = np.hypot(Xc / R, Yc / R).astype(np.float32) + th = np.arctan2(Yc, Xc).astype(np.float32) + mask = (r <= 1.0) & ~np.isnan(vals) + if not np.any(mask): + raise RuntimeError("No valid points inside unit disk.") + idx = np.nonzero(mask)[0] + m = int(n_modes) + G = np.zeros((m, m), dtype=np.float64) + h = np.zeros((m,), dtype=np.float64) + v = vals.astype(np.float64) + + for start in range(0, len(idx), chunk_size): + sl = idx[start:start + chunk_size] + r_b, th_b, v_b = r[sl], th[sl], v[sl] + Zb = np.column_stack([zernike_noll(j, r_b, th_b).astype(np.float32) + for j in range(1, m + 1)]) + G += (Zb.T @ Zb).astype(np.float64) + h += (Zb.T @ v_b).astype(np.float64) + try: + coeffs = np.linalg.solve(G, h) + except LinAlgError: + coeffs = np.linalg.lstsq(G, h, rcond=None)[0] + return coeffs, R + + +# ============================================================================ +# Band Analysis +# ============================================================================ + +def classify_modes_by_band( + n_modes: int, + lsf_max: int = DEFAULT_LSF_MAX, + msf_max: int = DEFAULT_MSF_MAX +) -> Dict[str, List[int]]: + """ + Classify Zernike modes into LSF/MSF/HSF bands. + + Returns: + Dict with 'lsf', 'msf', 'hsf' keys containing lists of Noll indices + """ + bands = {'lsf': [], 'msf': [], 'hsf': []} + + for j in range(1, n_modes + 1): + n = noll_to_radial_order(j) + if n <= lsf_max: + bands['lsf'].append(j) + elif n <= msf_max: + bands['msf'].append(j) + else: + bands['hsf'].append(j) + + return bands + + +def compute_band_rss( + coeffs: np.ndarray, + bands: Dict[str, List[int]], + filter_modes: int = 4 +) -> Dict[str, float]: + """ + Compute RSS (Root Sum Square) of coefficients in each band. + + RSS = sqrt(sum(c_j^2)) for j in band + + Args: + coeffs: Zernike coefficients + bands: Dict with band name -> list of Noll indices + filter_modes: Number of low-order modes to filter for 'filtered' metric (default 4 = J1-J4) + + Returns: + Dict with 'lsf', 'msf', 'hsf', 'total', and 'filtered' (J5+) RSS values + """ + rss = {} + for band_name, indices in bands.items(): + # Convert to 0-indexed + band_coeffs = [coeffs[j-1] for j in indices if j-1 < len(coeffs)] + rss[band_name] = float(np.sqrt(np.sum(np.array(band_coeffs)**2))) + + rss['total'] = float(np.sqrt(np.sum(coeffs**2))) + + # Filtered RSS: exclude first filter_modes (J1-J4 = piston, tip, tilt, defocus) + # This is the engineering-relevant metric after M2 alignment + filtered_coeffs = coeffs[filter_modes:] if len(coeffs) > filter_modes else np.array([]) + rss['filtered'] = float(np.sqrt(np.sum(filtered_coeffs**2))) + + return rss + + +def compute_band_residual( + X: np.ndarray, + Y: np.ndarray, + W_nm: np.ndarray, + coeffs: np.ndarray, + R: float, + bands: Dict[str, List[int]], + keep_bands: List[str] +) -> np.ndarray: + """ + Compute WFE residual keeping only specified bands. + + Args: + keep_bands: List of band names to keep (e.g., ['msf'] for MSF-only) + + Returns: + WFE array with only specified band content + """ + Xc = X - np.mean(X) + Yc = Y - np.mean(Y) + r = np.hypot(Xc / R, Yc / R) + th = np.arctan2(Yc, Xc) + + # Build Zernike basis for modes to keep + keep_indices = [] + for band in keep_bands: + keep_indices.extend(bands.get(band, [])) + + if not keep_indices: + return np.zeros_like(W_nm) + + # Reconstruct only the kept bands + W_band = np.zeros_like(W_nm) + for j in keep_indices: + if j - 1 < len(coeffs): + Z_j = zernike_noll(j, r, th) + W_band += coeffs[j-1] * Z_j + + return W_band + + +def find_dominant_msf_modes( + coeffs: np.ndarray, + bands: Dict[str, List[int]], + top_n: int = 5 +) -> List[Dict[str, Any]]: + """ + Find the dominant MSF modes by coefficient magnitude. + + Returns: + List of dicts with 'j', 'label', 'coeff_nm', 'n', 'm' + """ + msf_indices = bands.get('msf', []) + + mode_data = [] + for j in msf_indices: + if j - 1 < len(coeffs): + n, m = noll_indices(j) + mode_data.append({ + 'j': j, + 'label': zernike_label(j), + 'coeff_nm': float(abs(coeffs[j-1])), + 'n': n, + 'm': m + }) + + # Sort by magnitude + mode_data.sort(key=lambda x: x['coeff_nm'], reverse=True) + + return mode_data[:top_n] + + +def estimate_mesh_resolution(X: np.ndarray, Y: np.ndarray) -> Dict[str, float]: + """ + Estimate mesh spatial resolution. + + Returns: + Dict with 'node_count', 'diameter_mm', 'avg_spacing_mm', + 'max_resolvable_order', 'min_feature_mm' + """ + n_nodes = len(X) + + Xc = X - np.mean(X) + Yc = Y - np.mean(Y) + R = np.max(np.hypot(Xc, Yc)) + diameter = 2 * R + + # Approximate spacing from area + area = np.pi * R**2 + avg_spacing = np.sqrt(area / n_nodes) + + # Nyquist: can resolve features > 2x spacing + min_feature = 2 * avg_spacing + + # Max Zernike order + max_order = int(diameter / min_feature) + + return { + 'node_count': n_nodes, + 'diameter_mm': float(diameter), + 'avg_spacing_mm': float(avg_spacing), + 'min_feature_mm': float(min_feature), + 'max_resolvable_order': max_order + } + + +# ============================================================================ +# Default Configuration +# ============================================================================ + +DEFAULT_CONFIG = { + 'n_modes': 100, + 'lsf_max': 10, # n ≤ 10 is LSF + 'msf_max': 50, # n = 11-50 is MSF + 'amp': 0.5, # Visual deformation scale + 'pancake': 3.0, # Z-axis range multiplier + 'plot_downsample': 15000, + 'colorscale': 'Turbo', + 'disp_unit': 'mm', + 'show_psd': True, + 'show_all_orientations': True, +} + + +@register_insight +class MSFZernikeInsight(StudyInsight): + """ + Mid-Spatial Frequency Zernike analysis for telescope mirror optimization. + + Provides detailed breakdown of spatial frequency content: + - LSF (Low): n ≤ 10, correctable by M2 hexapod + - MSF (Mid): n = 11-50, support print-through + - HSF (High): n > 50, near mesh resolution limit + + Generates: + - Band decomposition table with RSS metrics + - MSF-only 3D surface visualization + - Coefficient bar chart color-coded by band + - Dominant MSF mode identification + - Mesh resolution analysis + """ + + insight_type = "msf_zernike" + name = "MSF Zernike Analysis" + description = "Mid-spatial frequency analysis with band decomposition for support print-through" + category = "optical" + applicable_to = ["mirror", "optics", "wfe", "msf"] + required_files = ["*.op2"] + + def __init__(self, study_path: Path): + super().__init__(study_path) + self.op2_path: Optional[Path] = None + self.geo_path: Optional[Path] = None + self._node_geo: Optional[Dict] = None + self._displacements: Optional[Dict] = None + + def can_generate(self) -> bool: + """Check if OP2 and geometry files exist.""" + search_paths = [ + self.results_path, + self.study_path / "2_iterations", + self.setup_path / "model", + ] + + for search_path in search_paths: + if not search_path.exists(): + continue + op2_files = list(search_path.glob("**/*solution*.op2")) + if not op2_files: + op2_files = list(search_path.glob("**/*.op2")) + if op2_files: + self.op2_path = max(op2_files, key=lambda p: p.stat().st_mtime) + break + + if self.op2_path is None: + return False + + try: + self.geo_path = self._find_geometry_file(self.op2_path) + return True + except FileNotFoundError: + return False + + def _find_geometry_file(self, op2_path: Path) -> Path: + """Find BDF/DAT geometry file for OP2.""" + folder = op2_path.parent + base = op2_path.stem + + for ext in ['.dat', '.bdf']: + cand = folder / (base + ext) + if cand.exists(): + return cand + + for f in folder.iterdir(): + if f.suffix.lower() in ['.dat', '.bdf']: + return f + + raise FileNotFoundError(f"No geometry file found for {op2_path}") + + def _load_data(self): + """Load geometry and displacement data.""" + if self._node_geo is not None: + return + + _load_dependencies() + + bdf = _BDF() + bdf.read_bdf(str(self.geo_path)) + self._node_geo = {int(nid): node.get_position() + for nid, node in bdf.nodes.items()} + + op2 = _OP2() + op2.read_op2(str(self.op2_path)) + + if not op2.displacements: + raise RuntimeError("No displacement data in OP2") + + self._displacements = {} + for key, darr in op2.displacements.items(): + data = darr.data + dmat = data[0] if data.ndim == 3 else (data if data.ndim == 2 else None) + if dmat is None: + continue + + ngt = darr.node_gridtype.astype(int) + node_ids = ngt if ngt.ndim == 1 else ngt[:, 0] + isubcase = getattr(darr, 'isubcase', None) + label = str(isubcase) if isubcase else str(key) + + self._displacements[label] = { + 'node_ids': node_ids.astype(int), + 'disp': dmat.copy() + } + + def _build_wfe_arrays( + self, + label: str, + disp_unit: str = 'mm' + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Build X, Y, WFE arrays for a subcase.""" + nm_per_unit = 1e6 if disp_unit == 'mm' else 1e9 + + data = self._displacements[label] + node_ids = data['node_ids'] + dmat = data['disp'] + + X, Y, WFE = [], [], [] + valid_nids = [] + for nid, vec in zip(node_ids, dmat): + geo = self._node_geo.get(int(nid)) + if geo is None: + continue + X.append(geo[0]) + Y.append(geo[1]) + wfe = vec[2] * 2.0 * nm_per_unit + WFE.append(wfe) + valid_nids.append(nid) + + return (np.array(X), np.array(Y), np.array(WFE), np.array(valid_nids)) + + def _compute_relative_wfe( + self, + X1, Y1, WFE1, nids1, + X2, Y2, WFE2, nids2 + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Compute WFE1 - WFE2 for common nodes.""" + ref_map = {int(nid): (x, y, w) for nid, x, y, w in zip(nids2, X2, Y2, WFE2)} + + X_rel, Y_rel, WFE_rel = [], [], [] + for nid, x, y, w in zip(nids1, X1, Y1, WFE1): + nid = int(nid) + if nid in ref_map: + _, _, w_ref = ref_map[nid] + X_rel.append(x) + Y_rel.append(y) + WFE_rel.append(w - w_ref) + + return np.array(X_rel), np.array(Y_rel), np.array(WFE_rel) + + def _analyze_subcase( + self, + X: np.ndarray, + Y: np.ndarray, + W_nm: np.ndarray, + cfg: Dict + ) -> Dict[str, Any]: + """ + Full MSF analysis for a single subcase. + + Returns: + Dict with coefficients, bands, RSS, dominant modes, mesh info + """ + n_modes = cfg['n_modes'] + lsf_max = cfg['lsf_max'] + msf_max = cfg['msf_max'] + + # Fit Zernike + coeffs, R = compute_zernike_coeffs(X, Y, W_nm, n_modes) + + # Band classification + bands = classify_modes_by_band(n_modes, lsf_max, msf_max) + + # RSS per band + rss = compute_band_rss(coeffs, bands) + + # Band percentages + total_rss = rss['total'] + pct = { + 'lsf': 100 * rss['lsf'] / total_rss if total_rss > 0 else 0, + 'msf': 100 * rss['msf'] / total_rss if total_rss > 0 else 0, + 'hsf': 100 * rss['hsf'] / total_rss if total_rss > 0 else 0, + } + + # Dominant MSF modes + dominant_msf = find_dominant_msf_modes(coeffs, bands, top_n=5) + + # MSF-only residual + W_msf = compute_band_residual(X, Y, W_nm, coeffs, R, bands, ['msf']) + + # Mesh resolution + mesh_info = estimate_mesh_resolution(X, Y) + + return { + 'coefficients': coeffs, + 'R': R, + 'bands': bands, + 'rss': rss, + 'rss_pct': pct, + 'dominant_msf': dominant_msf, + 'W_msf': W_msf, + 'mesh_info': mesh_info, + 'n_modes': n_modes, + 'lsf_max': lsf_max, + 'msf_max': msf_max, + } + + def _generate_html( + self, + analyses: Dict[str, Dict], + cfg: Dict, + X: np.ndarray, + Y: np.ndarray + ) -> str: + """Generate comprehensive HTML report.""" + _load_dependencies() + + n_modes = cfg['n_modes'] + amp = cfg.get('amp', 0.5) + pancake = cfg.get('pancake', 3.0) + downsample = cfg.get('plot_downsample', 15000) + colorscale = cfg.get('colorscale', 'Turbo') + + # Use first available analysis for visualization + primary_key = list(analyses.keys())[0] + primary = analyses[primary_key] + + coeffs = primary['coefficients'] + bands = primary['bands'] + W_msf = primary['W_msf'] + mesh_info = primary['mesh_info'] + + # Prepare band-colored bar chart data + bar_colors = [] + for j in range(1, n_modes + 1): + if j in bands['lsf']: + bar_colors.append(BAND_COLORS['lsf']) + elif j in bands['msf']: + bar_colors.append(BAND_COLORS['msf']) + else: + bar_colors.append(BAND_COLORS['hsf']) + + labels = [zernike_label(j) for j in range(1, n_modes + 1)] + coeff_abs = np.abs(coeffs) + + # Downsample for 3D plot + n = len(X) + if n > downsample: + rng = np.random.default_rng(42) + sel = rng.choice(n, size=downsample, replace=False) + Xp, Yp, Wp = X[sel], Y[sel], W_msf[sel] + else: + Xp, Yp, Wp = X, Y, W_msf + + res_amp = amp * Wp + max_amp = float(np.max(np.abs(res_amp))) if res_amp.size else 1.0 + + # Create figure with subplots + fig = _make_subplots( + rows=4, cols=1, + specs=[[{"type": "scene"}], [{"type": "table"}], + [{"type": "table"}], [{"type": "xy"}]], + row_heights=[0.40, 0.20, 0.15, 0.25], + vertical_spacing=0.03, + subplot_titles=[ + "MSF Content Only (n=11-50)", + "Band Decomposition (RSS)", + "Dominant MSF Modes", + "Coefficient Magnitude by Band" + ] + ) + + # 3D MSF surface + try: + tri = _Triangulation(Xp, Yp) + if tri.triangles is not None and len(tri.triangles) > 0: + i, j, k = tri.triangles.T + fig.add_trace(_go.Mesh3d( + x=Xp, y=Yp, z=res_amp, + i=i, j=j, k=k, + intensity=res_amp, + colorscale=colorscale, + opacity=1.0, + flatshading=False, + lighting=dict(ambient=0.4, diffuse=0.8, specular=0.3, + roughness=0.5, fresnel=0.2), + lightposition=dict(x=100, y=200, z=300), + showscale=True, + colorbar=dict(title=dict(text="MSF WFE (nm)", side="right"), + thickness=15, len=0.5, tickformat=".1f"), + hovertemplate="X: %{x:.1f}
Y: %{y:.1f}
MSF: %{z:.2f} nm" + ), row=1, col=1) + except Exception: + fig.add_trace(_go.Scatter3d( + x=Xp, y=Yp, z=res_amp, + mode='markers', + marker=dict(size=2, color=res_amp, colorscale=colorscale, showscale=True), + showlegend=False + ), row=1, col=1) + + # Configure 3D scene + fig.update_scenes( + camera=dict(eye=dict(x=1.2, y=1.2, z=0.8), up=dict(x=0, y=0, z=1)), + xaxis=dict(title="X (mm)", showgrid=True, + gridcolor='rgba(128,128,128,0.3)', + showbackground=True, backgroundcolor='rgba(240,240,240,0.9)'), + yaxis=dict(title="Y (mm)", showgrid=True, + gridcolor='rgba(128,128,128,0.3)', + showbackground=True, backgroundcolor='rgba(240,240,240,0.9)'), + zaxis=dict(title="MSF WFE (nm)", + range=[-max_amp * pancake, max_amp * pancake], + showgrid=True, gridcolor='rgba(128,128,128,0.3)', + showbackground=True, backgroundcolor='rgba(230,230,250,0.9)'), + aspectmode='manual', + aspectratio=dict(x=1, y=1, z=0.4) + ) + + # Band decomposition table + band_headers = ["Band", "Order Range", "Feature Size"] + band_values = [ + ["LSF (Low)", "MSF (Mid)", "HSF (High)", "Total", "Filtered (J5+)"], + [f"n ≤ {cfg['lsf_max']}", f"n = {cfg['lsf_max']+1}-{cfg['msf_max']}", + f"n > {cfg['msf_max']}", "All", "J5+ (excl. piston/tip/tilt/defocus)"], + [f"> {int(mesh_info['diameter_mm']/cfg['lsf_max'])} mm", + f"{int(mesh_info['diameter_mm']/cfg['msf_max'])}-{int(mesh_info['diameter_mm']/(cfg['lsf_max']+1))} mm", + f"< {int(mesh_info['diameter_mm']/cfg['msf_max'])} mm", "-", "M2 correctable removed"], + ] + + # Add columns for each analysis + for key, analysis in analyses.items(): + band_headers.append(f"{key} RSS (nm)") + rss = analysis['rss'] + pct = analysis['rss_pct'] + # Calculate filtered percentage + filtered_pct = 100 * rss['filtered'] / rss['total'] if rss['total'] > 0 else 0 + band_values.append([ + f"{rss['lsf']:.2f} ({pct['lsf']:.1f}%)", + f"{rss['msf']:.2f} ({pct['msf']:.1f}%)", + f"{rss['hsf']:.2f} ({pct['hsf']:.1f}%)", + f"{rss['total']:.2f}", + f"{rss['filtered']:.2f} ({filtered_pct:.1f}%)", + ]) + + fig.add_trace(_go.Table( + header=dict(values=band_headers, + align="left", fill_color='#1f2937', font=dict(color='white')), + cells=dict(values=band_values, + align="left", fill_color='#374151', font=dict(color='white')) + ), row=2, col=1) + + # Dominant MSF modes table + dom_modes = primary['dominant_msf'] + if dom_modes: + fig.add_trace(_go.Table( + header=dict(values=["Rank", "Mode", "|Coeff| (nm)", "Order n"], + align="left", fill_color='#1f2937', font=dict(color='white')), + cells=dict(values=[ + [f"#{i+1}" for i in range(len(dom_modes))], + [m['label'] for m in dom_modes], + [f"{m['coeff_nm']:.3f}" for m in dom_modes], + [str(m['n']) for m in dom_modes], + ], align="left", fill_color='#374151', font=dict(color='white')) + ), row=3, col=1) + + # Bar chart with band colors + fig.add_trace( + _go.Bar( + x=coeff_abs.tolist(), y=labels, + orientation='h', + marker_color=bar_colors, + hovertemplate="%{y}
|Coeff| = %{x:.3f} nm", + showlegend=False + ), + row=4, col=1 + ) + + # Add legend for bands + for band_name, color in BAND_COLORS.items(): + fig.add_trace(_go.Scatter( + x=[None], y=[None], + mode='markers', + marker=dict(size=10, color=color), + name=band_name.upper(), + showlegend=True + )) + + # Layout + fig.update_layout( + width=1400, height=1600, + margin=dict(t=60, b=20, l=20, r=20), + paper_bgcolor='#111827', plot_bgcolor='#1f2937', + font=dict(color='white'), + title=dict( + text=f"Atomizer MSF Analysis | {n_modes} modes | Mesh: {mesh_info['node_count']} nodes, {mesh_info['avg_spacing_mm']:.1f}mm spacing", + x=0.5, font=dict(size=16) + ), + legend=dict( + orientation="h", + yanchor="bottom", + y=1.02, + xanchor="right", + x=1 + ) + ) + + return fig.to_html(include_plotlyjs='cdn', full_html=True) + + def _generate(self, config: InsightConfig) -> InsightResult: + """Generate MSF analysis.""" + self._load_data() + + cfg = {**DEFAULT_CONFIG, **config.extra} + disp_unit = cfg['disp_unit'] + + # Map subcases + disps = self._displacements + if '1' in disps and '2' in disps: + sc_map = {'90°': '1', '20°': '2', '40°': '3', '60°': '4'} + elif '90' in disps and '20' in disps: + sc_map = {'90°': '90', '20°': '20', '40°': '40', '60°': '60'} + else: + available = sorted(disps.keys(), key=lambda x: int(x) if x.isdigit() else 0) + if len(available) >= 4: + sc_map = {'90°': available[0], '20°': available[1], + '40°': available[2], '60°': available[3]} + else: + return InsightResult(success=False, + error=f"Need 4 subcases, found: {available}") + + # Validate subcases + for angle, label in sc_map.items(): + if label not in disps: + return InsightResult(success=False, + error=f"Subcase '{label}' (angle {angle}) not found") + + # Load reference data + X_20, Y_20, WFE_20, nids_20 = self._build_wfe_arrays(sc_map['20°'], disp_unit) + X_90, Y_90, WFE_90, nids_90 = self._build_wfe_arrays(sc_map['90°'], disp_unit) + + analyses = {} + + # Analyze each key orientation + for angle in ['40°', '60°']: + label = sc_map[angle] + X, Y, WFE, nids = self._build_wfe_arrays(label, disp_unit) + + # Absolute + analyses[f"{angle} Abs"] = self._analyze_subcase(X, Y, WFE, cfg) + + # Relative to 20° + X_rel, Y_rel, WFE_rel = self._compute_relative_wfe( + X, Y, WFE, nids, X_20, Y_20, WFE_20, nids_20) + analyses[f"{angle} vs 20°"] = self._analyze_subcase(X_rel, Y_rel, WFE_rel, cfg) + + # Relative to 90° (manufacturing) + X_rel90, Y_rel90, WFE_rel90 = self._compute_relative_wfe( + X, Y, WFE, nids, X_90, Y_90, WFE_90, nids_90) + analyses[f"{angle} vs 90°"] = self._analyze_subcase(X_rel90, Y_rel90, WFE_rel90, cfg) + + # Also analyze 90° absolute (manufacturing baseline) + analyses["90° Mfg"] = self._analyze_subcase(X_90, Y_90, WFE_90, cfg) + + # Generate HTML using 40° vs 20° as primary visualization + primary_analysis = analyses["40° vs 20°"] + X_40, Y_40, _, _ = self._build_wfe_arrays(sc_map['40°'], disp_unit) + X_rel40, Y_rel40, _ = self._compute_relative_wfe( + X_40, Y_40, analyses["40° Abs"]['W_msf'], np.arange(len(X_40)), + X_20, Y_20, WFE_20, nids_20) + + html = self._generate_html(analyses, cfg, X_40, Y_40) + + timestamp = datetime.now().strftime("%Y%m%d_%H%M%S") + output_dir = config.output_dir or self.insights_path + output_dir.mkdir(parents=True, exist_ok=True) + + html_path = output_dir / f"msf_zernike_{timestamp}.html" + html_path.write_text(html, encoding='utf-8') + + # Build summary + summary = { + 'n_modes': cfg['n_modes'], + 'lsf_max_order': cfg['lsf_max'], + 'msf_max_order': cfg['msf_max'], + 'mesh_nodes': primary_analysis['mesh_info']['node_count'], + 'mesh_spacing_mm': primary_analysis['mesh_info']['avg_spacing_mm'], + 'max_resolvable_order': primary_analysis['mesh_info']['max_resolvable_order'], + } + + # Add key metrics for each analysis + for key, analysis in analyses.items(): + safe_key = key.replace('°', 'deg').replace(' ', '_') + summary[f'{safe_key}_lsf_rss'] = analysis['rss']['lsf'] + summary[f'{safe_key}_msf_rss'] = analysis['rss']['msf'] + summary[f'{safe_key}_hsf_rss'] = analysis['rss']['hsf'] + summary[f'{safe_key}_total_rss'] = analysis['rss']['total'] + summary[f'{safe_key}_filtered_rss'] = analysis['rss']['filtered'] + summary[f'{safe_key}_msf_pct'] = analysis['rss_pct']['msf'] + + return InsightResult( + success=True, + html_path=html_path, + summary=summary + ) diff --git a/optimization_engine/insights/zernike_dashboard.py b/optimization_engine/insights/zernike_dashboard.py new file mode 100644 index 00000000..55e00b88 --- /dev/null +++ b/optimization_engine/insights/zernike_dashboard.py @@ -0,0 +1,1304 @@ +""" +Zernike Dashboard - Comprehensive Mirror Wavefront Analysis + +A complete, professional dashboard for M1 mirror optimization that combines: +- Executive summary with all key metrics at a glance +- All orientation views (40°, 60°, 90°) on a single page +- MSF band analysis integrated +- Light/white theme for better readability +- Interactive comparison between subcases + +This is the unified replacement for the separate zernike_wfe and msf_zernike insights. + +Author: Atomizer Framework +Date: 2024-12-22 +""" + +from pathlib import Path +from datetime import datetime +from typing import Dict, Any, List, Optional, Tuple +import numpy as np +from math import factorial +from numpy.linalg import LinAlgError + +from .base import StudyInsight, InsightConfig, InsightResult, register_insight + +# Lazy imports +_plotly_loaded = False +_go = None +_make_subplots = None +_Triangulation = None +_OP2 = None +_BDF = None +_ZernikeOPDExtractor = None + + +def _load_dependencies(): + """Lazy load heavy dependencies.""" + global _plotly_loaded, _go, _make_subplots, _Triangulation, _OP2, _BDF, _ZernikeOPDExtractor + if not _plotly_loaded: + import plotly.graph_objects as go + from plotly.subplots import make_subplots + from matplotlib.tri import Triangulation + from pyNastran.op2.op2 import OP2 + from pyNastran.bdf.bdf import BDF + from ..extractors.extract_zernike_figure import ZernikeOPDExtractor + _go = go + _make_subplots = make_subplots + _Triangulation = Triangulation + _OP2 = OP2 + _BDF = BDF + _ZernikeOPDExtractor = ZernikeOPDExtractor + _plotly_loaded = True + + +# ============================================================================ +# Color Themes +# ============================================================================ + +LIGHT_THEME = { + 'bg': '#ffffff', + 'card_bg': '#f8fafc', + 'card_border': '#e2e8f0', + 'text': '#1e293b', + 'text_muted': '#64748b', + 'accent': '#3b82f6', + 'success': '#10b981', + 'warning': '#f59e0b', + 'danger': '#ef4444', + 'grid': 'rgba(100,116,139,0.2)', + 'surface_bg': '#f1f5f9', +} + +BAND_COLORS = { + 'lsf': '#3b82f6', # Blue + 'msf': '#f59e0b', # Amber + 'hsf': '#ef4444', # Red +} + +ANGLE_COLORS = { + '40': '#10b981', # Emerald + '60': '#3b82f6', # Blue + '90': '#8b5cf6', # Purple +} + + +# ============================================================================ +# Zernike Mathematics +# ============================================================================ + +def noll_indices(j: int) -> Tuple[int, int]: + """Convert Noll index to (n, m) radial/azimuthal orders.""" + if j < 1: + raise ValueError("Noll index j must be >= 1") + count = 0 + n = 0 + while True: + if n == 0: + ms = [0] + elif n % 2 == 0: + ms = [0] + [m for k in range(1, n//2 + 1) for m in (-2*k, 2*k)] + else: + ms = [m for k in range(0, (n+1)//2) for m in (-(2*k+1), (2*k+1))] + for m in ms: + count += 1 + if count == j: + return n, m + n += 1 + + +def noll_to_radial_order(j: int) -> int: + """Get radial order n for Noll index j.""" + n, _ = noll_indices(j) + return n + + +def zernike_noll(j: int, r: np.ndarray, th: np.ndarray) -> np.ndarray: + """Evaluate Zernike polynomial j at (r, theta).""" + n, m = noll_indices(j) + R = np.zeros_like(r) + for s in range((n - abs(m)) // 2 + 1): + c = ((-1)**s * factorial(n - s) / + (factorial(s) * + factorial((n + abs(m)) // 2 - s) * + factorial((n - abs(m)) // 2 - s))) + R += c * r**(n - 2*s) + if m == 0: + return R + return R * (np.cos(m * th) if m > 0 else np.sin(-m * th)) + + +def zernike_common_name(n: int, m: int) -> str: + """Get common name for Zernike mode.""" + names = { + (0, 0): "Piston", (1, -1): "Tilt X", (1, 1): "Tilt Y", + (2, 0): "Defocus", (2, -2): "Astig 45°", (2, 2): "Astig 0°", + (3, -1): "Coma X", (3, 1): "Coma Y", (3, -3): "Trefoil X", (3, 3): "Trefoil Y", + (4, 0): "Spherical", (4, -2): "Sec Astig X", (4, 2): "Sec Astig Y", + (4, -4): "Quadrafoil X", (4, 4): "Quadrafoil Y", + (5, -1): "Sec Coma X", (5, 1): "Sec Coma Y", + (6, 0): "Sec Spherical", + } + return names.get((n, m), f"Z({n},{m})") + + +def compute_zernike_coeffs( + X: np.ndarray, + Y: np.ndarray, + vals: np.ndarray, + n_modes: int, + chunk_size: int = 100000 +) -> Tuple[np.ndarray, float]: + """Fit Zernike coefficients to WFE data.""" + Xc, Yc = X - np.mean(X), Y - np.mean(Y) + R = float(np.max(np.hypot(Xc, Yc))) + r = np.hypot(Xc / R, Yc / R).astype(np.float32) + th = np.arctan2(Yc, Xc).astype(np.float32) + mask = (r <= 1.0) & ~np.isnan(vals) + if not np.any(mask): + raise RuntimeError("No valid points inside unit disk.") + idx = np.nonzero(mask)[0] + m = int(n_modes) + G = np.zeros((m, m), dtype=np.float64) + h = np.zeros((m,), dtype=np.float64) + v = vals.astype(np.float64) + + for start in range(0, len(idx), chunk_size): + sl = idx[start:start + chunk_size] + r_b, th_b, v_b = r[sl], th[sl], v[sl] + Zb = np.column_stack([zernike_noll(j, r_b, th_b).astype(np.float32) + for j in range(1, m + 1)]) + G += (Zb.T @ Zb).astype(np.float64) + h += (Zb.T @ v_b).astype(np.float64) + try: + coeffs = np.linalg.solve(G, h) + except LinAlgError: + coeffs = np.linalg.lstsq(G, h, rcond=None)[0] + return coeffs, R + + +def classify_modes_by_band( + n_modes: int, + lsf_max: int = 10, + msf_max: int = 50 +) -> Dict[str, List[int]]: + """Classify Zernike modes into LSF/MSF/HSF bands.""" + bands = {'lsf': [], 'msf': [], 'hsf': []} + for j in range(1, n_modes + 1): + n = noll_to_radial_order(j) + if n <= lsf_max: + bands['lsf'].append(j) + elif n <= msf_max: + bands['msf'].append(j) + else: + bands['hsf'].append(j) + return bands + + +def compute_band_rss(coeffs: np.ndarray, bands: Dict[str, List[int]]) -> Dict[str, float]: + """Compute RSS (Root Sum Square) of coefficients in each band.""" + rss = {} + for band_name, indices in bands.items(): + band_coeffs = [coeffs[j-1] for j in indices if j-1 < len(coeffs)] + rss[band_name] = float(np.sqrt(np.sum(np.array(band_coeffs)**2))) + rss['total'] = float(np.sqrt(np.sum(coeffs**2))) + # Filtered: exclude J1-J4 (piston, tip, tilt, defocus) + filtered_coeffs = coeffs[4:] if len(coeffs) > 4 else np.array([]) + rss['filtered'] = float(np.sqrt(np.sum(filtered_coeffs**2))) + return rss + + +# ============================================================================ +# Default Configuration +# ============================================================================ + +DEFAULT_CONFIG = { + 'n_modes': 50, + 'amp': 0.5, + 'pancake': 3.0, + 'plot_downsample': 8000, + 'filter_low_orders': 4, + 'lsf_max': 10, + 'msf_max': 50, + 'disp_unit': 'mm', + 'use_opd': True, # Use OPD method by default +} + + +@register_insight +class ZernikeDashboardInsight(StudyInsight): + """ + Comprehensive Zernike Dashboard for M1 Mirror Optimization. + + Generates a single-page dashboard with: + - Executive summary with key metrics + - All orientations (40°, 60°, 90°) in comparison view + - MSF band analysis + - Interactive surface plots + - Professional light theme + """ + + insight_type = "zernike_dashboard" + name = "Zernike Dashboard" + description = "Comprehensive mirror WFE dashboard with all orientations and MSF analysis" + category = "optical" + applicable_to = ["mirror", "optics", "wfe"] + required_files = ["*.op2"] + + def __init__(self, study_path: Path): + super().__init__(study_path) + self.op2_path: Optional[Path] = None + self.geo_path: Optional[Path] = None + self._node_geo: Optional[Dict] = None + self._displacements: Optional[Dict] = None + + def can_generate(self) -> bool: + """Check if OP2 and geometry files exist.""" + search_paths = [ + self.results_path, + self.study_path / "2_iterations", + self.setup_path / "model", + ] + + for search_path in search_paths: + if not search_path.exists(): + continue + op2_files = list(search_path.glob("**/*solution*.op2")) + if not op2_files: + op2_files = list(search_path.glob("**/*.op2")) + if op2_files: + self.op2_path = max(op2_files, key=lambda p: p.stat().st_mtime) + break + + if self.op2_path is None: + return False + + try: + self.geo_path = self._find_geometry_file(self.op2_path) + return True + except FileNotFoundError: + return False + + def _find_geometry_file(self, op2_path: Path) -> Path: + """Find BDF/DAT geometry file for OP2.""" + folder = op2_path.parent + base = op2_path.stem + + for ext in ['.dat', '.bdf']: + cand = folder / (base + ext) + if cand.exists(): + return cand + + for f in folder.iterdir(): + if f.suffix.lower() in ['.dat', '.bdf']: + return f + + raise FileNotFoundError(f"No geometry file found for {op2_path}") + + def _load_data(self): + """Load geometry and displacement data.""" + if self._node_geo is not None: + return + + _load_dependencies() + + bdf = _BDF() + bdf.read_bdf(str(self.geo_path)) + self._node_geo = {int(nid): node.get_position() + for nid, node in bdf.nodes.items()} + + op2 = _OP2() + op2.read_op2(str(self.op2_path)) + + if not op2.displacements: + raise RuntimeError("No displacement data in OP2") + + self._displacements = {} + for key, darr in op2.displacements.items(): + data = darr.data + dmat = data[0] if data.ndim == 3 else (data if data.ndim == 2 else None) + if dmat is None: + continue + + ngt = darr.node_gridtype.astype(int) + node_ids = ngt if ngt.ndim == 1 else ngt[:, 0] + isubcase = getattr(darr, 'isubcase', None) + label = str(isubcase) if isubcase else str(key) + + self._displacements[label] = { + 'node_ids': node_ids.astype(int), + 'disp': dmat.copy() + } + + def _build_wfe_arrays_standard( + self, + label: str, + disp_unit: str = 'mm' + ) -> Dict[str, np.ndarray]: + """Build arrays using standard Z-only method.""" + nm_per_unit = 1e6 if disp_unit == 'mm' else 1e9 + + data = self._displacements[label] + node_ids = data['node_ids'] + dmat = data['disp'] + + X, Y, WFE = [], [], [] + valid_nids = [] + dx_arr, dy_arr, dz_arr = [], [], [] + + for nid, vec in zip(node_ids, dmat): + geo = self._node_geo.get(int(nid)) + if geo is None: + continue + X.append(geo[0]) + Y.append(geo[1]) + wfe = vec[2] * 2.0 * nm_per_unit + WFE.append(wfe) + valid_nids.append(nid) + dx_arr.append(vec[0]) + dy_arr.append(vec[1]) + dz_arr.append(vec[2]) + + return { + 'X': np.array(X), + 'Y': np.array(Y), + 'WFE': np.array(WFE), + 'node_ids': np.array(valid_nids), + 'dx': np.array(dx_arr), + 'dy': np.array(dy_arr), + 'dz': np.array(dz_arr), + 'lateral_disp': np.sqrt(np.array(dx_arr)**2 + np.array(dy_arr)**2), + } + + def _build_wfe_arrays_opd( + self, + label: str, + disp_unit: str = 'mm' + ) -> Dict[str, np.ndarray]: + """Build arrays using OPD method (recommended).""" + _load_dependencies() + + extractor = _ZernikeOPDExtractor( + self.op2_path, + bdf_path=self.geo_path, + displacement_unit=disp_unit + ) + + opd_data = extractor._build_figure_opd_data(label) + + return { + 'X': opd_data['x_deformed'], + 'Y': opd_data['y_deformed'], + 'WFE': opd_data['wfe_nm'], + 'node_ids': opd_data['node_ids'], + 'dx': opd_data['dx'], + 'dy': opd_data['dy'], + 'dz': opd_data['dz'], + 'lateral_disp': opd_data['lateral_disp'], + 'x_original': opd_data['x_original'], + 'y_original': opd_data['y_original'], + } + + def _build_wfe_arrays( + self, + label: str, + disp_unit: str = 'mm', + use_opd: bool = True + ) -> Dict[str, np.ndarray]: + """Build X, Y, WFE arrays for a subcase.""" + if use_opd: + return self._build_wfe_arrays_opd(label, disp_unit) + else: + return self._build_wfe_arrays_standard(label, disp_unit) + + def _compute_relative_wfe( + self, + data1: Dict[str, np.ndarray], + data2: Dict[str, np.ndarray] + ) -> Dict[str, np.ndarray]: + """Compute relative WFE (target - reference) for common nodes.""" + X1, Y1, WFE1, nids1 = data1['X'], data1['Y'], data1['WFE'], data1['node_ids'] + X2, Y2, WFE2, nids2 = data2['X'], data2['Y'], data2['WFE'], data2['node_ids'] + + ref_map = {int(nid): (x, y, w) for nid, x, y, w in zip(nids2, X2, Y2, WFE2)} + + dx1_map = {int(nid): d for nid, d in zip(data1['node_ids'], data1['dx'])} + dy1_map = {int(nid): d for nid, d in zip(data1['node_ids'], data1['dy'])} + dz1_map = {int(nid): d for nid, d in zip(data1['node_ids'], data1['dz'])} + dx2_map = {int(nid): d for nid, d in zip(data2['node_ids'], data2['dx'])} + dy2_map = {int(nid): d for nid, d in zip(data2['node_ids'], data2['dy'])} + dz2_map = {int(nid): d for nid, d in zip(data2['node_ids'], data2['dz'])} + + X_rel, Y_rel, WFE_rel, nids_rel = [], [], [], [] + dx_rel, dy_rel, dz_rel = [], [], [] + + for nid, x, y, w in zip(nids1, X1, Y1, WFE1): + nid = int(nid) + if nid not in ref_map or nid not in dx2_map: + continue + + _, _, w_ref = ref_map[nid] + X_rel.append(x) + Y_rel.append(y) + WFE_rel.append(w - w_ref) + nids_rel.append(nid) + dx_rel.append(dx1_map[nid] - dx2_map[nid]) + dy_rel.append(dy1_map[nid] - dy2_map[nid]) + dz_rel.append(dz1_map[nid] - dz2_map[nid]) + + return { + 'X': np.array(X_rel), + 'Y': np.array(Y_rel), + 'WFE': np.array(WFE_rel), + 'node_ids': np.array(nids_rel), + 'dx': np.array(dx_rel), + 'dy': np.array(dy_rel), + 'dz': np.array(dz_rel), + 'lateral_disp': np.sqrt(np.array(dx_rel)**2 + np.array(dy_rel)**2) if dx_rel else np.array([]), + } + + def _compute_metrics( + self, + X: np.ndarray, + Y: np.ndarray, + W_nm: np.ndarray, + n_modes: int, + filter_orders: int, + cfg: Dict + ) -> Dict[str, Any]: + """Compute comprehensive metrics including band analysis.""" + coeffs, R = compute_zernike_coeffs(X, Y, W_nm, n_modes) + + Xc = X - np.mean(X) + Yc = Y - np.mean(Y) + r = np.hypot(Xc / R, Yc / R) + th = np.arctan2(Yc, Xc) + + Z = np.column_stack([zernike_noll(j, r, th) for j in range(1, n_modes + 1)]) + W_res_filt = W_nm - Z[:, :filter_orders].dot(coeffs[:filter_orders]) + W_res_filt_j1to3 = W_nm - Z[:, :3].dot(coeffs[:3]) + + # Band analysis + bands = classify_modes_by_band(n_modes, cfg['lsf_max'], cfg['msf_max']) + band_rss = compute_band_rss(coeffs, bands) + + # Aberration magnitudes + aberrations = { + 'defocus': float(abs(coeffs[3])) if len(coeffs) > 3 else 0.0, + 'astigmatism': float(np.sqrt(coeffs[4]**2 + coeffs[5]**2)) if len(coeffs) > 5 else 0.0, + 'coma': float(np.sqrt(coeffs[6]**2 + coeffs[7]**2)) if len(coeffs) > 7 else 0.0, + 'trefoil': float(np.sqrt(coeffs[8]**2 + coeffs[9]**2)) if len(coeffs) > 9 else 0.0, + 'spherical': float(abs(coeffs[10])) if len(coeffs) > 10 else 0.0, + } + + return { + 'coefficients': coeffs, + 'R': R, + 'global_rms': float(np.sqrt(np.mean(W_nm**2))), + 'filtered_rms': float(np.sqrt(np.mean(W_res_filt**2))), + 'rms_j1to3': float(np.sqrt(np.mean(W_res_filt_j1to3**2))), + 'W_res_filt': W_res_filt, + 'band_rss': band_rss, + 'aberrations': aberrations, + } + + def _generate_html( + self, + all_data: Dict[str, Dict], + cfg: Dict, + timestamp: str + ) -> str: + """Generate comprehensive HTML dashboard.""" + _load_dependencies() + + theme = LIGHT_THEME + + # Extract key metrics for executive summary + metrics_40 = all_data.get('40_vs_20', {}).get('metrics', {}) + metrics_60 = all_data.get('60_vs_20', {}).get('metrics', {}) + metrics_90 = all_data.get('90_abs', {}).get('metrics', {}) + + rms_40 = metrics_40.get('filtered_rms', 0) + rms_60 = metrics_60.get('filtered_rms', 0) + rms_90_mfg = metrics_90.get('rms_j1to3', 0) + + # Generate HTML + html = f''' + + + + + Zernike Dashboard - Atomizer + + + + +
+ +
+
+

Zernike WFE Dashboard

+
M1 Mirror Optimization Analysis • OPD Method (X,Y,Z corrected)
+
+
Generated: {timestamp}
+
+ + +
+
+
40° vs 20° RMS
+
+ {rms_40:.2f} nm +
+
Target: ≤ 4.0 nm
+
+ +
+
60° vs 20° RMS
+
+ {rms_60:.2f} nm +
+
Target: ≤ 10.0 nm
+
+ +
+
90° MFG (J1-J3)
+
+ {rms_90_mfg:.2f} nm +
+
Target: ≤ 20.0 nm
+
+ +
+
Weighted Score
+
+ {self._compute_weighted_score(rms_40, rms_60, rms_90_mfg):.1f} +
+
Lower is better
+
+
+ + +
+
+
+
Detailed Metrics Comparison
+
All values in nm • J1-J4 filtered (piston, tip, tilt, defocus removed)
+
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Metric40° vs 20°60° vs 20°90° MFGTarget
Filtered RMS (J1-J4){rms_40:.2f}{rms_60:.2f}—≤ 4 / 10 nm
Optician Workload (J1-J3)——{rms_90_mfg:.2f}≤ 20 nm
Astigmatism (J5+J6){metrics_40.get('aberrations', {}).get('astigmatism', 0):.2f}{metrics_60.get('aberrations', {}).get('astigmatism', 0):.2f}{metrics_90.get('aberrations', {}).get('astigmatism', 0):.2f}—
Coma (J7+J8){metrics_40.get('aberrations', {}).get('coma', 0):.2f}{metrics_60.get('aberrations', {}).get('coma', 0):.2f}{metrics_90.get('aberrations', {}).get('coma', 0):.2f}—
Trefoil (J9+J10){metrics_40.get('aberrations', {}).get('trefoil', 0):.2f}{metrics_60.get('aberrations', {}).get('trefoil', 0):.2f}{metrics_90.get('aberrations', {}).get('trefoil', 0):.2f}—
Spherical (J11){metrics_40.get('aberrations', {}).get('spherical', 0):.2f}{metrics_60.get('aberrations', {}).get('spherical', 0):.2f}{metrics_90.get('aberrations', {}).get('spherical', 0):.2f}—
+
+
+ + +
+
+
+
Spatial Frequency Band Analysis
+
RSS decomposition • LSF: n≤10 (M2 correctable) • MSF: n=11-50 (support print-through) • HSF: n>50
+
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Band40° vs 20°60° vs 20°90° MFGPhysics
● LSF (n ≤ 10){metrics_40.get('band_rss', {}).get('lsf', 0):.2f} nm{metrics_60.get('band_rss', {}).get('lsf', 0):.2f} nm{metrics_90.get('band_rss', {}).get('lsf', 0):.2f} nmM2 hexapod correctable
● MSF (n = 11-50){metrics_40.get('band_rss', {}).get('msf', 0):.2f} nm{metrics_60.get('band_rss', {}).get('msf', 0):.2f} nm{metrics_90.get('band_rss', {}).get('msf', 0):.2f} nmSupport print-through
● HSF (n > 50){metrics_40.get('band_rss', {}).get('hsf', 0):.2f} nm{metrics_60.get('band_rss', {}).get('hsf', 0):.2f} nm{metrics_90.get('band_rss', {}).get('hsf', 0):.2f} nmHigh frequency (mesh limit)
Total RSS{metrics_40.get('band_rss', {}).get('total', 0):.2f} nm{metrics_60.get('band_rss', {}).get('total', 0):.2f} nm{metrics_90.get('band_rss', {}).get('total', 0):.2f} nm
+
+
+ + +
+
+
+
Surface Residual Visualization
+
WFE residual after J1-J4 removal • Visual amplification applied for clarity
+
+
+
+
+
+
40° vs 20° (Relative)
+
+
+
+
60° vs 20° (Relative)
+
+
+
+
90° Manufacturing (Absolute)
+
+
+
+
+
+ + +
+
+
+
Zernike Coefficient Magnitudes
+
Absolute values • Color-coded by spatial frequency band
+
+
+
+
+
+
+
+ LSF (n ≤ 10) +
+
+
+ MSF (n = 11-50) +
+
+
+ HSF (n > 50) +
+
+
+
+
+ + + +''' + + return html + + def _compute_weighted_score(self, rms_40: float, rms_60: float, rms_90: float) -> float: + """Compute weighted optimization score.""" + # Standard weights from M1 optimization + w40 = 100 / 4.0 # Target 4 nm + w60 = 50 / 10.0 # Target 10 nm + w90 = 20 / 20.0 # Target 20 nm + return w40 * rms_40 + w60 * rms_60 + w90 * rms_90 + + def _get_status_class(self, value: float, target: float, warn: float) -> str: + """Get CSS class for value status.""" + if value <= target: + return ' pass' + elif value <= warn: + return ' warn' + else: + return ' fail' + + def _generate_plot_js(self, all_data: Dict[str, Dict], cfg: Dict) -> str: + """Generate JavaScript for Plotly plots.""" + _load_dependencies() + + amp = cfg.get('amp', 0.5) + downsample = cfg.get('plot_downsample', 8000) + + js_parts = [] + + # Generate 3D surface plots for each angle + for key, angle, div_id in [ + ('40_vs_20', '40', 'plot-40'), + ('60_vs_20', '60', 'plot-60'), + ('90_abs', '90', 'plot-90'), + ]: + data = all_data.get(key, {}) + X = data.get('X', np.array([])) + Y = data.get('Y', np.array([])) + W = data.get('metrics', {}).get('W_res_filt', np.array([])) + + if len(X) == 0: + continue + + # Downsample + n = len(X) + if n > downsample: + rng = np.random.default_rng(42) + sel = rng.choice(n, size=downsample, replace=False) + Xp, Yp, Wp = X[sel], Y[sel], W[sel] + else: + Xp, Yp, Wp = X, Y, W + + res_amp = amp * Wp + max_amp = float(np.max(np.abs(res_amp))) if res_amp.size else 1.0 + + # Build triangulation + try: + tri = _Triangulation(Xp, Yp) + i_arr = tri.triangles[:, 0].tolist() + j_arr = tri.triangles[:, 1].tolist() + k_arr = tri.triangles[:, 2].tolist() + except: + i_arr, j_arr, k_arr = [], [], [] + + color = ANGLE_COLORS[angle] + + js_parts.append(f''' + // {angle}° Plot + Plotly.newPlot('{div_id}', [{{ + type: 'mesh3d', + x: {Xp.tolist()}, + y: {Yp.tolist()}, + z: {res_amp.tolist()}, + i: {i_arr}, + j: {j_arr}, + k: {k_arr}, + intensity: {res_amp.tolist()}, + colorscale: 'RdBu_r', + showscale: true, + colorbar: {{ + title: 'WFE (nm)', + thickness: 12, + len: 0.6 + }}, + lighting: {{ + ambient: 0.5, + diffuse: 0.8, + specular: 0.3 + }}, + hovertemplate: 'X: %{{x:.1f}} mm
Y: %{{y:.1f}} mm
WFE: %{{z:.2f}} nm' + }}], {{ + margin: {{t: 10, b: 10, l: 10, r: 10}}, + scene: {{ + camera: {{eye: {{x: 1.2, y: 1.2, z: 0.8}}}}, + xaxis: {{title: 'X (mm)', showgrid: true, gridcolor: 'rgba(100,100,100,0.2)'}}, + yaxis: {{title: 'Y (mm)', showgrid: true, gridcolor: 'rgba(100,100,100,0.2)'}}, + zaxis: {{title: 'WFE (nm)', range: [{-max_amp * 3}, {max_amp * 3}], showgrid: true, gridcolor: 'rgba(100,100,100,0.2)'}}, + aspectmode: 'manual', + aspectratio: {{x: 1, y: 1, z: 0.4}}, + bgcolor: '#f8fafc' + }}, + paper_bgcolor: '#ffffff' + }}, {{responsive: true}}); + ''') + + # Generate coefficient bar chart + metrics_40 = all_data.get('40_vs_20', {}).get('metrics', {}) + coeffs = metrics_40.get('coefficients', np.array([])) + n_modes = len(coeffs) if len(coeffs) > 0 else 50 + + if len(coeffs) > 0: + bands = classify_modes_by_band(n_modes, cfg.get('lsf_max', 10), cfg.get('msf_max', 50)) + bar_colors = [] + labels = [] + for j in range(1, n_modes + 1): + n, m = noll_indices(j) + labels.append(f'J{j} {zernike_common_name(n, m)}') + if j in bands['lsf']: + bar_colors.append(BAND_COLORS['lsf']) + elif j in bands['msf']: + bar_colors.append(BAND_COLORS['msf']) + else: + bar_colors.append(BAND_COLORS['hsf']) + + coeff_abs = np.abs(coeffs).tolist() + + js_parts.append(f''' + // Coefficient Bar Chart + Plotly.newPlot('coeff-chart', [{{ + type: 'bar', + x: {coeff_abs}, + y: {labels}, + orientation: 'h', + marker: {{ + color: {bar_colors} + }}, + hovertemplate: '%{{y}}
|Coeff| = %{{x:.3f}} nm' + }}], {{ + margin: {{t: 20, b: 40, l: 200, r: 20}}, + height: 800, + xaxis: {{title: '|Coefficient| (nm)', gridcolor: 'rgba(100,100,100,0.2)'}}, + yaxis: {{autorange: 'reversed'}}, + paper_bgcolor: '#ffffff', + plot_bgcolor: '#f8fafc' + }}, {{responsive: true}}); + ''') + + return '\n'.join(js_parts) + + def _generate(self, config: InsightConfig) -> InsightResult: + """Generate comprehensive Zernike dashboard.""" + self._load_data() + + cfg = {**DEFAULT_CONFIG, **config.extra} + n_modes = cfg['n_modes'] + filter_orders = cfg['filter_low_orders'] + disp_unit = cfg['disp_unit'] + use_opd = cfg['use_opd'] + + # Map subcases + disps = self._displacements + if '1' in disps and '2' in disps: + sc_map = {'90': '1', '20': '2', '40': '3', '60': '4'} + elif '90' in disps and '20' in disps: + sc_map = {'90': '90', '20': '20', '40': '40', '60': '60'} + else: + available = sorted(disps.keys(), key=lambda x: int(x) if x.isdigit() else 0) + if len(available) >= 4: + sc_map = {'90': available[0], '20': available[1], + '40': available[2], '60': available[3]} + else: + return InsightResult(success=False, + error=f"Need 4 subcases, found: {available}") + + # Validate + for angle, label in sc_map.items(): + if label not in disps: + return InsightResult(success=False, + error=f"Subcase '{label}' (angle {angle}) not found") + + timestamp = datetime.now().strftime("%Y-%m-%d %H:%M:%S") + output_dir = config.output_dir or self.insights_path + output_dir.mkdir(parents=True, exist_ok=True) + + all_data = {} + + # Load reference data (20°) + data_ref = self._build_wfe_arrays(sc_map['20'], disp_unit, use_opd) + + # Process each angle + for angle in ['40', '60']: + data = self._build_wfe_arrays(sc_map[angle], disp_unit, use_opd) + rel_data = self._compute_relative_wfe(data, data_ref) + + metrics = self._compute_metrics( + rel_data['X'], rel_data['Y'], rel_data['WFE'], + n_modes, filter_orders, cfg + ) + + all_data[f'{angle}_vs_20'] = { + 'X': rel_data['X'], + 'Y': rel_data['Y'], + 'WFE': rel_data['WFE'], + 'metrics': metrics, + } + + # 90° absolute (manufacturing) + data_90 = self._build_wfe_arrays(sc_map['90'], disp_unit, use_opd) + metrics_90 = self._compute_metrics( + data_90['X'], data_90['Y'], data_90['WFE'], + n_modes, filter_orders, cfg + ) + all_data['90_abs'] = { + 'X': data_90['X'], + 'Y': data_90['Y'], + 'WFE': data_90['WFE'], + 'metrics': metrics_90, + } + + # Generate HTML + html = self._generate_html(all_data, cfg, timestamp) + + file_timestamp = datetime.now().strftime("%Y%m%d_%H%M%S") + html_path = output_dir / f"zernike_dashboard_{file_timestamp}.html" + html_path.write_text(html, encoding='utf-8') + + # Build summary + summary = { + '40_vs_20_filtered_rms': all_data['40_vs_20']['metrics']['filtered_rms'], + '60_vs_20_filtered_rms': all_data['60_vs_20']['metrics']['filtered_rms'], + '90_mfg_rms_j1to3': all_data['90_abs']['metrics']['rms_j1to3'], + '90_mfg_filtered_rms': all_data['90_abs']['metrics']['filtered_rms'], + 'weighted_score': self._compute_weighted_score( + all_data['40_vs_20']['metrics']['filtered_rms'], + all_data['60_vs_20']['metrics']['filtered_rms'], + all_data['90_abs']['metrics']['rms_j1to3'] + ), + 'method': 'OPD' if use_opd else 'Standard', + 'n_modes': n_modes, + } + + return InsightResult( + success=True, + html_path=html_path, + summary=summary + ) diff --git a/optimization_engine/insights/zernike_opd_comparison.py b/optimization_engine/insights/zernike_opd_comparison.py new file mode 100644 index 00000000..49a6c788 --- /dev/null +++ b/optimization_engine/insights/zernike_opd_comparison.py @@ -0,0 +1,625 @@ +""" +Zernike OPD Method Comparison Insight +===================================== + +Compares standard Z-only Zernike method against rigorous OPD method that +accounts for lateral (X, Y) displacements. Critical for understanding +whether your optimization objective is capturing true optical performance. + +WHY THIS MATTERS: +----------------- +When supports pinch or laterally deform the mirror, nodes shift in X and Y. +The standard Zernike method ignores this, potentially leading to: +- Optimized designs that "cheat" by distorting laterally +- False low WFE when the actual optical surface is degraded +- Optimizer convergence to non-optimal solutions + +This insight visualizes: +1. The difference between methods (where they disagree) +2. Lateral displacement magnitude map (where pinching occurs) +3. Which method you should be using for your study + +Applicable to: All mirror/optics studies, CRITICAL for lateral support optimization. +""" + +from pathlib import Path +from datetime import datetime +from typing import Dict, Any, Optional, Tuple +import numpy as np + +from .base import StudyInsight, InsightConfig, InsightResult, register_insight + +# Lazy imports +_plotly_loaded = False +_go = None +_make_subplots = None +_Triangulation = None + + +def _load_dependencies(): + """Lazy load heavy dependencies.""" + global _plotly_loaded, _go, _make_subplots, _Triangulation + if not _plotly_loaded: + import plotly.graph_objects as go + from plotly.subplots import make_subplots + from matplotlib.tri import Triangulation + _go = go + _make_subplots = make_subplots + _Triangulation = Triangulation + _plotly_loaded = True + + +# Default configuration +DEFAULT_CONFIG = { + 'amp': 0.5, # Visual deformation scale + 'pancake': 3.0, # Z-axis range multiplier + 'plot_downsample': 15000, + 'colorscale_wfe': 'RdBu', # Diverging for WFE difference + 'colorscale_lateral': 'Hot', # Sequential for lateral magnitude + 'disp_unit': 'mm', + 'n_modes': 50, + 'filter_orders': 4, +} + + +@register_insight +class ZernikeOPDComparisonInsight(StudyInsight): + """ + Compare standard vs OPD-rigorous Zernike methods. + + Generates visualizations showing: + - Lateral displacement magnitude map (where pinching occurs) + - WFE difference map (where methods disagree) + - Quantitative comparison tables + - Recommendation for which method to use + """ + + insight_type = "zernike_opd_comparison" + name = "Zernike Method Comparison (OPD vs Standard)" + description = "Compare Z-only vs rigorous OPD Zernike methods - reveals lateral displacement impact" + category = "optical" + applicable_to = ["mirror", "optics", "wfe", "lateral"] + required_files = ["*.op2"] + + def __init__(self, study_path: Path): + super().__init__(study_path) + self.op2_path: Optional[Path] = None + self.geo_path: Optional[Path] = None + self._node_geo: Optional[Dict] = None + self._displacements: Optional[Dict] = None + + def can_generate(self) -> bool: + """Check if OP2 and geometry files exist.""" + search_paths = [ + self.results_path, + self.study_path / "2_iterations", + self.setup_path / "model", + ] + + for search_path in search_paths: + if not search_path.exists(): + continue + op2_files = list(search_path.glob("**/*solution*.op2")) + if not op2_files: + op2_files = list(search_path.glob("**/*.op2")) + if op2_files: + self.op2_path = max(op2_files, key=lambda p: p.stat().st_mtime) + break + + if self.op2_path is None: + return False + + try: + self.geo_path = self._find_geometry_file(self.op2_path) + return True + except FileNotFoundError: + return False + + def _find_geometry_file(self, op2_path: Path) -> Path: + """Find BDF/DAT geometry file for OP2.""" + folder = op2_path.parent + base = op2_path.stem + + for ext in ['.dat', '.bdf']: + cand = folder / (base + ext) + if cand.exists(): + return cand + + for f in folder.iterdir(): + if f.suffix.lower() in ['.dat', '.bdf']: + return f + + raise FileNotFoundError(f"No geometry file found for {op2_path}") + + def _load_data(self): + """Load geometry and displacement data.""" + if self._node_geo is not None: + return + + _load_dependencies() + + from pyNastran.op2.op2 import OP2 + from pyNastran.bdf.bdf import BDF + + # Read geometry + bdf = BDF() + bdf.read_bdf(str(self.geo_path)) + self._node_geo = {int(nid): node.get_position() + for nid, node in bdf.nodes.items()} + + # Read displacements + op2 = OP2() + op2.read_op2(str(self.op2_path)) + + if not op2.displacements: + raise RuntimeError("No displacement data in OP2") + + self._displacements = {} + for key, darr in op2.displacements.items(): + data = darr.data + dmat = data[0] if data.ndim == 3 else (data if data.ndim == 2 else None) + if dmat is None: + continue + + ngt = darr.node_gridtype.astype(int) + node_ids = ngt if ngt.ndim == 1 else ngt[:, 0] + isubcase = getattr(darr, 'isubcase', None) + label = str(isubcase) if isubcase else str(key) + + self._displacements[label] = { + 'node_ids': node_ids.astype(int), + 'disp': dmat.copy() + } + + def _build_full_data( + self, + label: str, + disp_unit: str = 'mm' + ) -> Dict[str, np.ndarray]: + """Build complete data arrays including X, Y, Z displacements.""" + nm_per_unit = 1e6 if disp_unit == 'mm' else 1e9 + um_per_unit = nm_per_unit / 1000.0 + + data = self._displacements[label] + node_ids = data['node_ids'] + dmat = data['disp'] + + x0, y0, z0 = [], [], [] + dx, dy, dz = [], [], [] + valid_nids = [] + + for nid, vec in zip(node_ids, dmat): + geo = self._node_geo.get(int(nid)) + if geo is None: + continue + + x0.append(geo[0]) + y0.append(geo[1]) + z0.append(geo[2]) + dx.append(vec[0]) + dy.append(vec[1]) + dz.append(vec[2]) + valid_nids.append(nid) + + return { + 'x0': np.array(x0), + 'y0': np.array(y0), + 'z0': np.array(z0), + 'dx': np.array(dx), + 'dy': np.array(dy), + 'dz': np.array(dz), + 'nids': np.array(valid_nids), + 'nm_scale': nm_per_unit, + 'um_scale': um_per_unit, + } + + def _estimate_focal_length(self, x: np.ndarray, y: np.ndarray, z: np.ndarray) -> float: + """Estimate parabola focal length from geometry.""" + r_squared = x**2 + y**2 + mask = r_squared > 0 + if not np.any(mask): + return 5000.0 # Default fallback + + # Fit z = a * r² (assume centered parabola) + A = r_squared[mask].reshape(-1, 1) + b = z[mask] + + try: + a, _, _, _ = np.linalg.lstsq(A, b, rcond=None) + a = a[0] + except: + a = np.median(z[mask] / r_squared[mask]) + + if abs(a) < 1e-12: + return 5000.0 + + return 1.0 / (4.0 * abs(a)) + + def _compute_opd_comparison( + self, + data: Dict[str, np.ndarray], + n_modes: int, + filter_orders: int + ) -> Dict[str, Any]: + """Compute both methods and return comparison data.""" + from ..extractors.extract_zernike import compute_zernike_coefficients + + x0 = data['x0'] + y0 = data['y0'] + z0 = data['z0'] + dx = data['dx'] + dy = data['dy'] + dz = data['dz'] + nm_scale = data['nm_scale'] + um_scale = data['um_scale'] + + # Estimate focal length + focal = self._estimate_focal_length(x0, y0, z0) + + # === STANDARD METHOD (Z-only) === + # Uses original (x0, y0) with Z-displacement + wfe_standard = dz * 2.0 * nm_scale # WFE = 2 * surface error + + coeffs_std, R_std = compute_zernike_coefficients(x0, y0, wfe_standard, n_modes) + + # Compute filtered RMS + x_c = x0 - np.mean(x0) + y_c = y0 - np.mean(y0) + r = np.hypot(x_c / R_std, y_c / R_std) + th = np.arctan2(y_c, x_c) + + from ..extractors.extract_zernike import zernike_noll + Z_low = np.column_stack([zernike_noll(j, r, th) for j in range(1, filter_orders + 1)]) + wfe_std_filtered = wfe_standard - Z_low @ coeffs_std[:filter_orders] + rms_std_filtered = float(np.sqrt(np.mean(wfe_std_filtered**2))) + rms_std_global = float(np.sqrt(np.mean(wfe_standard**2))) + + # === RIGOROUS OPD METHOD === + # Uses deformed (x0+dx, y0+dy) and computes true surface error + # The CORRECT formulation uses DIFFERENTIAL approach: + # surface_error = dz - delta_z_parabola + # where delta_z_parabola is the Z change needed to stay on ideal parabola + + x_def = x0 + dx + y_def = y0 + dy + + # Compute the CHANGE in parabola Z due to lateral displacement + # For parabola: z = -r²/(4f), so: + # Δz_parabola = z(x_def, y_def) - z(x0, y0) + # = -[(x0+dx)² + (y0+dy)² - x0² - y0²] / (4f) + r0_sq = x0**2 + y0**2 + r_def_sq = x_def**2 + y_def**2 + delta_r_sq = r_def_sq - r0_sq + + # For concave mirror, z = -r²/(4f), so Δz = -Δr²/(4f) + delta_z_parabola = -delta_r_sq / (4.0 * focal) + + # TRUE surface error = actual dz - expected dz (to stay on parabola) + surface_error = dz - delta_z_parabola + wfe_opd = surface_error * 2.0 * nm_scale + + coeffs_opd, R_opd = compute_zernike_coefficients(x_def, y_def, wfe_opd, n_modes) + + x_c_def = x_def - np.mean(x_def) + y_c_def = y_def - np.mean(y_def) + r_def = np.hypot(x_c_def / R_opd, y_c_def / R_opd) + th_def = np.arctan2(y_c_def, x_c_def) + + Z_low_opd = np.column_stack([zernike_noll(j, r_def, th_def) for j in range(1, filter_orders + 1)]) + wfe_opd_filtered = wfe_opd - Z_low_opd @ coeffs_opd[:filter_orders] + rms_opd_filtered = float(np.sqrt(np.mean(wfe_opd_filtered**2))) + rms_opd_global = float(np.sqrt(np.mean(wfe_opd**2))) + + # === LATERAL DISPLACEMENT === + lateral_disp = np.sqrt(dx**2 + dy**2) * um_scale # in µm + + # === DIFFERENCE ANALYSIS === + # Where methods disagree (use original coords for visualization) + wfe_difference = wfe_opd - wfe_standard # Could be different due to coordinate system + # More meaningful: compare the residuals + + return { + # Coordinates for plotting + 'x_plot': x0, + 'y_plot': y0, + + # Standard method + 'wfe_standard': wfe_standard, + 'wfe_std_filtered': wfe_std_filtered, + 'rms_std_global': rms_std_global, + 'rms_std_filtered': rms_std_filtered, + + # OPD method + 'wfe_opd': wfe_opd, + 'wfe_opd_filtered': wfe_opd_filtered, + 'rms_opd_global': rms_opd_global, + 'rms_opd_filtered': rms_opd_filtered, + 'x_deformed': x_def, + 'y_deformed': y_def, + + # Lateral displacement + 'lateral_disp_um': lateral_disp, + 'max_lateral_um': float(np.max(lateral_disp)), + 'rms_lateral_um': float(np.sqrt(np.mean(lateral_disp**2))), + 'p99_lateral_um': float(np.percentile(lateral_disp, 99)), + + # Deltas + 'delta_global': rms_opd_global - rms_std_global, + 'delta_filtered': rms_opd_filtered - rms_std_filtered, + 'pct_diff_filtered': 100.0 * (rms_opd_filtered - rms_std_filtered) / rms_std_filtered if rms_std_filtered > 0 else 0.0, + + # Parabola info + 'focal_length': focal, + } + + def _get_recommendation(self, comparison: Dict[str, Any]) -> Tuple[str, str]: + """Generate recommendation based on comparison results.""" + max_lateral = comparison['max_lateral_um'] + pct_diff = abs(comparison['pct_diff_filtered']) + + if max_lateral > 10.0: + severity = "CRITICAL" + msg = (f"Large lateral displacements detected (max {max_lateral:.1f} µm). " + f"Methods differ by {pct_diff:.1f}%. " + f"You MUST use OPD method for accurate optimization.") + elif max_lateral > 1.0: + severity = "RECOMMENDED" + msg = (f"Significant lateral displacements (max {max_lateral:.2f} µm). " + f"Methods differ by {pct_diff:.1f}%. " + f"OPD method recommended for better accuracy.") + elif max_lateral > 0.1: + severity = "OPTIONAL" + msg = (f"Small lateral displacements (max {max_lateral:.3f} µm). " + f"Methods differ by {pct_diff:.2f}%. " + f"OPD method provides minor improvement.") + else: + severity = "EQUIVALENT" + msg = (f"Negligible lateral displacements (max {max_lateral:.4f} µm). " + f"Both methods give essentially identical results.") + + return severity, msg + + def _generate_html( + self, + title: str, + comparison: Dict[str, Any], + config: Dict[str, Any] + ) -> str: + """Generate comparison HTML with multiple views.""" + _load_dependencies() + + downsample = config.get('plot_downsample', 15000) + colorscale_wfe = config.get('colorscale_wfe', 'RdBu') + colorscale_lateral = config.get('colorscale_lateral', 'Hot') + amp = config.get('amp', 0.5) + + x = comparison['x_plot'] + y = comparison['y_plot'] + lateral = comparison['lateral_disp_um'] + wfe_std = comparison['wfe_std_filtered'] + wfe_opd = comparison['wfe_opd_filtered'] + + # Downsample if needed + n = len(x) + if n > downsample: + rng = np.random.default_rng(42) + sel = rng.choice(n, size=downsample, replace=False) + x, y = x[sel], y[sel] + lateral = lateral[sel] + wfe_std = wfe_std[sel] + wfe_opd = wfe_opd[sel] + + # Get recommendation + severity, recommendation = self._get_recommendation(comparison) + severity_color = { + 'CRITICAL': '#ef4444', + 'RECOMMENDED': '#f59e0b', + 'OPTIONAL': '#3b82f6', + 'EQUIVALENT': '#10b981' + }.get(severity, '#6b7280') + + # Create subplots: 2 rows, 2 cols + # Row 1: Lateral displacement map, WFE difference map + # Row 2: Comparison table, Recommendation + fig = _make_subplots( + rows=2, cols=2, + specs=[ + [{"type": "scene"}, {"type": "scene"}], + [{"type": "table", "colspan": 2}, None] + ], + row_heights=[0.65, 0.35], + column_widths=[0.5, 0.5], + vertical_spacing=0.08, + horizontal_spacing=0.05, + subplot_titles=[ + "Lateral Displacement Magnitude (µm)", + "Filtered WFE Comparison (nm)", + "Method Comparison" + ] + ) + + # === Plot 1: Lateral displacement map === + try: + tri = _Triangulation(x, y) + if tri.triangles is not None and len(tri.triangles) > 0: + i, j, k = tri.triangles.T + fig.add_trace(_go.Mesh3d( + x=x, y=y, z=lateral * amp * 100, # Scale for visibility + i=i, j=j, k=k, + intensity=lateral, + colorscale=colorscale_lateral, + opacity=1.0, + flatshading=False, + lighting=dict(ambient=0.4, diffuse=0.8, specular=0.3), + showscale=True, + colorbar=dict( + title=dict(text="Lateral (µm)", side="right"), + x=0.42, thickness=15, len=0.5 + ), + hovertemplate="X: %{x:.1f}
Y: %{y:.1f}
Lateral: %{intensity:.3f} µm" + ), row=1, col=1) + except: + fig.add_trace(_go.Scatter3d( + x=x, y=y, z=lateral * amp * 100, + mode='markers', + marker=dict(size=2, color=lateral, colorscale=colorscale_lateral, showscale=True) + ), row=1, col=1) + + # === Plot 2: WFE comparison (show OPD method) === + try: + tri = _Triangulation(x, y) + if tri.triangles is not None and len(tri.triangles) > 0: + i, j, k = tri.triangles.T + fig.add_trace(_go.Mesh3d( + x=x, y=y, z=wfe_opd * amp, + i=i, j=j, k=k, + intensity=wfe_opd, + colorscale='Turbo', + opacity=1.0, + flatshading=False, + lighting=dict(ambient=0.4, diffuse=0.8, specular=0.3), + showscale=True, + colorbar=dict( + title=dict(text="WFE (nm)", side="right"), + x=1.0, thickness=15, len=0.5 + ), + hovertemplate="X: %{x:.1f}
Y: %{y:.1f}
WFE: %{intensity:.2f} nm" + ), row=1, col=2) + except: + fig.add_trace(_go.Scatter3d( + x=x, y=y, z=wfe_opd * amp, + mode='markers', + marker=dict(size=2, color=wfe_opd, colorscale='Turbo', showscale=True) + ), row=1, col=2) + + # Configure 3D scenes + for scene_name in ['scene', 'scene2']: + fig.update_layout(**{ + scene_name: dict( + camera=dict(eye=dict(x=1.2, y=1.2, z=0.8)), + xaxis=dict(title="X (mm)", showgrid=True, gridcolor='rgba(128,128,128,0.3)'), + yaxis=dict(title="Y (mm)", showgrid=True, gridcolor='rgba(128,128,128,0.3)'), + zaxis=dict(title="", showgrid=True, gridcolor='rgba(128,128,128,0.3)'), + aspectmode='data', + ) + }) + + # === Comparison Table === + delta_sign = "+" if comparison['delta_filtered'] > 0 else "" + + fig.add_trace(_go.Table( + header=dict( + values=["Metric", "Standard (Z-only)", "OPD (Rigorous)", "Difference"], + align="left", + fill_color='#1f2937', + font=dict(color='white', size=12), + height=30 + ), + cells=dict( + values=[ + # Metric names + ["Global RMS (nm)", "Filtered RMS (nm)", "Max Lateral Disp (µm)", + "RMS Lateral Disp (µm)", "P99 Lateral Disp (µm)", "Focal Length (mm)", + f"RECOMMENDATION: {severity}"], + # Standard method + [f"{comparison['rms_std_global']:.2f}", f"{comparison['rms_std_filtered']:.2f}", + "N/A", "N/A", "N/A", f"{comparison['focal_length']:.1f}", ""], + # OPD method + [f"{comparison['rms_opd_global']:.2f}", f"{comparison['rms_opd_filtered']:.2f}", + f"{comparison['max_lateral_um']:.3f}", f"{comparison['rms_lateral_um']:.3f}", + f"{comparison['p99_lateral_um']:.3f}", f"{comparison['focal_length']:.1f}", ""], + # Difference + [f"{delta_sign}{comparison['delta_global']:.2f}", + f"{delta_sign}{comparison['delta_filtered']:.2f} ({delta_sign}{comparison['pct_diff_filtered']:.1f}%)", + "-", "-", "-", "-", recommendation], + ], + align="left", + fill_color=[['#374151']*7, ['#374151']*7, ['#374151']*7, + ['#374151']*6 + [severity_color]], + font=dict(color='white', size=11), + height=25 + ) + ), row=2, col=1) + + # Layout + fig.update_layout( + width=1600, + height=1000, + margin=dict(t=80, b=20, l=20, r=20), + paper_bgcolor='#111827', + plot_bgcolor='#1f2937', + font=dict(color='white'), + title=dict( + text=f"Atomizer Zernike Method Comparison - {title}
" + f"{severity}: {recommendation[:80]}...", + x=0.5, + font=dict(size=18) + ) + ) + + return fig.to_html(include_plotlyjs='cdn', full_html=True) + + def _generate(self, config: InsightConfig) -> InsightResult: + """Generate method comparison visualization.""" + self._load_data() + + cfg = {**DEFAULT_CONFIG, **config.extra} + n_modes = cfg['n_modes'] + filter_orders = cfg['filter_orders'] + disp_unit = cfg['disp_unit'] + + # Find subcases + disps = self._displacements + available = list(disps.keys()) + + if not available: + return InsightResult(success=False, error="No subcases found in OP2") + + timestamp = datetime.now().strftime("%Y%m%d_%H%M%S") + output_dir = config.output_dir or self.insights_path + output_dir.mkdir(parents=True, exist_ok=True) + + html_files = [] + all_comparisons = {} + + # Process each subcase + for label in available: + data = self._build_full_data(label, disp_unit) + comparison = self._compute_opd_comparison(data, n_modes, filter_orders) + + title = f"Subcase {label}" + html = self._generate_html(title, comparison, cfg) + + path = output_dir / f"zernike_opd_comparison_{timestamp}_{label}.html" + path.write_text(html, encoding='utf-8') + html_files.append(path) + + severity, _ = self._get_recommendation(comparison) + all_comparisons[label] = { + 'rms_std_filtered': comparison['rms_std_filtered'], + 'rms_opd_filtered': comparison['rms_opd_filtered'], + 'pct_diff': comparison['pct_diff_filtered'], + 'max_lateral_um': comparison['max_lateral_um'], + 'severity': severity, + } + + # Determine overall recommendation + severities = [c['severity'] for c in all_comparisons.values()] + if 'CRITICAL' in severities: + overall = "CRITICAL: Use OPD method for all optimization" + elif 'RECOMMENDED' in severities: + overall = "RECOMMENDED: Switch to OPD method for better accuracy" + elif 'OPTIONAL' in severities: + overall = "OPTIONAL: OPD method provides minor improvement" + else: + overall = "EQUIVALENT: Standard method is fine for this study" + + return InsightResult( + success=True, + html_path=html_files[0], + summary={ + 'html_files': [str(p) for p in html_files], + 'comparisons': all_comparisons, + 'overall_recommendation': overall, + } + ) diff --git a/optimization_engine/insights/zernike_wfe.py b/optimization_engine/insights/zernike_wfe.py index f3507187..a49887a4 100644 --- a/optimization_engine/insights/zernike_wfe.py +++ b/optimization_engine/insights/zernike_wfe.py @@ -7,6 +7,10 @@ Zernike polynomial decomposition. Generates three views: - 60 deg vs 20 deg (operational tilt comparison) - 90 deg Manufacturing (absolute with optician workload metrics) +Supports two WFE computation methods: +- Standard: Z-displacement only at original (x, y) coordinates +- OPD: Accounts for lateral (X, Y) displacement via interpolation (RECOMMENDED) + Applicable to: Mirror optimization studies with multi-subcase gravity loads. """ @@ -19,6 +23,9 @@ from numpy.linalg import LinAlgError from .base import StudyInsight, InsightConfig, InsightResult, register_insight +# Lazy import for OPD extractor +_ZernikeOPDExtractor = None + # Lazy imports to avoid startup overhead _plotly_loaded = False _go = None @@ -30,18 +37,20 @@ _BDF = None def _load_dependencies(): """Lazy load heavy dependencies.""" - global _plotly_loaded, _go, _make_subplots, _Triangulation, _OP2, _BDF + global _plotly_loaded, _go, _make_subplots, _Triangulation, _OP2, _BDF, _ZernikeOPDExtractor if not _plotly_loaded: import plotly.graph_objects as go from plotly.subplots import make_subplots from matplotlib.tri import Triangulation from pyNastran.op2.op2 import OP2 from pyNastran.bdf.bdf import BDF + from ..extractors.extract_zernike_figure import ZernikeOPDExtractor _go = go _make_subplots = make_subplots _Triangulation = Triangulation _OP2 = OP2 _BDF = BDF + _ZernikeOPDExtractor = ZernikeOPDExtractor _plotly_loaded = True @@ -144,10 +153,10 @@ def compute_zernike_coeffs( # Configuration Defaults # ============================================================================ DEFAULT_CONFIG = { - 'n_modes': 50, + 'n_modes': 36, # Reduced from 50 for faster computation (covers through 7th order) 'amp': 0.5, # Visual deformation scale 'pancake': 3.0, # Z-axis range multiplier - 'plot_downsample': 10000, + 'plot_downsample': 8000, # Reduced from 10000 for faster rendering 'filter_low_orders': 4, # Piston, tip, tilt, defocus 'colorscale': 'Turbo', 'disp_unit': 'mm', @@ -170,6 +179,7 @@ class ZernikeWFEInsight(StudyInsight): insight_type = "zernike_wfe" name = "Zernike WFE Analysis" description = "3D wavefront error surface with Zernike decomposition" + category = "optical" applicable_to = ["mirror", "optics", "wfe"] required_files = ["*.op2"] @@ -179,29 +189,59 @@ class ZernikeWFEInsight(StudyInsight): self.geo_path: Optional[Path] = None self._node_geo: Optional[Dict] = None self._displacements: Optional[Dict] = None + self._opd_extractor: Optional[Any] = None # Cached OPD extractor def can_generate(self) -> bool: - """Check if OP2 and geometry files exist.""" - # Look for OP2 in results or iterations - search_paths = [ - self.results_path, - self.study_path / "2_iterations", - self.setup_path / "model", + """Check if OP2 and geometry files exist. + + Uses fast non-recursive search to avoid slow glob operations. + """ + # Fast search order: best_design_archive first, then iterations + search_locations = [ + (self.study_path / "3_results" / "best_design_archive", False), # (path, is_iter_parent) + (self.study_path / "2_iterations", True), # iterations have subdirs + (self.setup_path / "model", False), ] - for search_path in search_paths: + op2_candidates = [] + + for search_path, is_iter_parent in search_locations: if not search_path.exists(): continue - op2_files = list(search_path.glob("**/*solution*.op2")) - if not op2_files: - op2_files = list(search_path.glob("**/*.op2")) - if op2_files: - self.op2_path = max(op2_files, key=lambda p: p.stat().st_mtime) - break - if self.op2_path is None: + if is_iter_parent: + # For iterations, check each subdir (one level only) + try: + for subdir in search_path.iterdir(): + if subdir.is_dir(): + for f in subdir.iterdir(): + if f.suffix.lower() == '.op2': + op2_candidates.append(f) + except (PermissionError, OSError): + continue + else: + # Direct check (non-recursive) + try: + for f in search_path.iterdir(): + if f.suffix.lower() == '.op2': + op2_candidates.append(f) + # Also check one level down for best_design_archive + elif f.is_dir(): + for sub in f.iterdir(): + if sub.suffix.lower() == '.op2': + op2_candidates.append(sub) + except (PermissionError, OSError): + continue + + if op2_candidates: + break # Found some, stop searching + + if not op2_candidates: return False + # Pick newest OP2 + self.op2_path = max(op2_candidates, key=lambda p: p.stat().st_mtime) + # Find geometry try: self.geo_path = self._find_geometry_file(self.op2_path) @@ -262,12 +302,16 @@ class ZernikeWFEInsight(StudyInsight): 'disp': dmat.copy() } - def _build_wfe_arrays( + def _build_wfe_arrays_standard( self, label: str, disp_unit: str = 'mm' - ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: - """Build X, Y, WFE arrays for a subcase.""" + ) -> Dict[str, np.ndarray]: + """Build X, Y, WFE arrays using standard Z-only method. + + This is the original method that uses only Z-displacement + at the original (x, y) coordinates. + """ nm_per_unit = 1e6 if disp_unit == 'mm' else 1e9 data = self._displacements[label] @@ -276,6 +320,8 @@ class ZernikeWFEInsight(StudyInsight): X, Y, WFE = [], [], [] valid_nids = [] + dx_arr, dy_arr, dz_arr = [], [], [] + for nid, vec in zip(node_ids, dmat): geo = self._node_geo.get(int(nid)) if geo is None: @@ -285,27 +331,148 @@ class ZernikeWFEInsight(StudyInsight): wfe = vec[2] * 2.0 * nm_per_unit # Z-disp to WFE WFE.append(wfe) valid_nids.append(nid) + dx_arr.append(vec[0]) + dy_arr.append(vec[1]) + dz_arr.append(vec[2]) - return (np.array(X), np.array(Y), np.array(WFE), np.array(valid_nids)) + return { + 'X': np.array(X), + 'Y': np.array(Y), + 'WFE': np.array(WFE), + 'node_ids': np.array(valid_nids), + 'dx': np.array(dx_arr), + 'dy': np.array(dy_arr), + 'dz': np.array(dz_arr), + 'lateral_disp': np.sqrt(np.array(dx_arr)**2 + np.array(dy_arr)**2), + 'method': 'standard', + } + + def _build_wfe_arrays_opd( + self, + label: str, + disp_unit: str = 'mm' + ) -> Dict[str, np.ndarray]: + """Build X, Y, WFE arrays using OPD method (accounts for lateral displacement). + + This is the RECOMMENDED method that: + 1. Uses deformed (x+dx, y+dy) coordinates for Zernike fitting + 2. Computes true surface error via interpolation of undeformed geometry + + Uses cached OPD extractor to avoid re-reading OP2/BDF files for each subcase. + """ + _load_dependencies() + + # Reuse cached extractor to avoid re-reading files + if self._opd_extractor is None: + self._opd_extractor = _ZernikeOPDExtractor( + self.op2_path, + bdf_path=self.geo_path, + displacement_unit=disp_unit + ) + + opd_data = self._opd_extractor._build_figure_opd_data(label) + + return { + 'X': opd_data['x_deformed'], + 'Y': opd_data['y_deformed'], + 'WFE': opd_data['wfe_nm'], + 'node_ids': opd_data['node_ids'], + 'dx': opd_data['dx'], + 'dy': opd_data['dy'], + 'dz': opd_data['dz'], + 'lateral_disp': opd_data['lateral_disp'], + 'x_original': opd_data['x_original'], + 'y_original': opd_data['y_original'], + 'method': 'opd', + } + + def _build_wfe_arrays( + self, + label: str, + disp_unit: str = 'mm', + method: str = 'opd' + ) -> Dict[str, np.ndarray]: + """Build X, Y, WFE arrays for a subcase. + + Args: + label: Subcase label + disp_unit: Displacement unit ('mm' or 'm') + method: 'opd' (recommended) or 'standard' + + Returns: + Dict with X, Y, WFE, node_ids, dx, dy, dz, lateral_disp arrays + """ + if method == 'opd': + return self._build_wfe_arrays_opd(label, disp_unit) + else: + return self._build_wfe_arrays_standard(label, disp_unit) def _compute_relative_wfe( self, - X1, Y1, WFE1, nids1, - X2, Y2, WFE2, nids2 - ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: - """Compute WFE1 - WFE2 for common nodes.""" + data1: Dict[str, np.ndarray], + data2: Dict[str, np.ndarray] + ) -> Dict[str, np.ndarray]: + """Compute relative displacement and WFE (target - reference) for common nodes. + + Args: + data1: Target subcase data dict + data2: Reference subcase data dict + + Returns: + Dict with relative WFE arrays AND relative displacement components (dx, dy, dz) + """ + X1, Y1, WFE1, nids1 = data1['X'], data1['Y'], data1['WFE'], data1['node_ids'] + X2, Y2, WFE2, nids2 = data2['X'], data2['Y'], data2['WFE'], data2['node_ids'] + + # Build reference maps for all data ref_map = {int(nid): (x, y, w) for nid, x, y, w in zip(nids2, X2, Y2, WFE2)} - X_rel, Y_rel, WFE_rel = [], [], [] + # Maps for displacement components + dx1_map = {int(nid): d for nid, d in zip(data1['node_ids'], data1['dx'])} + dy1_map = {int(nid): d for nid, d in zip(data1['node_ids'], data1['dy'])} + dz1_map = {int(nid): d for nid, d in zip(data1['node_ids'], data1['dz'])} + dx2_map = {int(nid): d for nid, d in zip(data2['node_ids'], data2['dx'])} + dy2_map = {int(nid): d for nid, d in zip(data2['node_ids'], data2['dy'])} + dz2_map = {int(nid): d for nid, d in zip(data2['node_ids'], data2['dz'])} + + X_rel, Y_rel, WFE_rel, nids_rel = [], [], [], [] + dx_rel, dy_rel, dz_rel = [], [], [] + lateral_rel = [] + for nid, x, y, w in zip(nids1, X1, Y1, WFE1): nid = int(nid) - if nid in ref_map: - _, _, w_ref = ref_map[nid] - X_rel.append(x) - Y_rel.append(y) - WFE_rel.append(w - w_ref) + if nid not in ref_map: + continue + if nid not in dx2_map: + continue - return np.array(X_rel), np.array(Y_rel), np.array(WFE_rel) + _, _, w_ref = ref_map[nid] + X_rel.append(x) + Y_rel.append(y) + WFE_rel.append(w - w_ref) + nids_rel.append(nid) + + # Compute relative displacements (target - reference) + dx_rel.append(dx1_map[nid] - dx2_map[nid]) + dy_rel.append(dy1_map[nid] - dy2_map[nid]) + dz_rel.append(dz1_map[nid] - dz2_map[nid]) + + # Relative lateral displacement magnitude + lat1 = np.sqrt(dx1_map[nid]**2 + dy1_map[nid]**2) + lat2 = np.sqrt(dx2_map[nid]**2 + dy2_map[nid]**2) + lateral_rel.append(lat1 - lat2) + + return { + 'X': np.array(X_rel), + 'Y': np.array(Y_rel), + 'WFE': np.array(WFE_rel), + 'node_ids': np.array(nids_rel), + 'dx': np.array(dx_rel), # Relative X displacement (mm) + 'dy': np.array(dy_rel), # Relative Y displacement (mm) + 'dz': np.array(dz_rel), # Relative Z displacement (mm) + 'lateral_disp': np.array(lateral_rel) if lateral_rel else np.zeros(len(X_rel)), + 'method': data1.get('method', 'unknown'), + } def _compute_metrics( self, @@ -585,8 +752,430 @@ class ZernikeWFEInsight(StudyInsight): return fig.to_html(include_plotlyjs='cdn', full_html=True) + def _generate_lateral_map_html( + self, + title: str, + data: Dict[str, np.ndarray], + config: Dict, + ) -> str: + """Generate HTML for lateral displacement visualization. + + Shows a 3D surface colored by lateral displacement magnitude, + with metrics table showing max/RMS/mean lateral displacement. + """ + _load_dependencies() + + X = data['X'] + Y = data['Y'] + lateral_disp = data['lateral_disp'] # in mm + downsample = config.get('plot_downsample', 10000) + + # Convert to µm for display + lateral_um = lateral_disp * 1000.0 # mm to µm + + # Downsample + n = len(X) + if n > downsample: + rng = np.random.default_rng(42) + sel = rng.choice(n, size=downsample, replace=False) + Xp, Yp, Lp = X[sel], Y[sel], lateral_um[sel] + else: + Xp, Yp, Lp = X, Y, lateral_um + + # Build mesh + mesh_traces = [] + try: + tri = _Triangulation(Xp, Yp) + if tri.triangles is not None and len(tri.triangles) > 0: + i, j, k = tri.triangles.T + mesh_traces.append(_go.Mesh3d( + x=Xp, y=Yp, z=Lp, + i=i, j=j, k=k, + intensity=Lp, + colorscale='Viridis', + opacity=1.0, + flatshading=False, + lighting=dict(ambient=0.4, diffuse=0.8, specular=0.3), + lightposition=dict(x=100, y=200, z=300), + showscale=True, + colorbar=dict(title=dict(text="Lateral (µm)", side="right"), + thickness=15, len=0.6, tickformat=".3f"), + hovertemplate="X: %{x:.1f}
Y: %{y:.1f}
Lateral: %{z:.4f} µm" + )) + except Exception: + pass + + if not mesh_traces: + mesh_traces.append(_go.Scatter3d( + x=Xp, y=Yp, z=Lp, + mode='markers', + marker=dict(size=2, color=Lp, colorscale='Viridis', showscale=True), + showlegend=False + )) + + # Create figure with subplots + fig = _make_subplots( + rows=2, cols=1, + specs=[[{"type": "scene"}], [{"type": "table"}]], + row_heights=[0.75, 0.25], + vertical_spacing=0.03, + subplot_titles=[ + f"Lateral Displacement Map - {title}", + "Lateral Displacement Statistics" + ] + ) + + for tr in mesh_traces: + fig.add_trace(tr, row=1, col=1) + + # Stats + max_lat = float(np.max(lateral_um)) + rms_lat = float(np.sqrt(np.mean(lateral_um**2))) + mean_lat = float(np.mean(lateral_um)) + min_lat = float(np.min(lateral_um)) + + fig.add_trace(_go.Table( + header=dict(values=["Statistic", "Value (µm)"], + align="left", fill_color='#1f2937', font=dict(color='white')), + cells=dict(values=[ + ["Max Lateral", "RMS Lateral", "Mean Lateral", "Min Lateral"], + [f"{max_lat:.4f}", f"{rms_lat:.4f}", f"{mean_lat:.4f}", f"{min_lat:.4f}"] + ], align="left", fill_color='#374151', font=dict(color='white')) + ), row=2, col=1) + + # Configure 3D scene + max_z = float(np.max(Lp)) if Lp.size else 1.0 + fig.update_scenes( + camera=dict(eye=dict(x=1.2, y=1.2, z=0.8), up=dict(x=0, y=0, z=1)), + xaxis=dict(title="X (mm)", showgrid=True, gridcolor='rgba(128,128,128,0.3)'), + yaxis=dict(title="Y (mm)", showgrid=True, gridcolor='rgba(128,128,128,0.3)'), + zaxis=dict(title="Lateral (µm)", + range=[0, max_z * 1.2], + showgrid=True, gridcolor='rgba(128,128,128,0.3)'), + aspectmode='manual', + aspectratio=dict(x=1, y=1, z=0.4) + ) + + fig.update_layout( + width=1200, height=900, + margin=dict(t=60, b=20, l=20, r=20), + paper_bgcolor='#111827', plot_bgcolor='#1f2937', + font=dict(color='white'), + title=dict(text=f"Lateral Displacement Analysis - {title}", + x=0.5, font=dict(size=18)) + ) + + return fig.to_html(include_plotlyjs='cdn', full_html=True) + + def _generate_dual_method_view_html( + self, + title: str, + data_std: Dict[str, np.ndarray], + data_opd: Dict[str, np.ndarray], + rms_std: Dict, + rms_opd: Dict, + config: Dict, + is_relative: bool = False, + ref_title: str = "20 deg", + ) -> str: + """Generate HTML with toggle between Standard/OPD methods AND X/Y/Z displacement components. + + For relative views, provides toggles to see: + - WFE (Z): The main wavefront error view (default) + - ΔX: Relative X displacement between subcases + - ΔY: Relative Y displacement between subcases + - ΔZ: Relative Z displacement between subcases + """ + _load_dependencies() + + n_modes = config.get('n_modes', 50) + amp = config.get('amp', 0.5) + pancake = config.get('pancake', 3.0) + downsample = config.get('plot_downsample', 10000) + colorscale = config.get('colorscale', 'Turbo') + + title_suffix = f" (relative to {ref_title})" if is_relative else " (absolute)" + + # Build traces for both methods (WFE view) + traces_std_wfe = [] + traces_opd_wfe = [] + + # Build displacement component traces (OPD method only, for relative views) + traces_dx = [] + traces_dy = [] + traces_dz = [] + + # Helper to build mesh trace + def build_mesh_trace(Xp, Yp, Zp, colorscale, label, unit, colorbar_title): + try: + tri = _Triangulation(Xp, Yp) + if tri.triangles is not None and len(tri.triangles) > 0: + i, j, k = tri.triangles.T + return _go.Mesh3d( + x=Xp, y=Yp, z=Zp, + i=i, j=j, k=k, + intensity=Zp, + colorscale=colorscale, + opacity=1.0, + flatshading=False, + lighting=dict(ambient=0.4, diffuse=0.8, specular=0.3), + lightposition=dict(x=100, y=200, z=300), + showscale=True, + colorbar=dict(title=dict(text=colorbar_title, side="right"), + thickness=15, len=0.6, tickformat=".2f" if 'µm' in unit else ".1f"), + hovertemplate=f"{label}
X: %{{x:.1f}}
Y: %{{y:.1f}}
Value: %{{z:.3f}} {unit}" + ) + except Exception: + pass + return _go.Scatter3d( + x=Xp, y=Yp, z=Zp, + mode='markers', + marker=dict(size=2, color=Zp, colorscale=colorscale, showscale=True), + showlegend=False + ) + + # Build WFE traces for both methods + for data, rms_data, traces, method_name in [ + (data_std, rms_std, traces_std_wfe, 'Standard'), + (data_opd, rms_opd, traces_opd_wfe, 'OPD') + ]: + X, Y = data['X'], data['Y'] + W_res_filt = rms_data['W_res_filt'] + + # Downsample + n = len(X) + if n > downsample: + rng = np.random.default_rng(42) + sel = rng.choice(n, size=downsample, replace=False) + Xp, Yp, Wp = X[sel], Y[sel], W_res_filt[sel] + else: + Xp, Yp, Wp = X, Y, W_res_filt + + res_amp = amp * Wp + traces.append(build_mesh_trace(Xp, Yp, res_amp, colorscale, method_name, 'nm', f'{method_name} WFE (nm)')) + + # Build displacement component traces (for relative views) + has_displacement_data = is_relative and 'dx' in data_opd and len(data_opd.get('dx', [])) > 0 + if has_displacement_data: + X, Y = data_opd['X'], data_opd['Y'] + dx_mm = data_opd['dx'] # mm + dy_mm = data_opd['dy'] # mm + dz_mm = data_opd['dz'] # mm + + # Convert to µm for display + dx_um = dx_mm * 1000.0 + dy_um = dy_mm * 1000.0 + dz_um = dz_mm * 1000.0 + + n = len(X) + if n > downsample: + rng = np.random.default_rng(42) + sel = rng.choice(n, size=downsample, replace=False) + Xp, Yp = X[sel], Y[sel] + dxp, dyp, dzp = dx_um[sel], dy_um[sel], dz_um[sel] + else: + Xp, Yp = X, Y + dxp, dyp, dzp = dx_um, dy_um, dz_um + + # Apply visual amplification for displacement views + disp_amp = amp * 1000.0 # Scale factor for µm display + + traces_dx.append(build_mesh_trace(Xp, Yp, dxp * amp, 'RdBu_r', 'ΔX Displacement', 'µm', 'ΔX (µm)')) + traces_dy.append(build_mesh_trace(Xp, Yp, dyp * amp, 'RdBu_r', 'ΔY Displacement', 'µm', 'ΔY (µm)')) + traces_dz.append(build_mesh_trace(Xp, Yp, dzp * amp, 'RdBu_r', 'ΔZ Displacement', 'µm', 'ΔZ (µm)')) + + # Create figure + fig = _make_subplots( + rows=2, cols=1, + specs=[[{"type": "scene"}], [{"type": "table"}]], + row_heights=[0.65, 0.35], + vertical_spacing=0.05, + subplot_titles=[ + f"Surface Analysis - {title}{title_suffix}", + "Metrics Comparison" + ] + ) + + # Add all traces in order: [std_wfe, opd_wfe, dx, dy, dz, table] + # Start with OPD WFE visible (default view) + for tr in traces_std_wfe: + tr.visible = False + fig.add_trace(tr, row=1, col=1) + for tr in traces_opd_wfe: + tr.visible = True # Default view + fig.add_trace(tr, row=1, col=1) + for tr in traces_dx: + tr.visible = False + fig.add_trace(tr, row=1, col=1) + for tr in traces_dy: + tr.visible = False + fig.add_trace(tr, row=1, col=1) + for tr in traces_dz: + tr.visible = False + fig.add_trace(tr, row=1, col=1) + + # Compute lateral stats (from OPD data) + lateral_um = data_opd.get('lateral_disp', np.zeros(1)) * 1000.0 + max_lat = float(np.max(np.abs(lateral_um))) + rms_lat = float(np.sqrt(np.mean(lateral_um**2))) + + # Compute % difference + std_filt = rms_std['filtered_rms'] + opd_filt = rms_opd['filtered_rms'] + pct_diff = 100.0 * (opd_filt - std_filt) / std_filt if std_filt > 0 else 0.0 + + # Displacement stats (for relative views) + disp_stats_rows = [] + if has_displacement_data: + dx_um = data_opd['dx'] * 1000.0 + dy_um = data_opd['dy'] * 1000.0 + dz_um = data_opd['dz'] * 1000.0 + disp_stats_rows = [ + "ΔX RMS (µm)", "ΔY RMS (µm)", "ΔZ RMS (µm)" + ] + disp_stats_values_std = ["—", "—", "—"] + disp_stats_values_opd = [ + f"{float(np.sqrt(np.mean(dx_um**2))):.4f}", + f"{float(np.sqrt(np.mean(dy_um**2))):.4f}", + f"{float(np.sqrt(np.mean(dz_um**2))):.4f}" + ] + disp_stats_diff = ["—", "—", "—"] + + # Comparison table + table_headers = ["Metric", "Standard (Z-only)", "OPD (X,Y,Z)", "Difference"] + table_rows = ["Global RMS (nm)", "Filtered RMS (nm)", "Method", "Max Lateral (µm)", "RMS Lateral (µm)"] + table_std = [f"{rms_std['global_rms']:.2f}", f"{std_filt:.2f}", "Z-displacement only", "—", "—"] + table_opd = [f"{rms_opd['global_rms']:.2f}", f"{opd_filt:.2f}", "Deformed coords + OPD", f"{max_lat:.3f}", f"{rms_lat:.3f}"] + table_diff = ["—", f"{pct_diff:+.2f}%", "← RECOMMENDED", "—", "—"] + + if has_displacement_data: + table_rows.extend(disp_stats_rows) + table_std.extend(disp_stats_values_std) + table_opd.extend(disp_stats_values_opd) + table_diff.extend(disp_stats_diff) + + fig.add_trace(_go.Table( + header=dict(values=table_headers, align="left", fill_color='#1f2937', font=dict(color='white')), + cells=dict(values=[table_rows, table_std, table_opd, table_diff], + align="left", fill_color='#374151', font=dict(color='white')) + ), row=2, col=1) + + # Build visibility arrays for toggles + n_std_wfe = len(traces_std_wfe) + n_opd_wfe = len(traces_opd_wfe) + n_dx = len(traces_dx) + n_dy = len(traces_dy) + n_dz = len(traces_dz) + # Total traces before table: n_std_wfe + n_opd_wfe + n_dx + n_dy + n_dz + # Then table is last + + def make_visibility(show_std_wfe=False, show_opd_wfe=False, show_dx=False, show_dy=False, show_dz=False): + vis = [] + vis.extend([show_std_wfe] * n_std_wfe) + vis.extend([show_opd_wfe] * n_opd_wfe) + vis.extend([show_dx] * n_dx) + vis.extend([show_dy] * n_dy) + vis.extend([show_dz] * n_dz) + vis.append(True) # Table always visible + return vis + + # Build button definitions + buttons_method = [ + dict(label="OPD Method (Recommended)", + method="update", + args=[{"visible": make_visibility(show_opd_wfe=True)}]), + dict(label="Standard Method (Z-only)", + method="update", + args=[{"visible": make_visibility(show_std_wfe=True)}]), + ] + + buttons_component = [ + dict(label="WFE (Z)", + method="update", + args=[{"visible": make_visibility(show_opd_wfe=True)}]), + ] + + if has_displacement_data: + buttons_component.extend([ + dict(label="ΔX Disp", + method="update", + args=[{"visible": make_visibility(show_dx=True)}]), + dict(label="ΔY Disp", + method="update", + args=[{"visible": make_visibility(show_dy=True)}]), + dict(label="ΔZ Disp", + method="update", + args=[{"visible": make_visibility(show_dz=True)}]), + ]) + + # Create update menus + updatemenus = [ + dict( + type="buttons", + direction="right", + x=0.0, y=1.15, + xanchor="left", + showactive=True, + buttons=buttons_method, + font=dict(size=11), + pad=dict(r=10, t=10), + ), + ] + + if has_displacement_data: + updatemenus.append( + dict( + type="buttons", + direction="right", + x=0.55, y=1.15, + xanchor="left", + showactive=True, + buttons=buttons_component, + font=dict(size=11), + pad=dict(r=10, t=10), + ) + ) + + fig.update_layout(updatemenus=updatemenus) + + # Configure 3D scene + max_amp_opd = float(np.max(np.abs(amp * rms_opd['W_res_filt']))) if rms_opd['W_res_filt'].size else 1.0 + fig.update_scenes( + camera=dict(eye=dict(x=1.2, y=1.2, z=0.8), up=dict(x=0, y=0, z=1)), + xaxis=dict(title="X (mm)", showgrid=True, gridcolor='rgba(128,128,128,0.3)'), + yaxis=dict(title="Y (mm)", showgrid=True, gridcolor='rgba(128,128,128,0.3)'), + zaxis=dict(title="Value", + range=[-max_amp_opd * pancake, max_amp_opd * pancake], + showgrid=True, gridcolor='rgba(128,128,128,0.3)'), + aspectmode='manual', + aspectratio=dict(x=1, y=1, z=0.4) + ) + + fig.update_layout( + width=1400, height=1100, + margin=dict(t=100, b=20, l=20, r=20), + paper_bgcolor='#111827', plot_bgcolor='#1f2937', + font=dict(color='white'), + title=dict(text=f"Atomizer Zernike Analysis - {title}", + x=0.5, font=dict(size=18)), + annotations=[ + dict(text="Method:", x=0.0, y=1.18, xref="paper", yref="paper", + showarrow=False, font=dict(size=12, color='white'), xanchor='left'), + dict(text="View:", x=0.55, y=1.18, xref="paper", yref="paper", + showarrow=False, font=dict(size=12, color='white'), xanchor='left') if has_displacement_data else {}, + ] + ) + + return fig.to_html(include_plotlyjs='cdn', full_html=True) + def _generate(self, config: InsightConfig) -> InsightResult: - """Generate all Zernike WFE views.""" + """Generate all Zernike WFE views with Standard/OPD toggle and lateral maps. + + Performance optimizations: + - Uses cached OPD extractor (reads OP2/BDF only once) + - Loads all subcase data upfront to minimize I/O + - Standard method reuses geometry from already-loaded data + """ self._load_data() # Merge config @@ -626,66 +1215,110 @@ class ZernikeWFEInsight(StudyInsight): html_files = [] summary = {} + # Load data for BOTH methods - OPD method shares cached extractor + # Standard method uses already-loaded _node_geo and _displacements # Reference: 20 deg - X_ref, Y_ref, WFE_ref, nids_ref = self._build_wfe_arrays(sc_map['20'], disp_unit) - rms_ref = self._compute_metrics(X_ref, Y_ref, WFE_ref, n_modes, filter_orders) + data_ref_std = self._build_wfe_arrays_standard(sc_map['20'], disp_unit) + data_ref_opd = self._build_wfe_arrays_opd(sc_map['20'], disp_unit) + + # 40 deg + data_40_std = self._build_wfe_arrays_standard(sc_map['40'], disp_unit) + data_40_opd = self._build_wfe_arrays_opd(sc_map['40'], disp_unit) + + # 60 deg + data_60_std = self._build_wfe_arrays_standard(sc_map['60'], disp_unit) + data_60_opd = self._build_wfe_arrays_opd(sc_map['60'], disp_unit) # 90 deg - X_90, Y_90, WFE_90, nids_90 = self._build_wfe_arrays(sc_map['90'], disp_unit) - rms_90 = self._compute_metrics(X_90, Y_90, WFE_90, n_modes, filter_orders) - mfg_metrics = self._compute_aberration_magnitudes(rms_90['coefficients']) + data_90_std = self._build_wfe_arrays_standard(sc_map['90'], disp_unit) + data_90_opd = self._build_wfe_arrays_opd(sc_map['90'], disp_unit) - # 40 deg vs 20 deg - X_40, Y_40, WFE_40, nids_40 = self._build_wfe_arrays(sc_map['40'], disp_unit) - X_40_rel, Y_40_rel, WFE_40_rel = self._compute_relative_wfe( - X_40, Y_40, WFE_40, nids_40, X_ref, Y_ref, WFE_ref, nids_ref) - rms_40_abs = self._compute_metrics(X_40, Y_40, WFE_40, n_modes, filter_orders) - rms_40_rel = self._compute_metrics(X_40_rel, Y_40_rel, WFE_40_rel, n_modes, filter_orders) + # ========================================= + # 40 deg vs 20 deg (with dual method toggle) + # ========================================= + rel_40_std = self._compute_relative_wfe(data_40_std, data_ref_std) + rel_40_opd = self._compute_relative_wfe(data_40_opd, data_ref_opd) - html_40 = self._generate_view_html( - "40 deg", X_40_rel, Y_40_rel, WFE_40_rel, rms_40_rel, cfg, - is_relative=True, ref_title="20 deg", - abs_pair=(rms_40_abs['global_rms'], rms_40_abs['filtered_rms'])) + rms_40_std = self._compute_metrics(rel_40_std['X'], rel_40_std['Y'], rel_40_std['WFE'], + n_modes, filter_orders) + rms_40_opd = self._compute_metrics(rel_40_opd['X'], rel_40_opd['Y'], rel_40_opd['WFE'], + n_modes, filter_orders) + + html_40 = self._generate_dual_method_view_html( + "40 deg", rel_40_std, rel_40_opd, rms_40_std, rms_40_opd, cfg, + is_relative=True, ref_title="20 deg") path_40 = output_dir / f"zernike_{timestamp}_40_vs_20.html" path_40.write_text(html_40, encoding='utf-8') html_files.append(path_40) - summary['40_vs_20_filtered_rms'] = rms_40_rel['filtered_rms'] - # 60 deg vs 20 deg - X_60, Y_60, WFE_60, nids_60 = self._build_wfe_arrays(sc_map['60'], disp_unit) - X_60_rel, Y_60_rel, WFE_60_rel = self._compute_relative_wfe( - X_60, Y_60, WFE_60, nids_60, X_ref, Y_ref, WFE_ref, nids_ref) - rms_60_abs = self._compute_metrics(X_60, Y_60, WFE_60, n_modes, filter_orders) - rms_60_rel = self._compute_metrics(X_60_rel, Y_60_rel, WFE_60_rel, n_modes, filter_orders) + # Lateral map for 40 deg + html_40_lat = self._generate_lateral_map_html("40 deg", data_40_opd, cfg) + path_40_lat = output_dir / f"zernike_{timestamp}_40_lateral.html" + path_40_lat.write_text(html_40_lat, encoding='utf-8') + html_files.append(path_40_lat) - html_60 = self._generate_view_html( - "60 deg", X_60_rel, Y_60_rel, WFE_60_rel, rms_60_rel, cfg, - is_relative=True, ref_title="20 deg", - abs_pair=(rms_60_abs['global_rms'], rms_60_abs['filtered_rms'])) + summary['40_vs_20_filtered_rms_std'] = rms_40_std['filtered_rms'] + summary['40_vs_20_filtered_rms_opd'] = rms_40_opd['filtered_rms'] + + # ========================================= + # 60 deg vs 20 deg (with dual method toggle) + # ========================================= + rel_60_std = self._compute_relative_wfe(data_60_std, data_ref_std) + rel_60_opd = self._compute_relative_wfe(data_60_opd, data_ref_opd) + + rms_60_std = self._compute_metrics(rel_60_std['X'], rel_60_std['Y'], rel_60_std['WFE'], + n_modes, filter_orders) + rms_60_opd = self._compute_metrics(rel_60_opd['X'], rel_60_opd['Y'], rel_60_opd['WFE'], + n_modes, filter_orders) + + html_60 = self._generate_dual_method_view_html( + "60 deg", rel_60_std, rel_60_opd, rms_60_std, rms_60_opd, cfg, + is_relative=True, ref_title="20 deg") path_60 = output_dir / f"zernike_{timestamp}_60_vs_20.html" path_60.write_text(html_60, encoding='utf-8') html_files.append(path_60) - summary['60_vs_20_filtered_rms'] = rms_60_rel['filtered_rms'] - # 90 deg Manufacturing - X_90_rel, Y_90_rel, WFE_90_rel = self._compute_relative_wfe( - X_90, Y_90, WFE_90, nids_90, X_ref, Y_ref, WFE_ref, nids_ref) - rms_90_rel = self._compute_metrics(X_90_rel, Y_90_rel, WFE_90_rel, n_modes, filter_orders) - corr_abr = self._compute_aberration_magnitudes(rms_90_rel['coefficients']) - correction_metrics = { - 'rms_filter_j1to3': rms_90_rel['rms_filter_j1to3'], - **corr_abr - } + # Lateral map for 60 deg + html_60_lat = self._generate_lateral_map_html("60 deg", data_60_opd, cfg) + path_60_lat = output_dir / f"zernike_{timestamp}_60_lateral.html" + path_60_lat.write_text(html_60_lat, encoding='utf-8') + html_files.append(path_60_lat) - html_90 = self._generate_view_html( - "90 deg (Manufacturing)", X_90, Y_90, WFE_90, rms_90, cfg, - is_relative=False, is_manufacturing=True, - mfg_metrics=mfg_metrics, correction_metrics=correction_metrics) + summary['60_vs_20_filtered_rms_std'] = rms_60_std['filtered_rms'] + summary['60_vs_20_filtered_rms_opd'] = rms_60_opd['filtered_rms'] + + # ========================================= + # 90 deg Manufacturing (absolute, with dual method toggle) + # ========================================= + rms_90_std = self._compute_metrics(data_90_std['X'], data_90_std['Y'], data_90_std['WFE'], + n_modes, filter_orders) + rms_90_opd = self._compute_metrics(data_90_opd['X'], data_90_opd['Y'], data_90_opd['WFE'], + n_modes, filter_orders) + + html_90 = self._generate_dual_method_view_html( + "90 deg (Manufacturing)", data_90_std, data_90_opd, rms_90_std, rms_90_opd, cfg, + is_relative=False) path_90 = output_dir / f"zernike_{timestamp}_90_mfg.html" path_90.write_text(html_90, encoding='utf-8') html_files.append(path_90) - summary['90_mfg_filtered_rms'] = rms_90['filtered_rms'] - summary['90_optician_workload'] = rms_90['rms_filter_j1to3'] + + # Lateral map for 90 deg + html_90_lat = self._generate_lateral_map_html("90 deg (Manufacturing)", data_90_opd, cfg) + path_90_lat = output_dir / f"zernike_{timestamp}_90_mfg_lateral.html" + path_90_lat.write_text(html_90_lat, encoding='utf-8') + html_files.append(path_90_lat) + + summary['90_mfg_filtered_rms_std'] = rms_90_std['filtered_rms'] + summary['90_mfg_filtered_rms_opd'] = rms_90_opd['filtered_rms'] + summary['90_optician_workload'] = rms_90_opd['rms_filter_j1to3'] + + # Lateral displacement summary + lateral_40 = data_40_opd.get('lateral_disp', np.zeros(1)) * 1000.0 # mm to µm + lateral_60 = data_60_opd.get('lateral_disp', np.zeros(1)) * 1000.0 + lateral_90 = data_90_opd.get('lateral_disp', np.zeros(1)) * 1000.0 + summary['lateral_40_max_um'] = float(np.max(lateral_40)) + summary['lateral_60_max_um'] = float(np.max(lateral_60)) + summary['lateral_90_max_um'] = float(np.max(lateral_90)) return InsightResult( success=True, diff --git a/tests/analyze_v11.py b/tests/analyze_v11.py new file mode 100644 index 00000000..b0841422 --- /dev/null +++ b/tests/analyze_v11.py @@ -0,0 +1,150 @@ +"""Analyze V11 optimization results.""" +import sqlite3 +from pathlib import Path +import json +import sys +sys.path.insert(0, '.') + +db_path = Path('studies/M1_Mirror/m1_mirror_cost_reduction_V11/3_results/study.db') + +conn = sqlite3.connect(db_path) +c = conn.cursor() + +# Get all completed trials with their objectives +c.execute(''' + SELECT t.trial_id + FROM trials t + WHERE t.state = 'COMPLETE' + ORDER BY t.trial_id +''') +completed_ids = [row[0] for row in c.fetchall()] +print(f'Completed trials: {len(completed_ids)}') + +# Build trial data +trials = [] +for tid in completed_ids: + c.execute(''' + SELECT key, value_json + FROM trial_user_attributes + WHERE trial_id = ? + ''', (tid,)) + attrs = {row[0]: json.loads(row[1]) for row in c.fetchall()} + + if 'rel_filtered_rms_40_vs_20' in attrs: + trials.append({ + 'id': tid, + 'rms_40': attrs.get('rel_filtered_rms_40_vs_20', 999), + 'rms_60': attrs.get('rel_filtered_rms_60_vs_20', 999), + 'mfg_90': attrs.get('mfg_90_optician_workload', 999), + 'ws': attrs.get('weighted_sum', 999), + 'lateral': attrs.get('lateral_rms_um', 0), + }) + +# Sort by weighted sum or calculate it +for t in trials: + if t['ws'] == 999: + # Calculate weighted sum + w40 = 100 / 4.0 # Target 4 nm + w60 = 50 / 10.0 # Target 10 nm + w90 = 20 / 20.0 # Target 20 nm + t['ws'] = w40 * t['rms_40'] + w60 * t['rms_60'] + w90 * t['mfg_90'] + +# Find best +trials.sort(key=lambda x: x['ws']) +best = trials[0] if trials else None + +print('\nV11 Optimization Summary') +print('=' * 70) +print(f'Completed trials: {len(trials)}') + +if trials: + rms40 = [t['rms_40'] for t in trials] + rms60 = [t['rms_60'] for t in trials] + mfg90 = [t['mfg_90'] for t in trials] + + print(f'\n40-20 RMS range: {min(rms40):.2f} - {max(rms40):.2f} nm (target: 4.0 nm)') + print(f'60-20 RMS range: {min(rms60):.2f} - {max(rms60):.2f} nm (target: 10.0 nm)') + print(f'90 MFG range: {min(mfg90):.2f} - {max(mfg90):.2f} nm (target: 20.0 nm)') + + print(f'\nBest Trial (#{best["id"]}):') + print(f' 40-20 RMS: {best["rms_40"]:.2f} nm') + print(f' 60-20 RMS: {best["rms_60"]:.2f} nm') + print(f' MFG 90: {best["mfg_90"]:.2f} nm') + print(f' Lateral: {best["lateral"]:.3f} um') + print(f' WS: {best["ws"]:.1f}') + + # Check if values make sense (not too good to be true) + print('\nSanity Check:') + if best['rms_40'] < 3.0: + print(' WARNING: 40-20 RMS suspiciously low!') + else: + print(f' 40-20 RMS {best["rms_40"]:.2f} nm - looks reasonable (expected ~6-7nm)') + + if best['rms_60'] < 8.0: + print(' WARNING: 60-20 RMS suspiciously low!') + else: + print(f' 60-20 RMS {best["rms_60"]:.2f} nm - looks reasonable (expected ~13-15nm)') + + # Compare to V7 baseline and V10 buggy values + print('\n' + '=' * 70) + print('Comparison to Known Values:') + print('=' * 70) + print('\nV7 Baseline (original Zernike, correct):') + print(' 40-20 RMS: 6.05 nm') + print(' 60-20 RMS: 13.03 nm') + print(' MFG 90: 26.34 nm') + + print('\nV10 BUGGY (abs(RMS_target - RMS_ref) - WRONG):') + print(' 40-20 RMS: ~1.99 nm <- WAY TOO LOW') + print(' 60-20 RMS: ~6.82 nm <- WAY TOO LOW') + + print('\nV11 Best (extract_relative() - CORRECT):') + print(f' 40-20 RMS: {best["rms_40"]:.2f} nm') + print(f' 60-20 RMS: {best["rms_60"]:.2f} nm') + print(f' MFG 90: {best["mfg_90"]:.2f} nm') + + # Verdict + print('\n' + '=' * 70) + print('VERDICT:') + print('=' * 70) + if best['rms_40'] > 5.0 and best['rms_40'] < 10.0: + print(' 40-20: CORRECT - matches expected V7 range') + else: + print(f' 40-20: INVESTIGATE - {best["rms_40"]:.2f} nm is outside expected range') + + if best['rms_60'] > 10.0 and best['rms_60'] < 20.0: + print(' 60-20: CORRECT - matches expected V7 range') + else: + print(f' 60-20: INVESTIGATE - {best["rms_60"]:.2f} nm is outside expected range') + +conn.close() + +# Now generate Zernike dashboard for visual verification +print('\n' + '=' * 70) +print('Generating Zernike Dashboard for Best Iteration...') +print('=' * 70) + +from optimization_engine.insights import ZernikeDashboardInsight, InsightConfig + +study_path = Path('studies/M1_Mirror/m1_mirror_cost_reduction_V11') + +# Find the iteration with the best design +best_archive = study_path / '3_results' / 'best_design_archive' +if best_archive.exists(): + op2_files = list(best_archive.glob('*.op2')) + if op2_files: + print(f'Found best design archive: {op2_files[0].name}') + +insight = ZernikeDashboardInsight(study_path) +if insight.can_generate(): + print('Generating dashboard...') + result = insight.generate(InsightConfig()) + if result.success: + print(f'\nDashboard generated: {result.html_path}') + print(f'\nDashboard Summary:') + for key, value in result.summary.items(): + print(f' {key}: {value}') + else: + print(f'Error: {result.error}') +else: + print('Cannot generate insight - no OP2 files found') diff --git a/tests/compare_v8_zernike_methods.py b/tests/compare_v8_zernike_methods.py new file mode 100644 index 00000000..5142014a --- /dev/null +++ b/tests/compare_v8_zernike_methods.py @@ -0,0 +1,252 @@ +""" +Compare V8 Best Candidate: OPD Method vs Standard Z-Only Method + +This script extracts WFE using both methods and compares the results +to quantify the difference between using full X,Y,Z displacement (OPD) +vs just Z displacement (Standard). +""" +import sys +sys.path.insert(0, '.') + +import sqlite3 +from pathlib import Path +import json +import numpy as np + +# Find V8 best trial +db_path = Path('studies/M1_Mirror/m1_mirror_cost_reduction_V8/3_results/study.db') +print(f"V8 Database: {db_path}") +print(f"Exists: {db_path.exists()}") + +if not db_path.exists(): + print("ERROR: V8 database not found!") + sys.exit(1) + +conn = sqlite3.connect(db_path) +c = conn.cursor() + +# Get all completed trials with their objectives +c.execute(''' + SELECT t.trial_id + FROM trials t + WHERE t.state = 'COMPLETE' + ORDER BY t.trial_id +''') +completed_ids = [row[0] for row in c.fetchall()] +print(f"Completed trials: {len(completed_ids)}") + +# Build trial data and find best +trials = [] +for tid in completed_ids: + c.execute(''' + SELECT key, value_json + FROM trial_user_attributes + WHERE trial_id = ? + ''', (tid,)) + attrs = {row[0]: json.loads(row[1]) for row in c.fetchall()} + + # V8 used different objective names - check what's available + rms_40 = attrs.get('rel_filtered_rms_40_vs_20', attrs.get('filtered_rms_40_20', None)) + rms_60 = attrs.get('rel_filtered_rms_60_vs_20', attrs.get('filtered_rms_60_20', None)) + mfg_90 = attrs.get('mfg_90_optician_workload', attrs.get('optician_workload_90', None)) + ws = attrs.get('weighted_sum', None) + + if rms_40 is not None: + trials.append({ + 'id': tid, + 'rms_40': rms_40, + 'rms_60': rms_60 if rms_60 else 999, + 'mfg_90': mfg_90 if mfg_90 else 999, + 'ws': ws if ws else 999, + }) + +conn.close() + +if not trials: + # Check what keys are available + conn = sqlite3.connect(db_path) + c = conn.cursor() + c.execute('SELECT DISTINCT key FROM trial_user_attributes LIMIT 20') + keys = [row[0] for row in c.fetchall()] + print(f"\nAvailable attribute keys: {keys}") + conn.close() + print("ERROR: No trials found with expected objective names!") + sys.exit(1) + +# Calculate weighted sum if not present +for t in trials: + if t['ws'] == 999: + w40 = 100 / 4.0 + w60 = 50 / 10.0 + w90 = 20 / 20.0 + t['ws'] = w40 * t['rms_40'] + w60 * t['rms_60'] + w90 * t['mfg_90'] + +# Find best +trials.sort(key=lambda x: x['ws']) +best = trials[0] + +print(f"\n{'='*70}") +print("V8 Best Trial Summary") +print('='*70) +print(f"Best Trial: #{best['id']}") +print(f" 40-20 RMS: {best['rms_40']:.2f} nm") +print(f" 60-20 RMS: {best['rms_60']:.2f} nm") +print(f" MFG 90: {best['mfg_90']:.2f} nm") +print(f" WS: {best['ws']:.1f}") + +# Now compare both extraction methods on this trial +iter_path = Path(f'studies/M1_Mirror/m1_mirror_cost_reduction_V8/2_iterations/iter{best["id"]}') +op2_path = iter_path / 'assy_m1_assyfem1_sim1-solution_1.op2' +geo_path = iter_path / 'assy_m1_assyfem1_sim1-solution_1.dat' + +print(f"\nIteration path: {iter_path}") +print(f"OP2 exists: {op2_path.exists()}") +print(f"Geometry exists: {geo_path.exists()}") + +if not op2_path.exists(): + print("ERROR: OP2 file not found for best trial!") + sys.exit(1) + +print(f"\n{'='*70}") +print("Comparing Zernike Methods: OPD vs Standard") +print('='*70) + +# Import extractors +from optimization_engine.extractors.extract_zernike_figure import ZernikeOPDExtractor +from optimization_engine.extractors.extract_zernike import ZernikeExtractor + +# Standard method (Z-only) +print("\n1. STANDARD METHOD (Z-only displacement)") +print("-" * 50) +try: + std_extractor = ZernikeExtractor( + op2_path, + bdf_path=geo_path, + n_modes=50, + filter_orders=4 + ) + + # Extract relative WFE + std_40_20 = std_extractor.extract_relative(target_subcase='3', reference_subcase='2') + std_60_20 = std_extractor.extract_relative(target_subcase='4', reference_subcase='2') + + # MFG uses J1-J3 filter - need new extractor instance + std_extractor_mfg = ZernikeExtractor( + op2_path, + bdf_path=geo_path, + n_modes=50, + filter_orders=3 # J1-J3 for manufacturing + ) + std_90 = std_extractor_mfg.extract_subcase(subcase_label='1') + + print(f" 40-20 Relative RMS: {std_40_20['relative_filtered_rms_nm']:.2f} nm") + print(f" 60-20 Relative RMS: {std_60_20['relative_filtered_rms_nm']:.2f} nm") + print(f" 90 MFG (J1-J3): {std_90['filtered_rms_nm']:.2f} nm") + + std_results = { + '40_20': std_40_20['relative_filtered_rms_nm'], + '60_20': std_60_20['relative_filtered_rms_nm'], + '90_mfg': std_90['filtered_rms_nm'], + } +except Exception as e: + print(f" ERROR: {e}") + import traceback + traceback.print_exc() + std_results = None + +# OPD method (X,Y,Z displacement with mesh interpolation) +print("\n2. OPD METHOD (X,Y,Z displacement with mesh interpolation)") +print("-" * 50) +try: + opd_extractor = ZernikeOPDExtractor( + op2_path, + bdf_path=geo_path, + n_modes=50, + filter_orders=4 + ) + + # Extract relative WFE using OPD method + opd_40_20 = opd_extractor.extract_relative(target_subcase='3', reference_subcase='2') + opd_60_20 = opd_extractor.extract_relative(target_subcase='4', reference_subcase='2') + + # MFG uses J1-J3 filter + opd_extractor_mfg = ZernikeOPDExtractor( + op2_path, + bdf_path=geo_path, + n_modes=50, + filter_orders=3 # J1-J3 for manufacturing + ) + opd_90 = opd_extractor_mfg.extract_subcase(subcase_label='1') + + print(f" 40-20 Relative RMS: {opd_40_20['relative_filtered_rms_nm']:.2f} nm") + print(f" 60-20 Relative RMS: {opd_60_20['relative_filtered_rms_nm']:.2f} nm") + print(f" 90 MFG (J1-J3): {opd_90['filtered_rms_nm']:.2f} nm") + + # Also get lateral displacement info + print(f"\n Lateral Displacement (40° vs 20°):") + print(f" Max: {opd_40_20.get('max_lateral_displacement_um', 'N/A')} µm") + print(f" RMS: {opd_40_20.get('rms_lateral_displacement_um', 'N/A')} µm") + + opd_results = { + '40_20': opd_40_20['relative_filtered_rms_nm'], + '60_20': opd_60_20['relative_filtered_rms_nm'], + '90_mfg': opd_90['filtered_rms_nm'], + 'lateral_max': opd_40_20.get('max_lateral_displacement_um', 0), + 'lateral_rms': opd_40_20.get('rms_lateral_displacement_um', 0), + } +except Exception as e: + print(f" ERROR: {e}") + import traceback + traceback.print_exc() + opd_results = None + +# Comparison +if std_results and opd_results: + print(f"\n{'='*70}") + print("COMPARISON: OPD vs Standard Method") + print('='*70) + + print(f"\n{'Metric':<25} {'Standard':<12} {'OPD':<12} {'Delta':<12} {'Delta %':<10}") + print("-" * 70) + + for key, label in [('40_20', '40-20 RMS (nm)'), ('60_20', '60-20 RMS (nm)'), ('90_mfg', '90 MFG (nm)')]: + std_val = std_results[key] + opd_val = opd_results[key] + delta = opd_val - std_val + delta_pct = 100 * delta / std_val if std_val > 0 else 0 + + print(f"{label:<25} {std_val:<12.2f} {opd_val:<12.2f} {delta:+<12.2f} {delta_pct:+.1f}%") + + print(f"\n{'='*70}") + print("INTERPRETATION") + print('='*70) + + delta_40 = opd_results['40_20'] - std_results['40_20'] + delta_60 = opd_results['60_20'] - std_results['60_20'] + + print(f""" +The OPD method accounts for lateral (X,Y) displacement when computing WFE. + +For telescope mirrors with lateral supports: +- Gravity causes the mirror to shift laterally (X,Y) as well as sag (Z) +- The Standard method ignores this lateral shift +- The OPD method interpolates the ideal surface at deformed (x+dx, y+dy) positions + +Key observations: +- 40-20 difference: {delta_40:+.2f} nm ({100*delta_40/std_results['40_20']:+.1f}%) +- 60-20 difference: {delta_60:+.2f} nm ({100*delta_60/std_results['60_20']:+.1f}%) +- Lateral displacement: Max {opd_results['lateral_max']:.3f} µm, RMS {opd_results['lateral_rms']:.3f} µm + +Significance: +""") + + if abs(delta_40) < 0.5 and abs(delta_60) < 0.5: + print(" -> SMALL DIFFERENCE: For this design, lateral displacement is minimal.") + print(" Both methods give similar results.") + else: + print(" -> SIGNIFICANT DIFFERENCE: Lateral displacement affects WFE computation.") + print(" OPD method is more physically accurate for this geometry.") + + if opd_results['lateral_rms'] > 0.1: + print(f"\n WARNING: Lateral RMS {opd_results['lateral_rms']:.3f} µm is notable.") + print(" OPD method recommended for accurate optimization.") diff --git a/tests/test_zernike_methods_comparison.py b/tests/test_zernike_methods_comparison.py new file mode 100644 index 00000000..46cd611f --- /dev/null +++ b/tests/test_zernike_methods_comparison.py @@ -0,0 +1,208 @@ +#!/usr/bin/env python +""" +Compare all Zernike extraction methods on M1 Mirror data. + +Methods compared: +1. Standard - Z-displacement only at original (x,y) +2. Parabola OPD - Uses parabola approximation with prescription focal length +3. Figure OPD - Uses actual figure.dat geometry (most rigorous) + +Run: + python tests/test_zernike_methods_comparison.py studies/M1_Mirror/m1_mirror_cost_reduction_V9 +""" + +from pathlib import Path +import sys + +sys.path.insert(0, str(Path(__file__).parent.parent)) + +import numpy as np + + +def run_full_comparison(study_dir: Path): + """Run comparison of all three Zernike methods.""" + from optimization_engine.extractors.extract_zernike import ZernikeExtractor + from optimization_engine.extractors.extract_zernike_opd import ZernikeAnalyticExtractor + from optimization_engine.extractors.extract_zernike_figure import ZernikeOPDExtractor + + # Find OP2 file + op2_files = list(study_dir.glob('3_results/best_design_archive/**/*.op2')) + if not op2_files: + op2_files = list(study_dir.glob('2_iterations/iter1/*.op2')) + if not op2_files: + raise FileNotFoundError(f"No OP2 file found in {study_dir}") + + op2_file = op2_files[0] + + # Figure file is optional - extractor will use BDF geometry filtered to OP2 nodes + figure_file = study_dir / '1_setup' / 'model' / 'figure.dat' + use_figure_file = figure_file.exists() + + print("=" * 80) + print("ZERNIKE METHODS COMPARISON - M1 MIRROR") + print("=" * 80) + print(f"\nOP2 file: {op2_file.name}") + if use_figure_file: + print(f"Figure file: {figure_file.name}") + else: + print("Figure file: Not found (using BDF geometry filtered to OP2 nodes)") + + # Initialize extractors + print("\nInitializing extractors...") + + std_extractor = ZernikeExtractor(op2_file) + print(f" Standard: {len(std_extractor.node_geometry)} nodes in BDF") + + analytic_extractor = ZernikeAnalyticExtractor(op2_file, focal_length=1445.0, concave=True) + print(f" Analytic (Parabola f=1445mm): {len(analytic_extractor.node_geometry)} nodes") + + # OPD extractor - ALWAYS use BDF geometry filtered to OP2 nodes (RECOMMENDED) + # (figure.dat may have mismatched coordinates from different model state) + opd_extractor = ZernikeOPDExtractor( + op2_file, + figure_path=None # Force use of BDF geometry + ) + print(f" OPD (BDF geometry): {len(opd_extractor.figure_geometry)} figure nodes") + + print() + + # Get available subcases + subcases = list(std_extractor.displacements.keys()) + print(f"Available subcases: {subcases}") + + # Results table header + print("\n" + "=" * 80) + print(f"{'Subcase':<10} {'Standard':<15} {'Analytic':<15} {'OPD':<15} {'Max Lat (um)':<12}") + print("-" * 80) + + all_results = {} + + for sc in subcases: + try: + # Extract with each method + std_res = std_extractor.extract_subcase(sc) + analytic_res = analytic_extractor.extract_subcase(sc) + opd_res = opd_extractor.extract_subcase(sc) + + std_rms = std_res['filtered_rms_nm'] + analytic_rms = analytic_res['filtered_rms_nm'] + opd_rms = opd_res['filtered_rms_nm'] + max_lat = opd_res.get('max_lateral_displacement_um', 0) + + print(f"{sc:<10} {std_rms:<15.2f} {analytic_rms:<15.2f} {opd_rms:<15.2f} {max_lat:<12.3f}") + + all_results[sc] = { + 'standard': std_res, + 'analytic': analytic_res, + 'opd': opd_res + } + + except Exception as e: + print(f"{sc:<10} ERROR: {e}") + + print("-" * 80) + + # Detailed comparison for each subcase + print("\n" + "=" * 80) + print("DETAILED COMPARISON") + print("=" * 80) + + for sc, results in all_results.items(): + print(f"\n{'-' * 40}") + print(f"SUBCASE {sc}") + print(f"{'-' * 40}") + + std = results['standard'] + analytic = results['analytic'] + opd = results['opd'] + + print(f"\n{'Metric':<25} {'Standard':<15} {'Analytic':<15} {'OPD':<15}") + print("-" * 70) + + # RMS metrics + print(f"{'Filtered RMS (nm)':<25} {std['filtered_rms_nm']:<15.2f} {analytic['filtered_rms_nm']:<15.2f} {opd['filtered_rms_nm']:<15.2f}") + print(f"{'Global RMS (nm)':<25} {std['global_rms_nm']:<15.2f} {analytic['global_rms_nm']:<15.2f} {opd['global_rms_nm']:<15.2f}") + + # Aberrations + for aberr in ['defocus', 'astigmatism', 'coma', 'trefoil', 'spherical']: + key = f'{aberr}_nm' + if key in std and key in analytic and key in opd: + print(f"{aberr.capitalize():<25} {std[key]:<15.2f} {analytic[key]:<15.2f} {opd[key]:<15.2f}") + + # Node count + print(f"{'Nodes':<25} {std['n_nodes']:<15} {analytic['n_nodes']:<15} {opd['n_nodes']:<15}") + + # Lateral displacement (only for OPD methods) + print() + print(f"Lateral displacement (OPD method):") + print(f" Max: {opd.get('max_lateral_displacement_um', 0):.3f} um") + print(f" RMS: {opd.get('rms_lateral_displacement_um', 0):.3f} um") + print(f" Mean: {opd.get('mean_lateral_displacement_um', 0):.3f} um") + + # Differences + print() + diff_std_opd = opd['filtered_rms_nm'] - std['filtered_rms_nm'] + diff_analytic_opd = opd['filtered_rms_nm'] - analytic['filtered_rms_nm'] + + print(f"Difference from Standard:") + print(f" OPD: {diff_std_opd:+.2f} nm ({100*diff_std_opd/std['filtered_rms_nm']:+.1f}%)") + print(f" Analytic: {analytic['filtered_rms_nm'] - std['filtered_rms_nm']:+.2f} nm") + print() + print(f"Difference OPD vs Analytic: {diff_analytic_opd:+.2f} nm") + + # Tracking WFE analysis + if '2' in all_results and '3' in all_results and '4' in all_results: + print("\n" + "=" * 80) + print("TRACKING WFE ANALYSIS (elevation changes)") + print("=" * 80) + + for method_name, method_key in [('Standard', 'standard'), ('Analytic', 'analytic'), ('OPD', 'opd')]: + print(f"\n{method_name}:") + + z20 = np.array(all_results['2'][method_key].get('coefficients', [0]*36)) + z40 = np.array(all_results['3'][method_key].get('coefficients', [0]*36)) + z60 = np.array(all_results['4'][method_key].get('coefficients', [0]*36)) + + if len(z20) > 4: + # Differential (J5+) + diff_40_20 = z40 - z20 + diff_60_20 = z60 - z20 + + rms_40_20 = np.sqrt(np.sum(diff_40_20[4:]**2)) + rms_60_20 = np.sqrt(np.sum(diff_60_20[4:]**2)) + + print(f" 40-20 deg tracking: {rms_40_20:.2f} nm RMS (J5+ filtered)") + print(f" 60-20 deg tracking: {rms_60_20:.2f} nm RMS (J5+ filtered)") + else: + print(f" (coefficients not available)") + + print("\n" + "=" * 80) + print("SUMMARY") + print("=" * 80) + print(""" +Key findings: +- Standard method: Uses Z-displacement only at original (x,y) - fast but ignores lateral shift +- Analytic method: Accounts for lateral shift using parabola formula (requires focal length) +- OPD method: Uses actual mesh geometry - MOST RIGOROUS, no shape assumption + +Recommendation: Use OPD method (ZernikeOPDExtractor) for all mirror optimization. +The Analytic method is useful for comparison against theoretical parabola. +""") + + +def main(): + import argparse + parser = argparse.ArgumentParser(description='Compare Zernike extraction methods') + parser.add_argument('study_dir', type=str, help='Path to study directory') + args = parser.parse_args() + + study_dir = Path(args.study_dir).resolve() + if not study_dir.exists(): + print(f"Study directory not found: {study_dir}") + sys.exit(1) + + run_full_comparison(study_dir) + + +if __name__ == '__main__': + main()