ms-11: طرح را بهبود دهید، سپس بهبود را به چالش بکشید

from pathlib import Path
import sys, json, hashlib, platform
from dataclasses import asdict, replace
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
ROOT = Path.cwd()
if not (ROOT/'simulations').is_dir(): ROOT = ROOT.parent
sys.path.insert(0, str(ROOT/'simulations'))
from design_studies import (Design, order_metrics, spectral_sweep, bandwidth,
                           local_prediction, optimize_widths, grcwa_check)
from unit_cell_te import Cell, solve_cell, sample_ey
plt.rcParams.update({'figure.dpi': 125, 'font.size': 10, 'axes.spines.top': False,
                     'axes.spines.right': False})
OUT = ROOT/'data/design-studies'; OUT.mkdir(parents=True, exist_ok=True)
def metadata(design):
    import scipy
    return dict(parameters=asdict(design), units='um, degrees; incident-normalized power',
      conventions='exp(-i omega t), TE Ey, r at z=0, t at z=height',
      sources={name:hashlib.sha256((ROOT/'simulations'/name).read_text(encoding='utf-8').encode()).hexdigest()
               for name in ['design_studies.py','unit_cell_te.py','metasurface_synthesis.py']},
      environment={'python':platform.python_version(),'numpy':np.__version__,'scipy':scipy.__version__})
def geometry(ax, design, title):
    L=len(design.widths_um)*design.pitch_um; h=design.height_um
    ax.add_patch(Rectangle((0,h),L,.32,fc='#c3ddeb'))
    for j,w in enumerate(design.widths_um):
        ax.add_patch(Rectangle(((j+.5)*design.pitch_um-w/2,0),w,h,fc='#167b9b',ec='white'))
        ax.text((j+.5)*design.pitch_um,-.07,str(j+1),ha='center',fontsize=9)
    ax.axvline(0,c='#9866aa',ls='--');ax.axvline(L,c='#9866aa',ls='--')
    ax.set(xlim=(-.1,L+.1),ylim=(h+.32,-.4),xlabel='x (um)',ylabel='z (um)',title=title)
    ax.text(.05,.04,'TE: E along y; k along +z',transform=ax.transAxes,fontsize=8)
    ax.text(.05,.13,'Repeated in x; invariant in y',transform=ax.transAxes,fontsize=8)

baseline=Design()
result=optimize_widths(baseline,starts=2,maxfev=1000,half_order=32,bounds=(.08,.57),seed=17)
selected=replace(baseline,widths_um=tuple(result['best']['widths_um']))
# Geometry is frozen here; subsequent checks are never fed back into this search.
checks=[]
for M in [32,64,96,144,192]:
    checks.append(dict(half_order=M,baseline=order_metrics(baseline.solve(M)),selected=order_metrics(selected.solve(M))))
assert checks[-1]['selected']['T']>checks[-1]['baseline']['T']
assert abs(checks[-1]['selected']['T']-checks[-2]['selected']['T'])<1e-4
print('Optimizer runs:',json.dumps(result['runs'],indent=2))
print('Baseline / selected, checked M=192:',checks[-1]['baseline']['T'],checks[-1]['selected']['T'])
print('Selected widths (um):',selected.widths_um)

fig,ax=plt.subplots(1,3,figsize=(13,3.8),layout='constrained')
geometry(ax[0],baseline,'Baseline geometry');geometry(ax[1],selected,'Selected geometry, same height/pitch')
history=result['history'];ax[2].plot([r['evaluation'] for r in history],[r['T'] for r in history],'.',ms=2,alpha=.35,label='Queried candidates')
ax[2].plot([r['evaluation'] for r in history],[r['best_T'] for r in history],label='Best seen',c='#128d85')
ax[2].axhline(checks[0]['baseline']['T'],ls='--',c='#888',label='Baseline at search resolution')
ax[2].set(xlabel='Full-wave objective evaluations',ylabel='T(+1) at M=32',title='Retain the best candidate');ax[2].legend(fontsize=7)
plt.show()

# Withheld operating points: diagnostics, not a broadband training objective.
held_out=[]
for wl in [1.45,1.50,1.60,1.65]:
    held_out.append(dict(wavelength_um=wl,baseline=order_metrics(replace(baseline,wavelength_um=wl).solve(144)),
                        selected=order_metrics(replace(selected,wavelength_um=wl).solve(144))))
# Deterministic fabrication scenarios, not a probability distribution or yield.
perturbations=[]
for j in [-1,0,1,2,3,4,5]:
    for delta in [-.005,.005]:
        w=np.array(selected.widths_um)
        if j<0: w+=delta
        else: w[j]+=delta
        perturbations.append(dict(ridge='all' if j<0 else j+1,delta_um=delta,
            T=order_metrics(replace(selected,widths_um=tuple(w)).solve(144))['T']))
# Optional separate implementation; no commercial solver is needed.
external=[]
try:
    import grcwa
except ImportError:
    external_status='not run: install requirements-validation.txt and rerun this notebook'
else:
    ours=selected.solve(48)
    for pixels in [4096,8192]:
        other=grcwa_check(selected,48,pixels);idx=other['orders']+48
        error=max(abs(ours['R'][idx]-other['R']).max(),abs(ours['T'][idx]-other['T']).max())
        external.append(dict(half_order=48,pixels=pixels,max_order_power_error=float(error)))
    assert external[-1]['max_order_power_error']<1e-4
    external_status='passed: separate RCWA implementation; not FEM/FDTD'
print('Independent implementation:',external_status,external)
print('Held-out wavelength checks:',held_out)
print('Worst of 14 specified perturbations:',min(r['T'] for r in perturbations))

fig,ax=plt.subplots(1,3,figsize=(13,3.5),layout='constrained')
for name,color in [('baseline','#777'),('selected','#128d85')]:
    ax[0].plot([c['half_order'] for c in checks],[c[name]['T'] for c in checks],'o-',label=name,c=color)
    ax[1].plot([c['wavelength_um'] for c in held_out],[c[name]['T'] for c in held_out],'o--',label=name,c=color)
ax[0].set(xlabel='Retained half-order M',ylabel='T(+1)',title='Refine after freezing geometry');ax[0].legend(fontsize=8)
ax[1].set(xlabel='Withheld wavelength (um)',ylabel='T(+1)',title='Off-design checks (sampled only)');ax[1].legend(fontsize=8)
ax[2].bar(np.arange(14),[r['T'] for r in perturbations],color='#128d85')
ax[2].axhline(checks[-1]['selected']['T'],c='#777',ls='--',label='Nominal')
ax[2].set(xticks=np.arange(14),xticklabels=[str(r['ridge'])+('−' if r['delta_um']<0 else '+') for r in perturbations],
    xlabel='Perturbed ridge; minus/plus 5 nm',ylabel='T(+1)',title='Specified width-error scenarios');ax[2].tick_params(axis='x',rotation=60);ax[2].legend(fontsize=8)
plt.show()
report=dict(metadata=metadata(baseline),optimization=result,selected_parameters=asdict(selected),
    refinement=checks,held_out=held_out,perturbations=perturbations,
    independent_implementation=dict(status=external_status,checks=external))
(OUT/'inverse-design.json').write_text(json.dumps(report,indent=2)+'\n')