{
"cells": [
{
"cell_type": "markdown",
"id": "fed083e4-9be1-480c-82e9-4bfd7dbccf1f",
"metadata": {},
"source": [
"# Adjoint Plugin: 2 Checking Gradients\n",
"\n",
"In this notebook, we will show how to use the adjoint plugin for `DiffractionMonitor` outputs and also check the gradient values against gradients obtained using transfer matrix method (TMM) to validate their accuracy for a multilayer slab problem."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "51963da1",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import jax.numpy as jnp\n",
"import jax\n",
"import tmm\n",
"import matplotlib.pyplot as plt\n",
"from typing import Tuple, List\n",
"\n",
"import tidy3d as td\n",
"from tidy3d.web import run as run_sim\n",
"from tidy3d.plugins.adjoint import (\n",
" JaxSimulation,\n",
" JaxBox,\n",
" JaxMedium,\n",
" JaxStructure,\n",
" JaxSimulationData,\n",
")\n",
"from tidy3d.plugins.adjoint.web import run as run_adjoint\n"
]
},
{
"cell_type": "markdown",
"id": "b33d0892-62d3-4329-810f-f29ce10056b4",
"metadata": {},
"source": [
"First, we define some global parameters describing the transmission through a multilayer slab with some spacing between each slab.\n",
"\n",
"The layout is diagrammed below.\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "ccbd2c6a-419d-4570-8f91-f57ad5bad813",
"metadata": {},
"outputs": [],
"source": [
"# frequency we want to simulate at\n",
"freq0 = 2.0e14\n",
"k0 = 2 * np.pi * freq0 / td.C_0\n",
"freqs = [freq0]\n",
"wavelength = td.C_0 / freq0\n",
"\n",
"# background permittivity\n",
"bck_eps = 1.3**2\n",
"\n",
"# space between each slab\n",
"spc = 0.0\n",
"\n",
"# slab permittivities and thicknesses\n",
"slab_eps0 = [2.0**2, 1.8**2, 1.5**2, 1.9**2]\n",
"slab_ds0 = [0.5, 0.25, 0.5, 0.5]\n",
"\n",
"# incidence angle\n",
"theta = 0 * np.pi / 8\n",
"\n",
"# resolution\n",
"dl = 0.01\n"
]
},
{
"cell_type": "markdown",
"id": "48dd37b0-ba97-45cd-a861-53ea6ae3fa74",
"metadata": {},
"source": [
"## Transfer Matrix Method (Ground Truth)\n",
"\n",
"Next we use the `tmm` package to write a function to return the transmission `T` of `p` polarized light given a set of slab permittivities and thicknesses. We'll also write a function to compute the numerical gradient using TMM and will take these to be our \"ground truths\" when evaluating the accuracy of our values obtained through FDTD and the adjoint plugin.\n",
"\n",
"### Transmission Calculation with TMM\n",
"\n",
"First, we write a function to compute transmission."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "cc4f6e91-6a53-476e-b4c6-e76e50246426",
"metadata": {},
"outputs": [],
"source": [
"def compute_T_tmm(slab_eps=slab_eps0, slab_ds=slab_ds0) -> float:\n",
" \"\"\"Get transmission as a function of slab permittivities and thicknesses.\"\"\"\n",
"\n",
" # construct lists of permittivities and thicknesses including spaces between\n",
" new_slab_eps = []\n",
" new_slab_ds = []\n",
" for eps, d in zip(slab_eps, slab_ds):\n",
" new_slab_eps.append(eps)\n",
" new_slab_eps.append(bck_eps)\n",
" new_slab_ds.append(d)\n",
" new_slab_ds.append(spc)\n",
" slab_eps = new_slab_eps[:-1]\n",
" slab_ds = new_slab_ds[:-1]\n",
"\n",
" # add the input and output spaces to the lists\n",
" eps_list = [bck_eps] + slab_eps + [bck_eps]\n",
" n_list = np.sqrt(eps_list)\n",
" d_list = [np.inf] + slab_ds + [np.inf]\n",
"\n",
" # compute transmission with TMM\n",
" return tmm.coh_tmm(\"p\", n_list, d_list, theta, wavelength)[\"T\"]\n"
]
},
{
"cell_type": "markdown",
"id": "c6c0da6b-2d2f-4be5-80e6-c84afd8d1a23",
"metadata": {},
"source": [
"We run this function with our starting parameters and see that we get a transmission of about 98% for the set of input parameters."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "c1743e54-d7f3-468a-8ac4-d42d4c2bc8dc",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"T (tmm) = 0.901\n"
]
}
],
"source": [
"T_tmm = compute_T_tmm(slab_eps=slab_eps0, slab_ds=slab_ds0)\n",
"print(f\"T (tmm) = {T_tmm:.3f}\")\n"
]
},
{
"cell_type": "markdown",
"id": "452db467-f486-468d-8ec9-b9e669755c92",
"metadata": {},
"source": [
"### Numerical Gradient with TMM\n",
"\n",
"Next, we will use our `compute_T_tmm()` function to compute the \"numerical\" gradient to use as comparison against our adjoint results with FDTD.\n",
"\n",
"The derivative of a function $f(x)$ w.r.t. $x$ can be approximated using finite differences as\n",
"\n",
"$$ \\frac{df}{dx}(x) \\approx \\frac{f(x+\\Delta) - f(x-\\Delta)}{2\\Delta}$$\n",
"\n",
"with a small step $\\Delta$.\n",
"\n",
"To compute the gradient of our transmission with respect to each of the slab thicknesses and permittivities, we need to repeat this step for each of the values. Luckily, since TMM is very fast, we can compute these quantities quite quickly compared to if we were using FDTD.\n",
"\n",
"> Important note: We assume in our TMM numerical gradient that when the slabs are touching (`spc=0`) and a slab thickness is modified, that the thicknesses of the neighboring slabs adjust to accomidate this change. For example, if slab `i` increases by `dt`, slab `i-1` and `i+1` each decrease by `dt/2`. We also account for this in our FDTD set up by keeping the centers of all boxes constant and not tracking the gradient through these quantities. The reason this is required is that `tidy3d.plugins.adjoint` does not recognize the space between touching `JaxBox` objects as a single interface and will instead \"double count\" the gradient contribution of the interface if they are placed right next to each other. One must therefore be careful about overlapping or touching two `JaxBox` or other geometries when computing gradients.\n",
"\n",
"Here we write the function to return the numerical gradient."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "15d5bae1-6fc7-41c8-8772-2e80505be17d",
"metadata": {},
"outputs": [],
"source": [
"def compute_grad_tmm(\n",
" slab_eps=slab_eps0, slab_ds=slab_ds0\n",
") -> Tuple[List[float], List[float]]:\n",
" \"\"\"Compute numerical gradient of transmission w.r.t. each of the slab permittivities and thicknesses using TMM.\"\"\"\n",
"\n",
" delta = 1e-4\n",
"\n",
" # set up containers to store gradient and perturbed arguments\n",
" num_slabs = len(slab_eps)\n",
" grad_tmm = np.zeros((2, num_slabs), dtype=float)\n",
" args = np.stack((slab_eps, slab_ds), axis=0)\n",
"\n",
" # loop through slab index and argument index (eps, d)\n",
" for arg_index in range(2):\n",
" for slab_index in range(num_slabs):\n",
" grad = 0.0\n",
"\n",
" # perturb the argument by delta in each + and - direction\n",
" for pm in (-1, +1):\n",
" args_num = args.copy()\n",
" args_num[arg_index][slab_index] += delta * pm\n",
"\n",
" # NEW: for slab thickness gradient, need to modify neighboring slabs too\n",
" if arg_index == 1 and spc == 0:\n",
" if slab_index > 0:\n",
" args_num[arg_index][slab_index - 1] -= delta * pm / 2\n",
" if slab_index < num_slabs - 1:\n",
" args_num[arg_index][slab_index + 1] -= delta * pm / 2\n",
"\n",
" # compute argument perturbed T and add to finite difference gradient contribution\n",
" T_tmm = compute_T_tmm(slab_eps=args_num[0], slab_ds=args_num[1])\n",
" grad += pm * T_tmm / 2 / delta\n",
"\n",
" grad_tmm[arg_index][slab_index] = grad\n",
" grad_eps, grad_ds = grad_tmm\n",
" return grad_eps, grad_ds\n"
]
},
{
"cell_type": "markdown",
"id": "d3c362ca-fb31-4ead-8e7c-fcc24d0faecc",
"metadata": {},
"source": [
"Let's run this function and observe the gradients. These will be saved later to compare against our adjoint plugin results."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "b5d1cc93-6d6e-43a7-8e5f-f1829770f3d5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"gradient w.r.t. eps (tmm) = [-0.15463022 0.0376046 -0.0850184 -0.15883104]\n",
"gradient w.r.t. ds (tmm) = [-0.86661161 -0.12292531 0.58010922 -1.05537497]\n"
]
}
],
"source": [
"grad_eps_tmm, grad_ds_tmm = compute_grad_tmm()\n",
"print(f\"gradient w.r.t. eps (tmm) = {grad_eps_tmm}\")\n",
"print(f\"gradient w.r.t. ds (tmm) = {grad_ds_tmm}\")\n"
]
},
{
"cell_type": "markdown",
"id": "50ca68ae-ca5a-49aa-9969-ac6d43324ea9",
"metadata": {},
"source": [
"## FDTD (Using adjoint plugin)\n",
"\n",
"Next, we will implement the same two functions using Tidy3D's adjoint plugin.\n",
"\n",
"### Transmission Calculation with FDTD\n",
"\n",
"We first write a function to compute the transmission of a multilayer slab using Tidy3D.\n",
"\n",
"As discussed in the previous adjoint tutorial notebook, we need to use `jax`-compatible components from the tidy3d subclass for any structures that may depend on the parameters. In this case, this means that the slabs must be `JaxStructures` containing `JaxBox` and `JaxMedium` and must be added to `JaxSimulation.input_structures`.\n",
"\n",
"We use a `DiffractionMonitor` to measure our transmission amplitudes. As the data corresponding to this monitor will be used in the differentiable function return value, we must add it to `JaxSimulation.output_monitors`.\n",
"\n",
"Below, we break up the transmission calculation into a few functions to make it easier to read and re-use later."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "4d12f133-a8d8-4e59-938f-fca4d2e2fa0b",
"metadata": {},
"outputs": [],
"source": [
"def make_sim(slab_eps=slab_eps0, slab_ds=slab_ds0) -> JaxSimulation:\n",
" \"\"\"Create a JaxSimulation given the slab permittivities and thicknesses.\"\"\"\n",
"\n",
" # frequency setup\n",
" wavelength = td.C_0 / freq0\n",
" fwidth = freq0 / 10.0\n",
" freqs = [freq0]\n",
"\n",
" # geometry setup\n",
" bck_medium = td.Medium(permittivity=bck_eps)\n",
"\n",
" space_above = 2\n",
" space_below = 2\n",
"\n",
" length_x = 0.1\n",
" length_y = 0.1\n",
" length_z = space_below + sum(slab_ds0) + space_above + (len(slab_ds0) - 1) * spc\n",
" sim_size = (length_x, length_y, length_z)\n",
"\n",
" # make structures\n",
" slabs = []\n",
" z_start = -length_z / 2 + space_below\n",
" for (d, eps) in zip(slab_ds, slab_eps):\n",
"\n",
" # dont track the gradient through the center of each slab\n",
" # as tidy3d doesn't have enough information to properly process the interface between touching JaxBox objects\n",
" z_center = jax.lax.stop_gradient(z_start + d / 2)\n",
" slab = JaxStructure(\n",
" geometry=JaxBox(center=[0, 0, z_center], size=[td.inf, td.inf, d]),\n",
" medium=JaxMedium(permittivity=eps),\n",
" )\n",
" slabs.append(slab)\n",
" z_start += d + spc\n",
"\n",
" # source setup\n",
" gaussian = td.GaussianPulse(freq0=freq0, fwidth=fwidth)\n",
" src_z = -length_z / 2 + space_below / 2.0\n",
"\n",
" source = td.PlaneWave(\n",
" center=(0, 0, src_z),\n",
" size=(td.inf, td.inf, 0),\n",
" source_time=gaussian,\n",
" direction=\"+\",\n",
" angle_theta=theta,\n",
" angle_phi=0,\n",
" pol_angle=0,\n",
" )\n",
"\n",
" # boundaries\n",
" boundary_x = td.Boundary.bloch_from_source(\n",
" source=source, domain_size=sim_size[0], axis=0, medium=bck_medium\n",
" )\n",
" boundary_y = td.Boundary.bloch_from_source(\n",
" source=source, domain_size=sim_size[1], axis=1, medium=bck_medium\n",
" )\n",
" boundary_spec = td.BoundarySpec(\n",
" x=boundary_x, y=boundary_y, z=td.Boundary.pml(num_layers=40)\n",
" )\n",
"\n",
" # monitors\n",
" mnt_z = length_z / 2 - space_above / 2.0\n",
" monitor_1 = td.DiffractionMonitor(\n",
" center=[0.0, 0.0, mnt_z],\n",
" size=[td.inf, td.inf, 0],\n",
" freqs=freqs,\n",
" name=\"diffraction\",\n",
" normal_dir=\"+\",\n",
" )\n",
"\n",
" # make simulation\n",
" return JaxSimulation(\n",
" size=sim_size,\n",
" grid_spec=td.GridSpec.auto(min_steps_per_wvl=100),\n",
" input_structures=slabs,\n",
" sources=[source],\n",
" output_monitors=[monitor_1],\n",
" run_time=10 / fwidth,\n",
" boundary_spec=boundary_spec,\n",
" medium=bck_medium,\n",
" subpixel=True,\n",
" shutoff=1e-8,\n",
" )\n"
]
},
{
"cell_type": "markdown",
"id": "9addea1f-4dd8-4e5e-a058-b94638fb2d70",
"metadata": {},
"source": [
"Let's generate a simulation and plot it to make sure it looks reasonable."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "bc0dc456-e742-4299-ae5e-22daa49e3997",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
[13:52:21] WARNING: 'JaxSimulation.input_structures' elements 0 and 1 are overlapping or touching. log.py:50\n",
"Geometric gradients for overlapping structures may contain errors. \n",
"
\n"
],
"text/plain": [
"\u001b[2;36m[13:52:21]\u001b[0m\u001b[2;36m \u001b[0mWARNING: \u001b[32m'JaxSimulation.input_structures'\u001b[0m elements \u001b[1;36m0\u001b[0m and \u001b[1;36m1\u001b[0m are overlapping or touching. \u001b]8;id=893464;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/log.py\u001b\\\u001b[2mlog.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=396661;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/log.py#50\u001b\\\u001b[2m50\u001b[0m\u001b]8;;\u001b\\\n",
"\u001b[2;36m \u001b[0mGeometric gradients for overlapping structures may contain errors. \u001b[2m \u001b[0m\n"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"
WARNING: 'JaxSimulation.input_structures' elements 1 and 2 are overlapping or touching. log.py:50\n",
"Geometric gradients for overlapping structures may contain errors. \n",
"
\n"
],
"text/plain": [
"\u001b[2;36m \u001b[0m\u001b[2;36m \u001b[0mWARNING: \u001b[32m'JaxSimulation.input_structures'\u001b[0m elements \u001b[1;36m1\u001b[0m and \u001b[1;36m2\u001b[0m are overlapping or touching. \u001b]8;id=976137;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/log.py\u001b\\\u001b[2mlog.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=469083;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/log.py#50\u001b\\\u001b[2m50\u001b[0m\u001b]8;;\u001b\\\n",
"\u001b[2;36m \u001b[0mGeometric gradients for overlapping structures may contain errors. \u001b[2m \u001b[0m\n"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"
WARNING: 'JaxSimulation.input_structures' elements 2 and 3 are overlapping or touching. log.py:50\n",
"Geometric gradients for overlapping structures may contain errors. \n",
"
\n"
],
"text/plain": [
"\u001b[2;36m \u001b[0m\u001b[2;36m \u001b[0mWARNING: \u001b[32m'JaxSimulation.input_structures'\u001b[0m elements \u001b[1;36m2\u001b[0m and \u001b[1;36m3\u001b[0m are overlapping or touching. \u001b]8;id=380366;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/log.py\u001b\\\u001b[2mlog.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=91250;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/log.py#50\u001b\\\u001b[2m50\u001b[0m\u001b]8;;\u001b\\\n",
"\u001b[2;36m \u001b[0mGeometric gradients for overlapping structures may contain errors. \u001b[2m \u001b[0m\n"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAANcAAANXCAYAAACmJcrwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAtiklEQVR4nO3deXiTdbr/8U/SNindwlYoS1kVioCCbIIryADqKJxRR+eoiHpQAfW4DnChoMxgx3HDFZdROQdlDqIH8egoshUcEUUGZ3ABBpRhkx3aUksLyff3h79mCN2d3k2C79d15ZI++Sa5m/ru8yRpWo9zzglAnfNGewDgREVcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcUZKXlyePx6O8vLxojwIjxGXs2Wef1cyZM6M9xo8ye/ZsTZ8+PdpjRHjppZfUpUsXJScn6+STT9ZTTz1V48uWlJRo/PjxatmypRo0aKB+/fpp4cKFFa5dsWKFzjrrLKWkpCgrK0u33XabDh06VLthHUx17drVnXvuueW2B4NBV1xc7ILBYP0PVUMXXXSRa9u2bbTHCHvuueecJHfppZe6F154wV1zzTVOkvvd735Xo8tfeeWVLjEx0d19993u+eefd/3793eJiYnuww8/jFi3Zs0al5yc7Hr27OlmzJjhJk2a5Px+vxs2bFit5q2XuA4dOlQfNxOTKosrHsRSXN9//71r0qSJu+iiiyK2X3XVVS41NdXt37+/yst/8sknTpJ7+OGHw9uKi4tdx44dXf/+/SPWXnDBBa5FixYuPz8/vO3FF190ktyCBQtqPHOt49q2bZu7/vrrXYsWLZzP53Pt2rVzN998syspKXHOOffKK684SS4vL8+NGTPGZWZmuoYNG4Yv/8wzz7hTTjnF+Xw+16JFCzd27Fh34MCBiNvYsGGD+8UvfuGaN2/u/H6/a9WqlbviiivcwYMHw2s++OADd+aZZ7pAIOBSU1Ndp06d3MSJE6udvyaXO3z4sJs8ebLr2LGj8/l8rnXr1u6ee+5xhw8fLnd9s2bNcn369HENGjRwDRs2dGeffXb4C9C2bVsnKeJUFtrSpUudJLd06dKI63v99dfd6aef7pKTk12TJk3cVVdd5bZt2xax5tprr3Wpqalu27Ztbvjw4S41NdU1bdrU3XXXXe7o0aPV3gdvvfWWu/DCC8Nfww4dOripU6dGXPbcc88tN3tVoZ1zzjnu1FNPrfC8Tp06uSFDhlQ7V1XeffddJ8m9++67EdtXrFjhJLlZs2ZVefl77rnHJSQkRATjnHMPPvigk+S2bNninHMuPz/fJSYmunvuuSdiXUlJiUtLS3M33HBDjWdOrM0h5I4dO9S3b18dPHhQN954o3JycrR9+3a98cYb+v777+Xz+cJrx44dq8zMTE2ePFlFRUWSpPvvv18PPPCABg8erDFjxmj9+vWaMWOGVq1apY8++khJSUkqLS3V0KFDVVJSoltvvVVZWVnavn273nnnHR08eFCBQEBffvmlfv7zn+vUU0/V1KlT5ff7tXHjRn300UdVzl+Ty4VCIV1yySX685//rBtvvFFdunTR2rVr9fjjj2vDhg166623wmsfeOAB3X///RowYICmTp0qn8+nTz75REuWLNGQIUM0ffp03XrrrUpLS9OkSZMkSc2bN690vpkzZ+q6665Tnz59lJubq127dumJJ57QRx99pDVr1qhhw4bhtcFgUEOHDlW/fv30yCOPaNGiRXr00UfVsWNHjRkzpsr7YebMmUpLS9Odd96ptLQ0LVmyRJMnT1ZBQYEefvhhSdKkSZOUn5+vbdu26fHHH5ckpaWlVXqd11xzjUaPHq0vvvhC3bp1C29ftWqVNmzYoHvvvTe87cCBAwoGg1XOKEkpKSlKSUmRJK1Zs0aS1Lt374g1vXr1ktfr1Zo1a3T11VdXel1r1qxRp06dlJGREbG9b9++kqTPP/9c2dnZWrt2rY4ePVrudnw+n3r06BGeo0ZqnKFzbuTIkc7r9bpVq1aVOy8UCjnn/rnnOuussyK+E+7evdv5fD43ZMiQiMcZTz/9tJPkXn75ZefcD8e7ktzcuXMrnePxxx93ktyePXtqM36NLjdr1izn9XrLHYeXHe9/9NFHzjnn/v73vzuv1+v+7d/+rdzjprL7wrnKDwuP33OVlpa6Zs2auW7durni4uLwunfeecdJcpMnTw5vu/baa50kN3Xq1Ijr7Nmzp+vVq1fVd4L74RDreDfddJNLSUmJ2DvX5rDw4MGDLjk52Y0fPz5i+2233eZSU1MjHhpUtEev6DRlypTwZcaNG+cSEhIqvO3MzEx35ZVXVjlf165d3aBBg8pt//LLL50k99xzzznnnJs7d66T5JYvX15u7eWXX+6ysrKqvJ1j1fjZwlAopLfeeksXX3xxuaolyePxRHw8evRoJSQkhD9etGiRSktLdfvtt8vr9Uasy8jI0LvvvitJCgQCkqQFCxbo+++/r3CWsu/g8+fPVygUqumnUKPLzZ07V126dFFOTo727t0bPg0aNEiStHTpUknSW2+9pVAopMmTJ0d8PlL5+6ImPvvsM+3evVtjx45VcnJyePtFF12knJyc8P1zrJtvvjni47PPPlvffPNNtbfVoEGD8L8LCwu1d+9enX322fr++++1bt26Ws8u/fB1Gz58uP74xz/K/f/33waDQc2ZM0cjRoxQampqeO1rr72mhQsXVnsaOXJk+DLFxcURR0bHSk5OVnFxcZXzFRcXy+/3V3jZsvOP/W9la6u7nWPV+LBwz549KigoiNjlV6V9+/YRH//jH/+QJHXu3Dliu8/nU4cOHcLnt2/fXnfeeacee+wxvfbaazr77LN1ySWX6Oqrrw6Hd8UVV+gPf/iD/uM//kMTJkzQ+eefr1/84he67LLLyv2PfqyaXO7vf/+7vv76a2VmZlZ4Hbt375Ykbdq0SV6vV6ecckqN7o/qVHb/SFJOTo7+/Oc/R2xLTk4uN2OjRo104MCBam/ryy+/1L333qslS5aooKAg4rz8/Pzajh42cuRIzZkzRx9++KHOOeccLVq0SLt27dI111wTse7MM8+s9XU3aNBApaWlFZ53+PDhiG8YlV2+pKSkwsuWnX/sfytbW93tHKtWj7lqozZDHO/RRx/VqFGjNH/+fH3wwQe67bbblJubq5UrV6p169Zq0KCBli9frqVLl+rdd9/V+++/rzlz5mjQoEH64IMPIvaYx89U3eVCoZC6d++uxx57rMLryM7O/tGfV12q7HOszsGDB3XuuecqIyNDU6dOVceOHZWcnKy//OUvGj9+fK2OBI43dOhQNW/eXK+++qrOOeccvfrqq8rKytLgwYMj1u3Zs6dGj7nS0tLCj/NatGihYDCo3bt3q1mzZuE1paWl2rdvn1q2bFnldbVo0ULbt28vt/27776TpPDlW7RoEbH9+LXV3U6Emh4/BoNBl5GR4YYPH17lurLHXMc/Lps9e7aT5P70pz9FbC8pKXGBQMBdeumllV7nRx995CS5SZMmVbpm2rRpTpJbuHBh9Z9MFZe78MILXatWrSIeN1Xk4YcfdpLcmjVrqlzXrVu3Gj3mKnvW69lnny23tkuXLhGPpcqeLTzelClTXHVf0nnz5jlJbtmyZRHbX3jhhXLPXv785z+v9VPxd9xxh2vUqJHbv3+/S0tLc3fccUe5NT/mMVfZY8/jny0s+3/jv//7v6uc6+67767w2cKyr3/Zs4UHDx6s8tnC66+/vsb3RY0fc3m9Xo0YMUL/93//p88++6yiSKu8/ODBg+Xz+fTkk09GrH3ppZeUn5+viy66SJJUUFCgo0ePRly2e/fu8nq94V31/v37y11/jx49JFW8Oy9Tk8v98pe/1Pbt2/Xiiy+WW1tcXBx+5nPEiBHyer2aOnVque/2x35+qampOnjwYKUzlendu7eaNWum5557LuJzeO+99/T111+H759/Vdke79gZS0tL9eyzz5Zbm5qaWuvDxGuuuUYHDhzQTTfdpEOHDlX4DN6Pecw1aNAgNW7cWDNmzIi4rhkzZiglJSXi/tm7d6/WrVsX8Zj9sssuUzAY1AsvvBDeVlJSoldeeUX9+vULH5EEAgENHjxYr776qgoLC8NrZ82apUOHDunyyy+v+Z1R4wzdD69xZWVluZSUFHf77be7559/3t1///2ua9eu4deqKttzOffP76xDhgxxTz/9tLv11ltdQkKC69OnjystLXXO/fCdtVWrVu722293zz77rHvyySddnz59XFJSkvv444+dc87953/+p+vZs6e799573YsvvuimTZvmWrVq5Vq3bh3xWtjxanK5YDDoLrzwQufxeNyVV17pnnrqKTd9+nR38803u8aNG0d8Xvfdd5+T5AYMGOAeeeQR99RTT7mRI0e6CRMmhNeMHTvWeTwe95vf/Mb98Y9/dIsXL3bOVfw6V9l9169fPzd9+nQ3ceJEl5KS4tq1axfxWuC/sufau3eva9SokWvbtq179NFH3WOPPeZ69uzpTjvttHLz/P73v3eS3B133OFmz57t3n777Sqvu0y3bt2cJNelS5cara+pZ555xklyl112mXvxxRfdyJEjnSQ3bdq0iHVl98PxryFefvnl4b3S888/7wYMGOASExPL7cVXr17t/H5/xE9oJCcn1/q1ulq/iPyPf/zDjRw50mVmZjq/3+86dOjgxo0bV+5F5Iricu6Hp95zcnJcUlKSa968uRszZkzE/zjffPONu/76613Hjh1dcnKya9y4sRs4cKBbtGhReM3ixYvd8OHDXcuWLZ3P53MtW7Z0v/rVr9yGDRuqnL2mlystLXUPPfSQ69q1q/P7/a5Ro0auV69e7oEHHih3WPHyyy+7nj17htede+65EYemO3fudBdddJFLT0+v0YvIc+bMCV9f48aNq3wR+Xg1icu5Hw6lzjjjDNegQQPXsmVL9+tf/9otWLCg3DyHDh1y//7v/+4aNmxY7YvIxyqL8sEHH6zR+tp44YUXXOfOnZ3P53MdO3Z0jz/+eLlD+MriKi4udnfffbfLyspyfr/f9enTx73//vsV3s6HH37oBgwY4JKTk11mZqYbN26cKygoqNWsHuf4vYWoW0888YTuuOMObd68WW3atIn2OFFDXKhTzjmddtppatKkSfg1wZ8qs6fi8dNSVFSkt99+W0uXLtXatWs1f/78aI8Udey5UCc2b96s9u3bq2HDhho7dqymTZsW7ZGijrgAI7wTGTBCXIARntCIQaFQSDt27FB6enqNfsLeOafCwkK1bNmyyh9cRv0irhi0Y8eOH/UDwlu3blXr1q0NJsKPQVwxKD09XZJ02223qUuXLtWu//rrr/Xkk0+GL4fYQFwxqOxQMDU1NfwetqqUvRHxx7xJE3Y4QI9hx787APGFuGJYYWFhpe++RewjrhiWkJCgvXv3ElicIq4Ylp6erqSkJAKLU8QVwzwej5o0aUJgcYq4YpzX6yWwOEVccYDA4hNxxQkCiz/EFUcILL4QV5whsPhBXHHo+MD4SY7YRFxx6tjAjv3llYgdxBXHygL7sb83HraIK855vV7eahKjiOsEwFtNYhNxAUaI6wTAb8eLTcQV50KhEM8WxijiimOhUEj79u2r0V9pRP0jrjhVFtaRI0d4tjBGEVccOjaspk2bKjGR3zMUi4grzhwfls/ni/ZIqARxxRHCii/EFScIK/4QVxwgrPhEXDGOsOIXccUw5xxhxTHiimGFhYWEFceIK4YFg0HCimPEFcPS09MJK44RVwzjJy/iG3EBRoirjs2YMUOnnnqqMjIylJGRof79++u9996L9liIAuKqY61bt9bvfvc7rV69Wp999pkGDRqk4cOH68svv6z1dfEmyPjGQX0du/jiiyM+njZtmmbMmKGVK1eqa9eutbquwsJChUIheb18D4xHxGUoGAxq7ty5KioqUv/+/StdV1JSopKSkvDHBQUF4cvv27dPTZo0IbA4xFfMwNq1a5WWlia/36+bb75Z8+bN0ymnnFLp+tzcXAUCgfApOztb0g9PxR85ckT79u1TKBSqr/FRR4jLQOfOnfX555/rk08+0ZgxY3Tttdfqq6++qnT9xIkTlZ+fHz5t3bpV0g9PxTdt2pTA4hSHhQZ8Pp9OOukkSVKvXr20atUqPfHEE3r++ecrXO/3++X3+yu9rqZNm2rv3r0cIsYZvkr1IBQKRTymqq2ywNiDxRf2XHVs4sSJuuCCC9SmTRsVFhZq9uzZysvL04IFC/6l62UPFn+Iq47t3r1bI0eO1HfffadAIKBTTz1VCxYs0M9+9rN/+boJLL4QVx176aWXTK+/osAQm/i2F4eOfwzGT3LEJuKKU8cGxq+zjk3EFcfKAuPXWccm4opzPp+PX2cdo4jrBMCbKmMTcQFG+JYXw4LBoEpLS6tdd/To0XqYBrVFXDHs6NGjOnz4cLXrahIg6h9xxbCBAwdqwIAB1a5bsWKFnnjiiXqYCLVBXDEsIyNDjRs3rtE6xB6e0ACMEBdghLgAI8QFGCEuwAhxAUaICzBCXIAR4gKMEBdghLgAI8QFGCEuwAhxAUaICzBCXIAR4gKMEBdghLgAI8QFGCEuwAhxAUaICzBCXIAR4gKMEBdghLgAI8QFGCEuwAhxAUaICzBCXIAR4gKMEBdghLgAI8QFGCEuwAhxAUaICzBCXIAR4gKMEBdghLgAI8QFGCEuwAhxAUaICzBCXIAR4gKMEBdghLgAI8QFGCEuwAhxAUaICzBCXIAR4gKMEBdghLgAI8QFGCEuwAhxAUaICzBCXIAR4gKMEBdghLgAI8QFGCEuwAhxAUaICzBCXIAR4gKMEBdghLgAI8QFGCEuwAhxAUaICzBCXIAR4gKMEBdghLgAI8QFGCEuwAhxAUaICzBCXIAR4gKMEBdghLgAI8QFGCEuwAhxAUaIq47l5uaqT58+Sk9PV7NmzTRixAitX78+2mMhCoirji1btkzjxo3TypUrtXDhQh05ckRDhgxRUVFRtEdDPUuM9gAnmvfffz/i45kzZ6pZs2ZavXq1zjnnnChNhWggLmP5+fmSpMaNG1e6pqSkRCUlJeGPCwoKzOeCPQ4LDYVCId1+++0688wz1a1bt0rX5ebmKhAIhE/Z2dn1OCWsEJehcePG6YsvvtD//M//VLlu4sSJys/PD5+2bt1aTxPCEoeFRm655Ra98847Wr58uVq3bl3lWr/fL7/fX0+Tob4QVx1zzunWW2/VvHnzlJeXp/bt20d7JEQJcdWxcePGafbs2Zo/f77S09O1c+dOSVIgEFCDBg2iPB3qE4+56tiMGTOUn5+v8847Ty1atAif5syZE+3RUM/Yc9Ux51y0R0CMYM8FGCEuwAhxAUaICzDCExqx7EihVHqgZusQc4grlu1aKm3ZXoN1G+xnQa1xWAgYIS7ACHEBRnjMFcuaD5TaDKh+3T9WSHrUfBzUDnHFsqR0ydeoZusQczgsBIwQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEZWL58uS6++GK1bNlSHo9Hb731VrRHQhQQl4GioiKddtppeuaZZ6I9CqIoMdoDnIguuOACXXDBBdEeA1FGXDGgpKREJSUl4Y8LCgqiOA3qCoeFMSA3N1eBQCB8ys7OjvZIqAPEFQMmTpyo/Pz88Gnr1q3RHgl1gMPCGOD3++X3+6M9BuoYey7ACHsuA4cOHdLGjRvDH3/77bf6/PPP1bhxY7Vp0yaKk6E+EZeBzz77TAMHDgx/fOedd0qSrr32Ws2cOTNKU6G+EZeB8847T865aI+BKOMxF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAozwInIM+/LKiUpq0LT6dcV762Ea1BZ7LsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRhJjPYAqNxb/S/QytanVLtu57avpP9dUQ8ToTbYcwFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLiPPPPOM2rVrp+TkZPXr10+ffvpptEdCPSMuA3PmzNGdd96pKVOm6C9/+YtOO+00DR06VLt37472aKhHJ3RcgwYN0gMPPFBu+4EDBzRo0CCz233sscc0evRoXXfddTrllFP03HPPKSUlRS+//LLZbSL2nNBx5eXl6emnn9aIESNUVFQU3l5aWqply5aZ3GZpaalWr16twYMHh7d5vV4NHjxYH3/8cYWXKSkpUUFBQcQJ8e+EjkuSFi1apJ07d+qMM87Q5s2bzW9v7969CgaDat68ecT25s2ba+fOnRVeJjc3V4FAIHzKzs42nxP2Tvi4WrRooWXLlql79+7q06eP8vLyoj1SORMnTlR+fn74tHXr1miPhDqQGO0BLHk8HkmS3+/X7Nmz9dvf/lbDhg3T+PHjzW6zadOmSkhI0K5duyK279q1S1lZWRVexu/3y+/3m82E6Dih91zOuYiP7733Xr322mt69NFHzW7T5/OpV69eWrx4cXhbKBTS4sWL1b9/f7PbRew5ofdc3377rTIzMyO2XXrppcrJydFnn31mdrt33nmnrr32WvXu3Vt9+/bV9OnTVVRUpOuuu87sNhF7Tui42rZtW+H2rl27qmvXrma3e8UVV2jPnj2aPHmydu7cqR49euj9998v9yQHTmwndFzRdMstt+iWW26J9hiIohP6MRcQTcQFGCEuwAhxAUaICzBCXIAR4gKM8DpXDLsiuFO9jqRWu251cKc+qId5UDvsuQAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRvhDDDFsfk5Hrep0arXrdmyQNN9+HtQOey7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ASGK0B0Dl7rqwl84555xq1y1fnqS3H6qHgVAr7LkAI8QFGCEuwAhxAUaICzBCXIAR4gKMEBdghLgAI8RVx6ZNm6YBAwYoJSVFDRs2jPY4iCLiqmOlpaW6/PLLNWbMmGiPgijjZwvr2AMPPCBJmjlzZnQHQdQRVwwoKSlRSUlJ+OOCgoIoToO6wmFhDMjNzVUgEAifsrOzoz0S6gBx1cCECRPk8XiqPK1bt+5HX//EiROVn58fPm3durUOp0e0cFhYA3fddZdGjRpV5ZoOHTr86Ov3+/3y+/0/+vKITcRVA5mZmcrMzIz2GIgzxFXHtmzZov3792vLli0KBoP6/PPPJUknnXSS0tLSojsc6hVx1bHJkyfrv/7rv8If9+zZU5K0dOlSnXfeeVGaCtHAExp1bObMmXLOlTsR1k8PcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlyAEeICjBAXYIS4ACPEBRghLsAIcQFGiAswQlxxZNZfZ2nWX2dFewzUUGK0B0DNHD56WNM/mS5J+mXXX8qf6I/uQKgWe644kbc5TzsLd2pn4U7lbc6L9jioAeKKE/O+nqeQCynkQpq3bl60x0ENEFcc2HVol5ZuXqpUX6pSfala8u0S7S7aHe2xUA3iigPv/v1dFZYUKsOfoQx/hgpLC/XuhnejPRaqQVwxzjmnN756Qx6PR16PV16PVx55NPeruXLORXs8VIG4Ytza3Wv11Z6vFPAHwtsC/oC+2vOV1u5eG8XJUB3iinFvr39bh48eVkpSSnhbSlKKDh89rLfXvx3FyVAd4ophpcFSzVs3T/4EvzweT3i7x+ORP8GveevmqeRoSRQnRFWIK4Z9vvtz7SzcqUByoNx5geQAr3nFOOKKYcu3LlfIheRL8JU7z5fg4zWvGEdcMWzNzjVK9aVWen7Za14Hig/U41SoKeKKYUVHipThz6j0/LLXvD7e8XE9ToWaIq4Y5tEPr21Vpuw1r7wtefU3FGqMuGJYmi+t2jUBf0Cb8zfbD4NaI64Y5vV6VXK0RMFQsNx5wVBQJUdLlOhN1OHg4ShMh+rwfq4YVlRapGDJDxFlB7IjzttRuCP8ni6ft/yziYg+9lwxLPe8XN3c+2YleBPKnZfgTdDNvW/W2796W7nn5UZhOlSHPVcMO7nRyfKmVf79r0VaC/XI6qGCRgX1OBVqij0XYIS4ACPEBRghLsAIT2jEgZAL6eDhgxHbmu/zyPuXTdp2cKn2/+1v0RkMVSKuGFYULFIgJaCGDRqWeyH58lUBJeV9oJWJy7SpoDBKE6IqxBXDVh9arVOanKLf/vK35c47svxFuWChSgN+HRE/oRGLiCsOVPrDux7PD+9QPuZdyogdPKEBGCEuwAiHhTGsV1ovndnwzArP+yjpdR06fEjeYKl8hUfqeTLUBHHFsNSEVGUkVvxO5O7/caMKNm+WJJVs2CCt+bz+BkONEFecaj1oUPjfe5cvlx56KIrToCI85gKMEBdghLgAI8QFGCEuwAhxAUaIqw5t3rxZN9xwg9q3b68GDRqoY8eOmjJlikpLS6M9GqKA17nq0Lp16xQKhfT888/rpJNO0hdffKHRo0erqKhIjzzySLTHQz0jrjo0bNgwDRs2LPxxhw4dtH79es2YMYO4foKIy1h+fr4aN25c5ZqSkhKVlPzzj9gVFPCr0k4EPOYytHHjRj311FO66aabqlyXm5urQCAQPmVnZ1e5HvGBuGpgwoQJ8vz/NyZWdlq3bl3EZbZv365hw4bp8ssv1+jRo6u8/okTJyo/Pz982rp1q+Wng3rCYWEN3HXXXRo1alSVazp06BD+944dOzRw4EANGDBAL7zwQrXX7/f75ff7/9UxEWOIqwYyMzOVmZlZo7Xbt2/XwIED1atXL73yyivyejk4+Kkirjq0fft2nXfeeWrbtq0eeeQR7dmzJ3xeVlZWFCdDNBBXHVq4cKE2btyojRs3qnXr1hHnOeeiNBWihWOWOjRq1Cg55yo84aeHuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRhKjPQAqV1BQoP3799doHWIPccWwpUuXavv27dWu27BhQz1Mg9oirhiWmJio5OTkGq1D7OGrEsMSEhLk8/lqtA6xhyc0TgDOuWiPgAoQV5wLhUIqLCyM9hioAHHFsVAopH379ikYDEZ7FFSAuOJUWVhHjhxRenp6tMdBBYgrDh0bVtOmTXm2MEYRV5w5PqyaPJuI6CCuOEJY8YW44gRhxR/iqmOXXHKJ2rRpo+TkZLVo0ULXXHONduzY8S9dJ2HFJ+KqYwMHDtTrr7+u9evX680339SmTZt02WWX/ejrI6z4xdNMdeyOO+4I/7tt27aaMGGCRowYoSNHjigpKalW1+WcI6w4RlyG9u/fr9dee00DBgyoMqySkhKVlJSEPy57C0lhYaGaNm1KWHGKw0ID48ePV2pqqpo0aaItW7Zo/vz5Va7Pzc1VIBAIn7KzsyVJwWCQsOIYcdXAhAkT5PF4qjytW7cuvP6ee+7RmjVr9MEHHyghIUEjR46s8odrJ06cqPz8/PBp69atkqT09HTCimMcFtbAXXfdpVGjRlW5pkOHDuF/lx3KderUSV26dFF2drZWrlyp/v37V3hZv98vv99fbjs/eRHf+OrVQGZmpjIzM3/UZUOhkCRFPKbCTwNx1aFPPvlEq1at0llnnaVGjRpp06ZNuu+++9SxY8dK91o4cfGYqw6lpKTof//3f3X++eerc+fOuuGGG3Tqqadq2bJlFR72VYc3QcY39lx1qHv37lqyZEmdXV9hYaFCoZC8Xr4HxiO+ajEsGAxq37594cdtiC/EFcPS09N15MgRAotTxBXDEhMT1bRpUwKLU8QV43w+H4HFKeKKAwQWn4grThBY/CGuOEJg8YW44gyBxQ/iikPHB8ZPcsQm4opTxwbGr7OOTcQVx8oC49dZxybiinM+n49fZx2jiOsEwJsqYxNxAUaICzBCXCeA4uLiaI+AChBXnCssLCSuGEVccaywsFAFBQVq0KBBtEdBBYgrTpWFlZGRQVwxirji0LFh8RpX7CKuOENY8YO44ghhxRfiihOEFX+IKw4QVnwirhhHWPGLuGJYcXExYcUx4ophxcXFhBXHiCuGNWjQgLDiGHHFMH7yIr4RF2CEuAAjvD88BpX9qrSioiLl5+dXu76oqCjicogNHsdXJOZ888036tixY60vt2nTpog/fI7oYs8Vgxo3bixJ2rJliwKBQHh7QUGBsrOztXXrVmVkZIS35+fnq02bNuHLITYQVwwq+zOtgUAgIqIyGRkZFW7nz7vGFr4agBHiAowQVwzy+/2aMmWK/H7/v7Qd0cWzhYAR9lyAEeICjBAXYIS4ACPEFSX79+/XVVddpeTkZHm9XiUmJqpPnz769NNP9cwzz6hdu3ZKTk5Wv3799Omnn0qSDh8+rKFDhyohIUEej0eBQECvvfaa5s6dq5ycHCUnJ8vj8ZQ79e3bt9y2YcOGRfke+AlwiIphw4a5Nm3auKSkJDdhwgTXpk0b17FjR5eSkuJ8Pp97+eWX3ZdffulGjx7tGjZs6Hbt2uVGjBjhJLkbb7zRzZ0717Vq1cpJcgkJCe73v/+9++qrr5wkJ8n95je/cYsXL3ZDhw51qampbsiQIe67774Ln/bv3x/tu+CER1xRUBZB165d3bhx45xzzr333ntOkktKSnJnnHFGeG0wGHQtW7Z0U6ZMcR6Px51++unh877++msnyTVv3tw551woFHKSXNOmTd1NN93knHPu4MGDzuv1ut69e9fjZwjnnOOwMAo+/vhjBQIBrVu3ToMHD5YkDR48WF6vV0eOHIn4G8der1eDBw/WwoUL5ZzTZZddFj4vJydHHo8n/FdOvv32W0nS999/rz/84Q/q27ev3nzzTWVmZupvf/ubmjVrps6dO2vMmDHat29fPX7GP03EFQU7d+5UkyZNFAwG1bx5c0k//OnVhg0bSpIOHToUsb558+batWuXJKldu3YR5znndPjw4fD1StLVV1+tQCCgSy+9VGPHjlVycrL69OmjxYsX66GHHtKyZct0wQUX8IfKjRFXHZowYUKFTygce1q3bp35HO3atVNiYqLGjx+vX//619q9e7datWql7t27a8SIEXrnnXe0atUq5eXlmc/yU8ZbTurQXXfdpVGjRlW5pkOHDsrKytK+ffuUkJAQ3iMdPXpUBw8elCSlpaVFXGbXrl1q3ry5Nm3apM2bN0ec5/F4lJycLEnKysqSJG3evDn87379+qm4uFiZmZkRMzRt2lQbN27U+eef/2M/XVSDPVcdyszMVE5OTpUnn8+n/v37Kz8/Xzk5OVq8eLEkacmSJQqFQkpKSlJi4j+/54VCIS1evFg/+9nP5PF49Oabb4bPW79+vZxz4d8S1b59e2VlZWnBggXq37+/JGnlypWSpLPOOit8uW3btmnfvn1q0aKF+X3ykxblJ1R+soYNG+batWvnfD6fmzRpkmvbtm34qXi/3+8ef/xxl5qa6k4++WTXsGFDt3PnzvBT8T179nRvvPGGa926tZPkEhMT3SOPPOKeffZZl5OT4yS5++67z02aNMl5vV7n9/tdXl6e+/bbb92iRYvc6aef7k4++WR3+PDhaN8NJzTiipJ9+/a5X/3qV87n8zmPxxN+unzlypXuqaeeci1btnSSXCAQcCtXrnTOOVdcXOwaNWoUfi0rPT3dvfrqq+711193nTp1comJiS45OdklJCQ4j8fjPB6P69SpkzvzzDNdZmamS0pKcm3btnWjR492O3fujPI9cOLjLSeAER5zAUaICzBCXIAR4gKMEBdghLgAI8QFGCEuwAhxAUaICzBCXIAR4opxe/bsUVZWlh588MHwthUrVsjn84XfroLYxA/uxoE//elPGjFihFasWKHOnTurR48eGj58uB577LFoj4YqEFecGDdunBYtWqTevXtr7dq1WrVqFX/VJMYRV5woLi5Wt27dtHXrVq1evVrdu3eP9kioBo+54sSmTZu0Y8cOhUKhcr9HA7GJPVccKC0tVd++fdWjRw917txZ06dP19q1a9WsWbNoj4YqEFccuOeee/TGG2/or3/9q9LS0nTuuecqEAjonXfeifZoqAKHhTEuLy9P06dP16xZs5SRkSGv16tZs2bpww8/1IwZM6I9HqrAngswwp4LMEJcgBHiAowQF2CEuAAjxAUYIS7ACHEBRogLMEJcgBHiAoz8P59XYjHQQJ3MAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"sim = make_sim()\n",
"f, ax = plt.subplots(1, 1, figsize=(10, 10))\n",
"sim.plot(y=0, ax=ax)\n",
"plt.show()\n"
]
},
{
"cell_type": "markdown",
"id": "e4c9a742-5988-444d-9ac4-34b612d42289",
"metadata": {},
"source": [
"Now we write a function to post process some run results to get the transmission we are after."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "1917a8b6-8ba9-4767-90ca-06cd6daadcd3",
"metadata": {},
"outputs": [],
"source": [
"def post_process_T(sim_data: JaxSimulationData) -> float:\n",
" \"\"\"Given some JaxSimulationData from the run, return the transmission of \"p\" polarized light.\"\"\"\n",
" amps = sim_data.output_monitor_data[\"diffraction\"].amps.sel(polarization=\"p\")\n",
" return jnp.sum(abs(amps.values) ** 2)\n"
]
},
{
"cell_type": "markdown",
"id": "b65e4d49-e8b1-4e05-9fce-57bd63b3f2cf",
"metadata": {},
"source": [
"And finally, put everything together in a single function that relates the permittivities and thicknesses of each slab to the transmission, through a `JaxSimulation` run using the adjoint plugin."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "f64ac290-b67d-4a9a-88e8-504de931fef8",
"metadata": {},
"outputs": [],
"source": [
"def compute_T_fdtd(slab_eps=slab_eps0, slab_ds=slab_ds0) -> float:\n",
" \"\"\"Given the slab permittivities and thicknesses, compute T, making sure to use `tidy3d.plugins.adjoint.web.run_adjoint`.\"\"\"\n",
" sim = make_sim(slab_eps=slab_eps, slab_ds=slab_ds)\n",
" sim_data = run_adjoint(sim, task_name=\"slab\", verbose=True)\n",
" return post_process_T(sim_data)\n"
]
},
{
"cell_type": "markdown",
"id": "c4c19d59-d9a3-4935-b785-5fb6b3b3a819",
"metadata": {},
"source": [
"### Computing T and Gradient with FDTD\n",
"\n",
"Now that we have this function defined, we are ready to compute our transmission and gradients using Tidy3d.\n",
"\n",
"We first call `jax.value_and_grad()` on our transmission calculation function, which returns a function that will give us both `T` and the gradient of `T` with respect to the input parameters in one shot. For more details, see the previous tutorial."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "59647cdf-8c2a-40d6-bbc9-b7537616938a",
"metadata": {},
"outputs": [],
"source": [
"compute_T_and_grad_fdtd = jax.value_and_grad(compute_T_fdtd, argnums=(0, 1))\n"
]
},
{
"cell_type": "markdown",
"id": "81caed90-8554-441b-8948-b240e181e98c",
"metadata": {},
"source": [
"Next, we call this function on our starting parameters, which will kick off the original (`fwd`) T transmission simulation and then the reverse (`adj`) simulation, which is used in combination with `fwd` for the gradient calculation."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "ef6716b3-79a3-40cc-ac2a-5056b5ed2ae9",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
[13:52:22] Created task 'slab_fwd' with task_id 'fdve-f4fca853-eac5-4078-8934-3dfd73521510v1'. webapi.py:139\n",
"
[13:52:34] Maximum FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost webapi.py:286\n",
"after a simulation run. \n",
"
\n"
],
"text/plain": [
"\u001b[2;36m[13:52:34]\u001b[0m\u001b[2;36m \u001b[0mMaximum FlexCredit cost: \u001b[1;36m0.025\u001b[0m. Use \u001b[32m'web.real_cost\u001b[0m\u001b[32m(\u001b[0m\u001b[32mtask_id\u001b[0m\u001b[32m)\u001b[0m\u001b[32m'\u001b[0m to get the billed FlexCredit cost \u001b]8;id=904256;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/web/webapi.py\u001b\\\u001b[2mwebapi.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=550109;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/web/webapi.py#286\u001b\\\u001b[2m286\u001b[0m\u001b]8;;\u001b\\\n",
"\u001b[2;36m \u001b[0mafter a simulation run. \u001b[2m \u001b[0m\n"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"
[13:53:11] Maximum FlexCredit cost: 0.025. Use 'web.real_cost(task_id)' to get the billed FlexCredit cost webapi.py:286\n",
"after a simulation run. \n",
"
\n"
],
"text/plain": [
"\u001b[2;36m[13:53:11]\u001b[0m\u001b[2;36m \u001b[0mMaximum FlexCredit cost: \u001b[1;36m0.025\u001b[0m. Use \u001b[32m'web.real_cost\u001b[0m\u001b[32m(\u001b[0m\u001b[32mtask_id\u001b[0m\u001b[32m)\u001b[0m\u001b[32m'\u001b[0m to get the billed FlexCredit cost \u001b]8;id=282903;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/web/webapi.py\u001b\\\u001b[2mwebapi.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=581119;file:///Users/twhughes/Documents/Flexcompute/tidy3d-docs/tidy3d/tidy3d/web/webapi.py#286\u001b\\\u001b[2m286\u001b[0m\u001b]8;;\u001b\\\n",
"\u001b[2;36m \u001b[0mafter a simulation run. \u001b[2m \u001b[0m\n"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"
![](\"img/AdjointMultilayer.png\")
\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "ccbd2c6a-419d-4570-8f91-f57ad5bad813",
"metadata": {},
"outputs": [],
"source": [
"# frequency we want to simulate at\n",
"freq0 = 2.0e14\n",
"k0 = 2 * np.pi * freq0 / td.C_0\n",
"freqs = [freq0]\n",
"wavelength = td.C_0 / freq0\n",
"\n",
"# background permittivity\n",
"bck_eps = 1.3**2\n",
"\n",
"# space between each slab\n",
"spc = 0.0\n",
"\n",
"# slab permittivities and thicknesses\n",
"slab_eps0 = [2.0**2, 1.8**2, 1.5**2, 1.9**2]\n",
"slab_ds0 = [0.5, 0.25, 0.5, 0.5]\n",
"\n",
"# incidence angle\n",
"theta = 0 * np.pi / 8\n",
"\n",
"# resolution\n",
"dl = 0.01\n"
]
},
{
"cell_type": "markdown",
"id": "48dd37b0-ba97-45cd-a861-53ea6ae3fa74",
"metadata": {},
"source": [
"## Transfer Matrix Method (Ground Truth)\n",
"\n",
"Next we use the `tmm` package to write a function to return the transmission `T` of `p` polarized light given a set of slab permittivities and thicknesses. We'll also write a function to compute the numerical gradient using TMM and will take these to be our \"ground truths\" when evaluating the accuracy of our values obtained through FDTD and the adjoint plugin.\n",
"\n",
"### Transmission Calculation with TMM\n",
"\n",
"First, we write a function to compute transmission."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "cc4f6e91-6a53-476e-b4c6-e76e50246426",
"metadata": {},
"outputs": [],
"source": [
"def compute_T_tmm(slab_eps=slab_eps0, slab_ds=slab_ds0) -> float:\n",
" \"\"\"Get transmission as a function of slab permittivities and thicknesses.\"\"\"\n",
"\n",
" # construct lists of permittivities and thicknesses including spaces between\n",
" new_slab_eps = []\n",
" new_slab_ds = []\n",
" for eps, d in zip(slab_eps, slab_ds):\n",
" new_slab_eps.append(eps)\n",
" new_slab_eps.append(bck_eps)\n",
" new_slab_ds.append(d)\n",
" new_slab_ds.append(spc)\n",
" slab_eps = new_slab_eps[:-1]\n",
" slab_ds = new_slab_ds[:-1]\n",
"\n",
" # add the input and output spaces to the lists\n",
" eps_list = [bck_eps] + slab_eps + [bck_eps]\n",
" n_list = np.sqrt(eps_list)\n",
" d_list = [np.inf] + slab_ds + [np.inf]\n",
"\n",
" # compute transmission with TMM\n",
" return tmm.coh_tmm(\"p\", n_list, d_list, theta, wavelength)[\"T\"]\n"
]
},
{
"cell_type": "markdown",
"id": "c6c0da6b-2d2f-4be5-80e6-c84afd8d1a23",
"metadata": {},
"source": [
"We run this function with our starting parameters and see that we get a transmission of about 98% for the set of input parameters."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "c1743e54-d7f3-468a-8ac4-d42d4c2bc8dc",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"T (tmm) = 0.901\n"
]
}
],
"source": [
"T_tmm = compute_T_tmm(slab_eps=slab_eps0, slab_ds=slab_ds0)\n",
"print(f\"T (tmm) = {T_tmm:.3f}\")\n"
]
},
{
"cell_type": "markdown",
"id": "452db467-f486-468d-8ec9-b9e669755c92",
"metadata": {},
"source": [
"### Numerical Gradient with TMM\n",
"\n",
"Next, we will use our `compute_T_tmm()` function to compute the \"numerical\" gradient to use as comparison against our adjoint results with FDTD.\n",
"\n",
"The derivative of a function $f(x)$ w.r.t. $x$ can be approximated using finite differences as\n",
"\n",
"$$ \\frac{df}{dx}(x) \\approx \\frac{f(x+\\Delta) - f(x-\\Delta)}{2\\Delta}$$\n",
"\n",
"with a small step $\\Delta$.\n",
"\n",
"To compute the gradient of our transmission with respect to each of the slab thicknesses and permittivities, we need to repeat this step for each of the values. Luckily, since TMM is very fast, we can compute these quantities quite quickly compared to if we were using FDTD.\n",
"\n",
"> Important note: We assume in our TMM numerical gradient that when the slabs are touching (`spc=0`) and a slab thickness is modified, that the thicknesses of the neighboring slabs adjust to accomidate this change. For example, if slab `i` increases by `dt`, slab `i-1` and `i+1` each decrease by `dt/2`. We also account for this in our FDTD set up by keeping the centers of all boxes constant and not tracking the gradient through these quantities. The reason this is required is that `tidy3d.plugins.adjoint` does not recognize the space between touching `JaxBox` objects as a single interface and will instead \"double count\" the gradient contribution of the interface if they are placed right next to each other. One must therefore be careful about overlapping or touching two `JaxBox` or other geometries when computing gradients.\n",
"\n",
"Here we write the function to return the numerical gradient."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "15d5bae1-6fc7-41c8-8772-2e80505be17d",
"metadata": {},
"outputs": [],
"source": [
"def compute_grad_tmm(\n",
" slab_eps=slab_eps0, slab_ds=slab_ds0\n",
") -> Tuple[List[float], List[float]]:\n",
" \"\"\"Compute numerical gradient of transmission w.r.t. each of the slab permittivities and thicknesses using TMM.\"\"\"\n",
"\n",
" delta = 1e-4\n",
"\n",
" # set up containers to store gradient and perturbed arguments\n",
" num_slabs = len(slab_eps)\n",
" grad_tmm = np.zeros((2, num_slabs), dtype=float)\n",
" args = np.stack((slab_eps, slab_ds), axis=0)\n",
"\n",
" # loop through slab index and argument index (eps, d)\n",
" for arg_index in range(2):\n",
" for slab_index in range(num_slabs):\n",
" grad = 0.0\n",
"\n",
" # perturb the argument by delta in each + and - direction\n",
" for pm in (-1, +1):\n",
" args_num = args.copy()\n",
" args_num[arg_index][slab_index] += delta * pm\n",
"\n",
" # NEW: for slab thickness gradient, need to modify neighboring slabs too\n",
" if arg_index == 1 and spc == 0:\n",
" if slab_index > 0:\n",
" args_num[arg_index][slab_index - 1] -= delta * pm / 2\n",
" if slab_index < num_slabs - 1:\n",
" args_num[arg_index][slab_index + 1] -= delta * pm / 2\n",
"\n",
" # compute argument perturbed T and add to finite difference gradient contribution\n",
" T_tmm = compute_T_tmm(slab_eps=args_num[0], slab_ds=args_num[1])\n",
" grad += pm * T_tmm / 2 / delta\n",
"\n",
" grad_tmm[arg_index][slab_index] = grad\n",
" grad_eps, grad_ds = grad_tmm\n",
" return grad_eps, grad_ds\n"
]
},
{
"cell_type": "markdown",
"id": "d3c362ca-fb31-4ead-8e7c-fcc24d0faecc",
"metadata": {},
"source": [
"Let's run this function and observe the gradients. These will be saved later to compare against our adjoint plugin results."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "b5d1cc93-6d6e-43a7-8e5f-f1829770f3d5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"gradient w.r.t. eps (tmm) = [-0.15463022 0.0376046 -0.0850184 -0.15883104]\n",
"gradient w.r.t. ds (tmm) = [-0.86661161 -0.12292531 0.58010922 -1.05537497]\n"
]
}
],
"source": [
"grad_eps_tmm, grad_ds_tmm = compute_grad_tmm()\n",
"print(f\"gradient w.r.t. eps (tmm) = {grad_eps_tmm}\")\n",
"print(f\"gradient w.r.t. ds (tmm) = {grad_ds_tmm}\")\n"
]
},
{
"cell_type": "markdown",
"id": "50ca68ae-ca5a-49aa-9969-ac6d43324ea9",
"metadata": {},
"source": [
"## FDTD (Using adjoint plugin)\n",
"\n",
"Next, we will implement the same two functions using Tidy3D's adjoint plugin.\n",
"\n",
"### Transmission Calculation with FDTD\n",
"\n",
"We first write a function to compute the transmission of a multilayer slab using Tidy3D.\n",
"\n",
"As discussed in the previous adjoint tutorial notebook, we need to use `jax`-compatible components from the tidy3d subclass for any structures that may depend on the parameters. In this case, this means that the slabs must be `JaxStructures` containing `JaxBox` and `JaxMedium` and must be added to `JaxSimulation.input_structures`.\n",
"\n",
"We use a `DiffractionMonitor` to measure our transmission amplitudes. As the data corresponding to this monitor will be used in the differentiable function return value, we must add it to `JaxSimulation.output_monitors`.\n",
"\n",
"Below, we break up the transmission calculation into a few functions to make it easier to read and re-use later."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "4d12f133-a8d8-4e59-938f-fca4d2e2fa0b",
"metadata": {},
"outputs": [],
"source": [
"def make_sim(slab_eps=slab_eps0, slab_ds=slab_ds0) -> JaxSimulation:\n",
" \"\"\"Create a JaxSimulation given the slab permittivities and thicknesses.\"\"\"\n",
"\n",
" # frequency setup\n",
" wavelength = td.C_0 / freq0\n",
" fwidth = freq0 / 10.0\n",
" freqs = [freq0]\n",
"\n",
" # geometry setup\n",
" bck_medium = td.Medium(permittivity=bck_eps)\n",
"\n",
" space_above = 2\n",
" space_below = 2\n",
"\n",
" length_x = 0.1\n",
" length_y = 0.1\n",
" length_z = space_below + sum(slab_ds0) + space_above + (len(slab_ds0) - 1) * spc\n",
" sim_size = (length_x, length_y, length_z)\n",
"\n",
" # make structures\n",
" slabs = []\n",
" z_start = -length_z / 2 + space_below\n",
" for (d, eps) in zip(slab_ds, slab_eps):\n",
"\n",
" # dont track the gradient through the center of each slab\n",
" # as tidy3d doesn't have enough information to properly process the interface between touching JaxBox objects\n",
" z_center = jax.lax.stop_gradient(z_start + d / 2)\n",
" slab = JaxStructure(\n",
" geometry=JaxBox(center=[0, 0, z_center], size=[td.inf, td.inf, d]),\n",
" medium=JaxMedium(permittivity=eps),\n",
" )\n",
" slabs.append(slab)\n",
" z_start += d + spc\n",
"\n",
" # source setup\n",
" gaussian = td.GaussianPulse(freq0=freq0, fwidth=fwidth)\n",
" src_z = -length_z / 2 + space_below / 2.0\n",
"\n",
" source = td.PlaneWave(\n",
" center=(0, 0, src_z),\n",
" size=(td.inf, td.inf, 0),\n",
" source_time=gaussian,\n",
" direction=\"+\",\n",
" angle_theta=theta,\n",
" angle_phi=0,\n",
" pol_angle=0,\n",
" )\n",
"\n",
" # boundaries\n",
" boundary_x = td.Boundary.bloch_from_source(\n",
" source=source, domain_size=sim_size[0], axis=0, medium=bck_medium\n",
" )\n",
" boundary_y = td.Boundary.bloch_from_source(\n",
" source=source, domain_size=sim_size[1], axis=1, medium=bck_medium\n",
" )\n",
" boundary_spec = td.BoundarySpec(\n",
" x=boundary_x, y=boundary_y, z=td.Boundary.pml(num_layers=40)\n",
" )\n",
"\n",
" # monitors\n",
" mnt_z = length_z / 2 - space_above / 2.0\n",
" monitor_1 = td.DiffractionMonitor(\n",
" center=[0.0, 0.0, mnt_z],\n",
" size=[td.inf, td.inf, 0],\n",
" freqs=freqs,\n",
" name=\"diffraction\",\n",
" normal_dir=\"+\",\n",
" )\n",
"\n",
" # make simulation\n",
" return JaxSimulation(\n",
" size=sim_size,\n",
" grid_spec=td.GridSpec.auto(min_steps_per_wvl=100),\n",
" input_structures=slabs,\n",
" sources=[source],\n",
" output_monitors=[monitor_1],\n",
" run_time=10 / fwidth,\n",
" boundary_spec=boundary_spec,\n",
" medium=bck_medium,\n",
" subpixel=True,\n",
" shutoff=1e-8,\n",
" )\n"
]
},
{
"cell_type": "markdown",
"id": "9addea1f-4dd8-4e5e-a058-b94638fb2d70",
"metadata": {},
"source": [
"Let's generate a simulation and plot it to make sure it looks reasonable."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "bc0dc456-e742-4299-ae5e-22daa49e3997",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"