From a0bc06deeb4e58cd4f8bb53264ecdca88c23fec2 Mon Sep 17 00:00:00 2001
From: Diptodip
Date: Mon, 25 Sep 2023 13:45:41 -0400
Subject: [PATCH 01/76] Add quantization to spatial light modulator
---
src/chromatix/elements/phase_masks.py | 14 +++++++++++---
1 file changed, 11 insertions(+), 3 deletions(-)
diff --git a/src/chromatix/elements/phase_masks.py b/src/chromatix/elements/phase_masks.py
index 2828148..2d61a70 100644
--- a/src/chromatix/elements/phase_masks.py
+++ b/src/chromatix/elements/phase_masks.py
@@ -3,9 +3,10 @@
from chex import Array, PRNGKey
from jax.scipy.ndimage import map_coordinates
from typing import Callable, Optional, Tuple, Union
-from ..field import Field
-from ..functional import wrap_phase, phase_change
-from ..utils import seidel_aberrations, zernike_aberrations
+from chromatix.field import Field
+from chromatix.functional import wrap_phase, phase_change
+from chromatix.utils import seidel_aberrations, zernike_aberrations
+from chromatix.ops import quantize
from chromatix.elements.utils import register
__all__ = [
@@ -107,6 +108,10 @@ class SpatialLightModulator(nn.Module):
spacing: The pitch of the SLM pixels.
phase_range: The phase range that the SLM can simulate, provided as
(min, max).
+ num_bits: The number of bits of precision the phase pixels should be
+ quantized to. Defaults to None, in which case no quantization is
+ applied. Otherwise, the phase will be quantized to have
+ ``2.0 ** num_bits`` values within ``phase_range``.
interpolation_order: The order of interpolation for the SLM pixels to
the shape of the incoming ``Field``. Can be 0 or 1. Defaults to 0.
f: Focal length of the system's objective. Defaults to None.
@@ -118,6 +123,7 @@ class SpatialLightModulator(nn.Module):
shape: Tuple[int, int]
spacing: float
phase_range: Tuple[float, float]
+ num_bits: Optional[Union[int, float]] = None
interpolation_order: int = 0
f: Optional[float] = None
n: Optional[float] = None
@@ -143,6 +149,8 @@ def __call__(self, field: Field) -> Field:
phase.shape == self.shape
), "Provided phase shape should match provided SLM shape"
phase = wrap_phase(phase, self.phase_range)
+ if self.num_bits is not None:
+ phase = quantize(phase, 2.0**self.num_bits, range=self.phase_range)
field_pixel_grid = jnp.meshgrid(
jnp.linspace(0, self.shape[0] - 1, num=field.spatial_shape[0]) + 0.5,
jnp.linspace(0, self.shape[1] - 1, num=field.spatial_shape[1]) + 0.5,
From 30b4d9c7edcd06aaf3c31ec9a4dc67ec2fae1854 Mon Sep 17 00:00:00 2001
From: Eric Bezzam
Date: Fri, 24 May 2024 16:25:11 +0000
Subject: [PATCH 02/76] Add support for bandlimited angular spectrum.
---
examples/bandlimited_angular_spectrum.py | 115 +++++++++++++++++++++++
src/chromatix/functional/propagation.py | 64 ++++++++++++-
2 files changed, 177 insertions(+), 2 deletions(-)
create mode 100644 examples/bandlimited_angular_spectrum.py
diff --git a/examples/bandlimited_angular_spectrum.py b/examples/bandlimited_angular_spectrum.py
new file mode 100644
index 0000000..50ee476
--- /dev/null
+++ b/examples/bandlimited_angular_spectrum.py
@@ -0,0 +1,115 @@
+"""
+Example of "Band-Limited Angular Spectrum Method for Numerical Simulation
+of Free-Space Propagation in Far and Near Fields" (2010) by Matsushima and
+Shimobaba.
+
+Specifically trying to replicate Fig 9a from the paper for a rectangular
+aperture.
+
+TODO: implement numerical integration for comparison?
+Something like this: https://github.com/ebezzam/waveprop/blob/a2d65116336bfb6e95732fd982e5c3ec2109cff3/waveprop/rs.py#L33
+
+"""
+from functools import partial
+import numpy as np
+import jax.numpy as jnp
+from scipy.special import fresnel
+import chromatix.functional as cf
+import matplotlib.pyplot as plt
+
+
+# setting like in BLAS paper (Fig 9) https://opg.optica.org/oe/fulltext.cfm?uri=oe-17-22-19662&id=186848
+shape = (1024, 1024)
+N_pad = (512, 512)
+spectrum = 0.532 # wavelength in microns
+dxi = 2 * spectrum
+D = dxi * shape[0] # field shape in microns
+w = D / 2
+z = 100 * D
+
+dxi = D / np.array(shape)
+spacing = dxi[..., np.newaxis]
+n = 1 # refractive index of medium
+
+# # setting like https://github.com/chromatix-team/chromatix/blob/7304cd312b28eebc2f15c3c466e53074141d553b/tests/test_propagate.py#L34C1-L52C28
+# D = 40 # microns
+# z = 100 # microns
+# spectrum = 0.532 # microns
+# shape = (512, 512)
+# N_pad = (512, 512)
+# n = 1 # refractive index of medium
+# dxi = D / np.array(shape)
+# spacing = dxi[..., np.newaxis]
+# w = dxi[1] * shape[1] # width of aperture in microns
+
+print("Field shape [um]: ", D)
+print("Width of aperture [um]: ", w)
+print("Propagation distance [um]: ", z)
+
+
+def analytical_result_square_aperture(x, z, D, spectrum, n):
+ # TODO: this uses Fresnel approximation
+ Nf = (D / 2) ** 2 / (spectrum / n * z)
+
+ def I(x):
+ Smin, Cmin = fresnel(jnp.sqrt(2 * Nf) * (1 - 2 * x / D))
+ Splus, Cplus = fresnel(jnp.sqrt(2 * Nf) * (1 + 2 * x / D))
+
+ return 1 / jnp.sqrt(2) * (Cmin + Cplus) + 1j / jnp.sqrt(2) * (Smin + Splus)
+
+ U = jnp.exp(1j * 2 * jnp.pi * z * n / spectrum) / 1j * I(x[0]) * I(x[1])
+ # Return U/l as the input field has area l^2
+ return U / D
+
+# Input field
+field = cf.plane_wave(
+ shape=shape,
+ dx=spacing,
+ spectrum=spectrum,
+ spectral_density=1.0,
+ pupil=partial(cf.square_pupil, w=w)
+)
+
+# # Fresnel
+# out_field_fresnel = cf.transform_propagate(field, z, n, N_pad=N_pad)
+# I_fresnel = out_field_fresnel.intensity.squeeze()
+
+# # Analytical (Fresnel)
+# xi = np.array(out_field_fresnel.grid.squeeze())
+# U_analytical = analytical_result_square_aperture(xi, z, D, spectrum, n)
+# I_analytical = jnp.abs(U_analytical) ** 2
+
+# Angular spectrum
+out_field_asm = cf.asm_propagate(field, z, n, N_pad=N_pad, mode="same")
+I_asm = out_field_asm.intensity.squeeze()
+
+# Angular spectrum (bandlimited)
+out_field_blas = cf.asm_propagate(field, z, n, N_pad=N_pad, mode="same", bandlimit=True)
+I_blas = out_field_blas.intensity.squeeze()
+
+# Compare
+# -- compute error
+intensities = [
+ ["Input", field.intensity.squeeze()],
+ # ["Analytical (Fresnel)", I_analytical],
+ # ["Fresnel", I_fresnel],
+ ["ASM", I_asm],
+ ["BLAS", I_blas],
+]
+# for approach, intensity in intensities[2:]:
+# rel_error = jnp.mean((I_analytical - intensity) ** 2) / jnp.mean(
+# I_analytical**2
+# )
+# print(f"{approach} error: ", rel_error)
+
+# -- plot
+fig, axs = plt.subplots(1, len(intensities), figsize=(15, 4))
+axs[0].set_ylabel("y (microns)")
+for ax, (title, intensity) in zip(axs, intensities):
+ ax.imshow(intensity, cmap="gray", extent=[-D/2, D/2, -D/2, D/2])
+ ax.set_title(title)
+ ax.set_xlabel("x (microns)")
+
+plot_fn = "propagation_comparison.png"
+plt.savefig(plot_fn)
+print(f"Saved plot to {plot_fn}")
diff --git a/src/chromatix/functional/propagation.py b/src/chromatix/functional/propagation.py
index b48989e..dc93ae8 100644
--- a/src/chromatix/functional/propagation.py
+++ b/src/chromatix/functional/propagation.py
@@ -141,6 +141,7 @@ def asm_propagate(
N_pad: int,
cval: float = 0,
kykx: Union[Array, Tuple[float, float]] = (0.0, 0.0),
+ bandlimit: bool = False,
mode: Literal["full", "same"] = "full",
) -> Field:
"""
@@ -161,12 +162,15 @@ def asm_propagate(
for zero padding.
kykx: If provided, defines the orientation of the propagation. Should
be an array of shape `[2,]` in the format [ky, kx].
+ bandlimit: If provided, bandlimited the kernel according to "Band-Limited
+ Angular Spectrum Method for Numerical Simulation of Free-Space
+ Propagation in Far and Near Fields" (2009) by Matsushima and Shimobaba.
mode: Either "full" or "same". If "same", the shape of the output
``Field`` will match the shape of the incoming ``Field``. Defaults
to "full", in which case the output shape will include padding.
"""
field = pad(field, N_pad, cval=cval)
- propagator = compute_asm_propagator(field, z, n, kykx)
+ propagator = compute_asm_propagator(field, z, n, kykx, bandlimit)
field = kernel_propagate(field, propagator)
if mode == "same":
field = crop(field, N_pad)
@@ -243,6 +247,7 @@ def compute_asm_propagator(
z: Union[float, Array],
n: float,
kykx: Union[Array, Tuple[float, float]] = (0.0, 0.0),
+ bandlimit: bool = False,
) -> Array:
"""
Compute propagation kernel for propagation with no Fresnel approximation.
@@ -266,7 +271,62 @@ def compute_asm_propagator(
delay = jnp.sqrt(jnp.abs(kernel))
delay = jnp.where(kernel >= 0, delay, 1j * delay) # keep evanescent modes
phase = 2 * jnp.pi * (z * n / field.spectrum) * delay
- return jnp.fft.ifftshift(jnp.exp(1j * phase), axes=field.spatial_dims)
+ kernel_field = jnp.exp(1j * phase)
+
+ if bandlimit:
+ Sy, Sx = (1 / field.dk).squeeze() # spatial dimension in microns
+ y0, x0 = (kykx / field.dk).squeeze() # spatial shift in microns, TODO check
+ z0 = z.squeeze() # propagation distance in microns
+ wv = field.spectrum.squeeze() # wavelength in microns
+
+ dfX = 1.0 / Sx
+ dfY = 1.0 / Sy
+ N_y, N_x = field.spatial_shape
+ fX = np.linspace(-N_x // 2, N_x // 2 - 1, num=N_x)[np.newaxis, :] * dfX
+ fY = np.linspace(-N_y // 2, N_y // 2 - 1, num=N_y)[:, np.newaxis] * dfY
+
+ # Table 1 of "Shifted angular spectrum method for off-axis numerical
+ # propagation" (2010) by Matsushima
+ du = 1 / (2 * Sx)
+ u_limit_p = ((x0 + 1 / (2 * du)) ** (-2) * z0**2 + 1) ** (-1 / 2) / wv
+ u_limit_n = ((x0 - 1 / (2 * du)) ** (-2) * z0**2 + 1) ** (-1 / 2) / wv
+ if Sx < x0:
+ u0 = (u_limit_p + u_limit_n) / 2
+ u_width = u_limit_p - u_limit_n
+ elif x0 <= -Sx:
+ u0 = -(u_limit_p + u_limit_n) / 2
+ u_width = u_limit_n - u_limit_p
+ else:
+ u0 = (u_limit_p - u_limit_n) / 2
+ u_width = u_limit_p + u_limit_n
+
+ dv = 1 / (2 * Sy)
+ v_limit_p = ((y0 + 1 / (2 * dv)) ** (-2) * z0**2 + 1) ** (-1 / 2) / wv
+ v_limit_n = ((y0 - 1 / (2 * dv)) ** (-2) * z0**2 + 1) ** (-1 / 2) / wv
+ if Sy < y0:
+ v0 = (v_limit_p + v_limit_n) / 2
+ v_width = v_limit_p - v_limit_n
+ elif y0 <= -Sy:
+ v0 = -(v_limit_p + v_limit_n) / 2
+ v_width = v_limit_n - v_limit_p
+ else:
+ v0 = (v_limit_p - v_limit_n) / 2
+ v_width = v_limit_p + v_limit_n
+
+ fx_max = u_width / 2
+ fy_max = v_width / 2
+
+ # bandlimit
+ H_filter = (np.abs(fX - u0) <= fx_max) * (np.abs(fY - v0) < fy_max)
+
+ # to jax
+ H_filter = jnp.array(H_filter)
+ H_filter = jnp.expand_dims(H_filter, (0, 3, 4))
+
+ # apply filter
+ kernel_field = kernel_field * H_filter
+
+ return jnp.fft.ifftshift(kernel_field, axes=field.spatial_dims)
def compute_padding_transform(height: int, spectrum: float, dx: float, z: float) -> int:
From 43cc2bb6a872e5ee69444ddf23834bcbc6dd5b2b Mon Sep 17 00:00:00 2001
From: Geneva Schlafly
Date: Fri, 24 May 2024 13:51:28 -0500
Subject: [PATCH 03/76] Add universal compensator polarizer
---
src/chromatix/functional/polarizers.py | 18 ++++++++++++++----
1 file changed, 14 insertions(+), 4 deletions(-)
diff --git a/src/chromatix/functional/polarizers.py b/src/chromatix/functional/polarizers.py
index a9eec4c..bfa2f9c 100644
--- a/src/chromatix/functional/polarizers.py
+++ b/src/chromatix/functional/polarizers.py
@@ -16,6 +16,7 @@
"linear_polarizer",
"left_circular_polarizer",
"right_circular_polarizer",
+ "universal_compensator",
# Waveplates
"wave_plate",
"halfwave_plate",
@@ -112,8 +113,8 @@ def linear_polarizer(field: VectorField, angle: float) -> VectorField:
"""
c, s = jnp.cos(angle), jnp.sin(angle)
- J00 = c**2
- J11 = s**2
+ J00 = c ** 2
+ J11 = s ** 2
J01 = s * c
J10 = J01
return polarizer(field, J00, J01, J10, J11)
@@ -169,8 +170,8 @@ def phase_retarder(
"""
s, c = jnp.sin(theta), jnp.cos(theta)
scale = jnp.exp(-1j * eta / 2)
- J00 = scale * (c**2 + jnp.exp(1j * eta) * s**2)
- J11 = scale * (s**2 + jnp.exp(1j * eta) * c**2)
+ J00 = scale * (c ** 2 + jnp.exp(1j * eta) * s ** 2)
+ J11 = scale * (s ** 2 + jnp.exp(1j * eta) * c ** 2)
J01 = scale * (1 - jnp.exp(1j * eta)) * jnp.exp(-1j * phi) * s * c
J10 = scale * (1 - jnp.exp(1j * eta)) * jnp.exp(1j * phi) * s * c
return polarizer(field, J00, J01, J10, J11)
@@ -214,3 +215,12 @@ def quarterwave_plate(field: VectorField, theta: float) -> VectorField:
VectorField: outgoing field.
"""
return phase_retarder(field, theta, eta=jnp.pi / 2, phi=0)
+
+
+def universal_compensator(field, retA, retB):
+ """Universal Polarizer for the LC-PolScope"""
+ field_LP = linear_polarizer(field, 0)
+ field_retA = phase_retarder(field_LP, -jnp.pi / 4, retA, 0)
+ field_retB = phase_retarder(field_retA, 0, retB, 0)
+ x = 5
+ return field_retB
From 512186dd46034a093b7e0576cf784cc54ce55170 Mon Sep 17 00:00:00 2001
From: Geneva Schlafly
Date: Fri, 24 May 2024 13:56:36 -0500
Subject: [PATCH 04/76] Document universal compensator
---
src/chromatix/functional/polarizers.py | 22 +++++++++++++++-------
1 file changed, 15 insertions(+), 7 deletions(-)
diff --git a/src/chromatix/functional/polarizers.py b/src/chromatix/functional/polarizers.py
index bfa2f9c..93967e3 100644
--- a/src/chromatix/functional/polarizers.py
+++ b/src/chromatix/functional/polarizers.py
@@ -113,8 +113,8 @@ def linear_polarizer(field: VectorField, angle: float) -> VectorField:
"""
c, s = jnp.cos(angle), jnp.sin(angle)
- J00 = c ** 2
- J11 = s ** 2
+ J00 = c**2
+ J11 = s**2
J01 = s * c
J10 = J01
return polarizer(field, J00, J01, J10, J11)
@@ -170,8 +170,8 @@ def phase_retarder(
"""
s, c = jnp.sin(theta), jnp.cos(theta)
scale = jnp.exp(-1j * eta / 2)
- J00 = scale * (c ** 2 + jnp.exp(1j * eta) * s ** 2)
- J11 = scale * (s ** 2 + jnp.exp(1j * eta) * c ** 2)
+ J00 = scale * (c**2 + jnp.exp(1j * eta) * s**2)
+ J11 = scale * (s**2 + jnp.exp(1j * eta) * c**2)
J01 = scale * (1 - jnp.exp(1j * eta)) * jnp.exp(-1j * phi) * s * c
J10 = scale * (1 - jnp.exp(1j * eta)) * jnp.exp(1j * phi) * s * c
return polarizer(field, J00, J01, J10, J11)
@@ -217,10 +217,18 @@ def quarterwave_plate(field: VectorField, theta: float) -> VectorField:
return phase_retarder(field, theta, eta=jnp.pi / 2, phi=0)
-def universal_compensator(field, retA, retB):
- """Universal Polarizer for the LC-PolScope"""
+def universal_compensator(field: VectorField, retA: float, retB: float) -> VectorField:
+ """Applies the Universal Polarizer for the LC-PolScope to the incoming field.
+
+ Args:
+ field (VectorField): incoming field.
+ retA (float): retardance induces at a 45 deg angle.
+ retB (float): retardance induces at a 0 deg angle.
+
+ Returns:
+ VectorField: outgoing field.
+ """
field_LP = linear_polarizer(field, 0)
field_retA = phase_retarder(field_LP, -jnp.pi / 4, retA, 0)
field_retB = phase_retarder(field_retA, 0, retB, 0)
- x = 5
return field_retB
From c65709513d68fc5d833db2471425562bfec1e0c4 Mon Sep 17 00:00:00 2001
From: Geneva Schlafly
Date: Fri, 24 May 2024 14:19:23 -0500
Subject: [PATCH 05/76] Black format multiplication
---
src/chromatix/functional/polarizers.py | 8 ++++----
1 file changed, 4 insertions(+), 4 deletions(-)
diff --git a/src/chromatix/functional/polarizers.py b/src/chromatix/functional/polarizers.py
index 93967e3..8966254 100644
--- a/src/chromatix/functional/polarizers.py
+++ b/src/chromatix/functional/polarizers.py
@@ -113,8 +113,8 @@ def linear_polarizer(field: VectorField, angle: float) -> VectorField:
"""
c, s = jnp.cos(angle), jnp.sin(angle)
- J00 = c**2
- J11 = s**2
+ J00 = c ** 2
+ J11 = s ** 2
J01 = s * c
J10 = J01
return polarizer(field, J00, J01, J10, J11)
@@ -170,8 +170,8 @@ def phase_retarder(
"""
s, c = jnp.sin(theta), jnp.cos(theta)
scale = jnp.exp(-1j * eta / 2)
- J00 = scale * (c**2 + jnp.exp(1j * eta) * s**2)
- J11 = scale * (s**2 + jnp.exp(1j * eta) * c**2)
+ J00 = scale * (c ** 2 + jnp.exp(1j * eta) * s ** 2)
+ J11 = scale * (s ** 2 + jnp.exp(1j * eta) * c ** 2)
J01 = scale * (1 - jnp.exp(1j * eta)) * jnp.exp(-1j * phi) * s * c
J10 = scale * (1 - jnp.exp(1j * eta)) * jnp.exp(1j * phi) * s * c
return polarizer(field, J00, J01, J10, J11)
From 6de1c252e5a895be28856bbfe4f975269214f9a9 Mon Sep 17 00:00:00 2001
From: Geneva Schlafly
Date: Fri, 24 May 2024 14:33:30 -0500
Subject: [PATCH 06/76] Create permittivity tensor phantom
---
src/chromatix/data/permittivity_tensors.py | 70 ++++++++++++++++++++++
src/chromatix/functional/polarizers.py | 8 +--
src/chromatix/utils/data.py | 3 +-
3 files changed, 75 insertions(+), 6 deletions(-)
create mode 100644 src/chromatix/data/permittivity_tensors.py
diff --git a/src/chromatix/data/permittivity_tensors.py b/src/chromatix/data/permittivity_tensors.py
new file mode 100644
index 0000000..2802fd4
--- /dev/null
+++ b/src/chromatix/data/permittivity_tensors.py
@@ -0,0 +1,70 @@
+import jax.numpy as jnp
+
+
+def generate_permittivity_tensor(n_o, n_e, extraordinary_axis="z"):
+ """
+ Generate the permittivity tensor for a uniaxial anisotropic material.
+
+ Args:
+ n_o (float): Ordinary refractive index
+ n_e (float): Extraordinary refractive index
+ extraordinary_axis (str): Axis which is extraordinary ('x', 'y', or 'z')
+
+ Returns:
+ jnp.ndarray: Permittivity tensor with the order of axes as zyx
+ """
+ epsilon_o = n_o**2
+ epsilon_e = n_e**2
+
+ if extraordinary_axis == "x":
+ epsilon_tensor = jnp.array(
+ [[epsilon_e, 0, 0], [0, epsilon_o, 0], [0, 0, epsilon_o]]
+ )
+ elif extraordinary_axis == "y":
+ epsilon_tensor = jnp.array(
+ [[epsilon_o, 0, 0], [0, epsilon_e, 0], [0, 0, epsilon_o]]
+ )
+ elif extraordinary_axis == "z":
+ epsilon_tensor = jnp.array(
+ [[epsilon_o, 0, 0], [0, epsilon_o, 0], [0, 0, epsilon_e]]
+ )
+ else:
+ raise ValueError("extraordinary_axis must be one of 'x', 'y', or 'z'")
+
+ return epsilon_tensor
+
+
+def create_homogeneous_phantom(shape, n_o, n_e, extraordinary_axis="z"):
+ """
+ Create a homogeneous uniaxial anisotropic phantom.
+
+ Args:
+ shape (tuple): Shape of the phantom (z, y, x)
+ n_o (float): Ordinary refractive index
+ n_e (float): Extraordinary refractive index
+ extraordinary_axis (str): Axis which is extraordinary ('x', 'y', or 'z')
+
+ Returns:
+ jnp.ndarray: 4D array representing the phantom with the
+ permittivity tensor at each voxel
+ """
+ epsilon_tensor = generate_permittivity_tensor(n_o, n_e, extraordinary_axis)
+ phantom = jnp.tile(epsilon_tensor, (*shape, 1, 1))
+ return phantom
+
+
+def create_calcite_crystal(shape, extraordinary_axis="z"):
+ """
+ Create a calcite crystal phantom.
+
+ Args:
+ shape (tuple): Shape of the phantom (z, y, x)
+ extraordinary_axis (str): Axis which is extraordinary ('x', 'y', or 'z')
+
+ Returns:
+ jnp.ndarray: 4D array representing the phantom with the
+ permittivity tensor at each voxel
+ """
+ n_o = 1.658
+ n_e = 1.486
+ return create_homogeneous_phantom(shape, n_o, n_e, extraordinary_axis)
diff --git a/src/chromatix/functional/polarizers.py b/src/chromatix/functional/polarizers.py
index 8966254..93967e3 100644
--- a/src/chromatix/functional/polarizers.py
+++ b/src/chromatix/functional/polarizers.py
@@ -113,8 +113,8 @@ def linear_polarizer(field: VectorField, angle: float) -> VectorField:
"""
c, s = jnp.cos(angle), jnp.sin(angle)
- J00 = c ** 2
- J11 = s ** 2
+ J00 = c**2
+ J11 = s**2
J01 = s * c
J10 = J01
return polarizer(field, J00, J01, J10, J11)
@@ -170,8 +170,8 @@ def phase_retarder(
"""
s, c = jnp.sin(theta), jnp.cos(theta)
scale = jnp.exp(-1j * eta / 2)
- J00 = scale * (c ** 2 + jnp.exp(1j * eta) * s ** 2)
- J11 = scale * (s ** 2 + jnp.exp(1j * eta) * c ** 2)
+ J00 = scale * (c**2 + jnp.exp(1j * eta) * s**2)
+ J11 = scale * (s**2 + jnp.exp(1j * eta) * c**2)
J01 = scale * (1 - jnp.exp(1j * eta)) * jnp.exp(-1j * phi) * s * c
J10 = scale * (1 - jnp.exp(1j * eta)) * jnp.exp(1j * phi) * s * c
return polarizer(field, J00, J01, J10, J11)
diff --git a/src/chromatix/utils/data.py b/src/chromatix/utils/data.py
index f08d2e6..a8680d8 100644
--- a/src/chromatix/utils/data.py
+++ b/src/chromatix/utils/data.py
@@ -76,8 +76,7 @@ def draw_disks(
image = np.zeros([s + radius * 2 for s in shape], dtype=np.uint8)
_samples = np.linspace(-radius, radius, num=radius * 2, dtype=np.float32)
circle = color * np.uint8(
- np.sum(np.array(np.meshgrid(_samples, _samples)) ** 2, axis=0)
- <= radius**2
+ np.sum(np.array(np.meshgrid(_samples, _samples)) ** 2, axis=0) <= radius**2
)
for c in coordinates:
slices = (slice(c[0], c[0] + radius * 2), slice(c[1], c[1] + radius * 2))
From 807ea0961821d35c06ccba8119e8dac99b0abfee Mon Sep 17 00:00:00 2001
From: Geneva Schlafly
Date: Fri, 24 May 2024 15:15:41 -0500
Subject: [PATCH 07/76] Add typing to parameters
---
src/chromatix/data/permittivity_tensors.py | 22 +++++++++++++++-------
1 file changed, 15 insertions(+), 7 deletions(-)
diff --git a/src/chromatix/data/permittivity_tensors.py b/src/chromatix/data/permittivity_tensors.py
index 2802fd4..52e1591 100644
--- a/src/chromatix/data/permittivity_tensors.py
+++ b/src/chromatix/data/permittivity_tensors.py
@@ -1,7 +1,10 @@
import jax.numpy as jnp
+from typing import Optional, Tuple
-def generate_permittivity_tensor(n_o, n_e, extraordinary_axis="z"):
+def generate_permittivity_tensor(
+ n_o: float, n_e: float, extraordinary_axis: Optional[str] = "x"
+):
"""
Generate the permittivity tensor for a uniaxial anisotropic material.
@@ -15,8 +18,7 @@ def generate_permittivity_tensor(n_o, n_e, extraordinary_axis="z"):
"""
epsilon_o = n_o**2
epsilon_e = n_e**2
-
- if extraordinary_axis == "x":
+ if extraordinary_axis == "z":
epsilon_tensor = jnp.array(
[[epsilon_e, 0, 0], [0, epsilon_o, 0], [0, 0, epsilon_o]]
)
@@ -24,17 +26,21 @@ def generate_permittivity_tensor(n_o, n_e, extraordinary_axis="z"):
epsilon_tensor = jnp.array(
[[epsilon_o, 0, 0], [0, epsilon_e, 0], [0, 0, epsilon_o]]
)
- elif extraordinary_axis == "z":
+ elif extraordinary_axis == "x":
epsilon_tensor = jnp.array(
[[epsilon_o, 0, 0], [0, epsilon_o, 0], [0, 0, epsilon_e]]
)
else:
raise ValueError("extraordinary_axis must be one of 'x', 'y', or 'z'")
-
return epsilon_tensor
-def create_homogeneous_phantom(shape, n_o, n_e, extraordinary_axis="z"):
+def create_homogeneous_phantom(
+ shape: Tuple[int, int, int],
+ n_o: float,
+ n_e: float,
+ extraordinary_axis: Optional[str] = "x",
+):
"""
Create a homogeneous uniaxial anisotropic phantom.
@@ -53,7 +59,9 @@ def create_homogeneous_phantom(shape, n_o, n_e, extraordinary_axis="z"):
return phantom
-def create_calcite_crystal(shape, extraordinary_axis="z"):
+def create_calcite_crystal(
+ shape: Tuple[int, int, int], extraordinary_axis: Optional[str] = "z"
+):
"""
Create a calcite crystal phantom.
From f275e0b971ef193e25ec5fbbaff97efea5c68321 Mon Sep 17 00:00:00 2001
From: Geneva Schlafly
Date: Fri, 24 May 2024 15:27:31 -0500
Subject: [PATCH 08/76] Use waveplate instead of general phase retarder
---
src/chromatix/functional/polarizers.py | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/src/chromatix/functional/polarizers.py b/src/chromatix/functional/polarizers.py
index 93967e3..577b73c 100644
--- a/src/chromatix/functional/polarizers.py
+++ b/src/chromatix/functional/polarizers.py
@@ -229,6 +229,6 @@ def universal_compensator(field: VectorField, retA: float, retB: float) -> Vecto
VectorField: outgoing field.
"""
field_LP = linear_polarizer(field, 0)
- field_retA = phase_retarder(field_LP, -jnp.pi / 4, retA, 0)
- field_retB = phase_retarder(field_retA, 0, retB, 0)
+ field_retA = wave_plate(field_LP, -jnp.pi / 4, retA)
+ field_retB = wave_plate(field_retA, 0, retB)
return field_retB
From 7d74fd0f97651b67ad0a2843a37824ae63c4dc44 Mon Sep 17 00:00:00 2001
From: Geneva Schlafly
Date: Fri, 24 May 2024 17:57:30 -0400
Subject: [PATCH 09/76] Add matrix mult Field class method
Essential for creating a field from the output of one field times a jones or permittivity matrix
---
src/chromatix/field.py | 3 +++
1 file changed, 3 insertions(+)
diff --git a/src/chromatix/field.py b/src/chromatix/field.py
index cc0a552..ef7fca9 100644
--- a/src/chromatix/field.py
+++ b/src/chromatix/field.py
@@ -230,6 +230,9 @@ def __mul__(self, other: Union[Number, jnp.ndarray, Field]) -> Field:
else:
return NotImplemented
+ def __matmul__(self, other: jnp.array) -> Field:
+ return self.replace(u=jnp.matmul(self.u, other))
+
def __rmul__(self, other: Any) -> Field:
return self * other
From 41917e7538961a3049449bf0944164fb517cfa78 Mon Sep 17 00:00:00 2001
From: Rainer Heintzmann
Date: Sat, 25 May 2024 13:55:07 +0200
Subject: [PATCH 10/76] towards sas propagation
---
docs/examples/sas_propagation_chr.ipynb | 479 ++++++++++++++++++++++++
docs/examples/sas_propagation_jax.ipynb | 301 +++++++++++++++
docs/examples/seidel_fitting.ipynb | 96 ++---
docs/examples/zernike_fitting.ipynb | 85 ++---
src/chromatix/functional/propagation.py | 77 +++-
5 files changed, 939 insertions(+), 99 deletions(-)
create mode 100644 docs/examples/sas_propagation_chr.ipynb
create mode 100644 docs/examples/sas_propagation_jax.ipynb
diff --git a/docs/examples/sas_propagation_chr.ipynb b/docs/examples/sas_propagation_chr.ipynb
new file mode 100644
index 0000000..6315373
--- /dev/null
+++ b/docs/examples/sas_propagation_chr.ipynb
@@ -0,0 +1,479 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import jax.numpy as jnp\n",
+ "import jax as jax\n",
+ "import matplotlib.pyplot as plt\n",
+ "from colorsys import hls_to_rgb\n",
+ "import matplotlib.pyplot as plt\n",
+ "from jax.numpy import pi\n",
+ "import chromatix.functional as cx"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# CC nadapez: from https://stackoverflow.com/a/20958684\n",
+ "def colorize(z):\n",
+ " r = np.abs(z)\n",
+ " arg = np.angle(z) \n",
+ "\n",
+ " h = (arg + pi) / (2 * pi) + 0.5\n",
+ " l = 1.0 - 1.0/(1.0 + r**0.3)\n",
+ " s = 0.8\n",
+ "\n",
+ " c = np.vectorize(hls_to_rgb) (h,l,s) # --> tuple\n",
+ " c = np.array(c) # --> array of (3,n,m) shape, but need (n,m,3)\n",
+ " c = c.swapaxes(0,2) \n",
+ " return c"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Scalable Angular Spectrum Method\n",
+ "General hint, we assume the array to be three dimensional. The first dimension is a batch dimension.\n",
+ "\n",
+ "License\n",
+ "If you copy this code, include this LICENSE statement:\n",
+ "\n",
+ "MIT License. Copyright (c) 2023 Felix Wechsler (info@felixwechsler.science), Rainer Heintzmann, Lars Lötgering\n",
+ "\n",
+ "This notebook is based on https://github.com/bionanoimaging/Scalable-Angular-Spectrum-Method-SAS/blob/main/SAS_pytorch.ipynb"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def zero_pad(arr):\n",
+ " '''\n",
+ " Pad arr with zeros to double the size. First dim is assumed to be batch dim which\n",
+ " won't be changed\n",
+ " '''\n",
+ " N_pad = ((0,0), (0, arr.shape[1]), (0, arr.shape[2])) # expands the shape *2 in X and Y\n",
+ " return jnp.pad(arr, N_pad, constant_values=0)\n",
+ " # out_arr = jnp.zeros((arr.shape[0], arr.shape[1] * 2, arr.shape[2] * 2), dtype=arr.dtype)\n",
+ " \n",
+ " # as1 = (arr.shape[1] + 1) // 2\n",
+ " # as2 = (arr.shape[2] + 1) // 2\n",
+ " # out_arr[:, as1:as1 + arr.shape[1], as2:as2 + arr.shape[2]] = arr\n",
+ " # return out_arr\n",
+ "\n",
+ "def zero_unpad(arr, original_shape):\n",
+ " '''\n",
+ " Strip off padding of arr with zeros to halve the size. First dim is assumed to be batch dim which\n",
+ " won't be changed\n",
+ " '''\n",
+ " as1 = (original_shape[1] + 1) // 2\n",
+ " as2 = (original_shape[2] + 1) // 2\n",
+ " return arr[:, as1:as1 + original_shape[1], as2:as2 + original_shape[2]]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "@jax.jit\n",
+ "def scalable_angular_spectrum(psi, z, lbd, L, skip_final_phase=True):\n",
+ " '''\n",
+ " Returns the complex electrical field psi propagated with the Scalable Angular Spectrum Method.\n",
+ " \n",
+ " Parameters:\n",
+ " psi (torch.tensor): the quadratically shaped input field, with leading batch dimension\n",
+ " z (number): propagation distance\n",
+ " lbd (number): vacuum wavelength\n",
+ " L (number): physical sidelength of the input field\n",
+ " skip_final_phase=True: Skip final multiplication of phase factor. For M>2 undersampled,\n",
+ " \n",
+ " Returns:\n",
+ " psi_final (torch.tensor): Propagated field\n",
+ " Q (number): Output field size, corresponds to magnificiation * L\n",
+ " \n",
+ " '''\n",
+ " N = psi.shape[-1]\n",
+ " z_limit = (- 4 * L * jnp.sqrt(8*L**2 / N**2 + lbd**2) * jnp.sqrt(L**2 * 1 / (8 * L**2 + N**2 * lbd**2))\\\n",
+ " / (lbd * (-1+2 * jnp.sqrt(2) * jnp.sqrt(L**2 * 1 / (8 * L**2 + N**2 * lbd**2)))))\n",
+ " \n",
+ " # assert z <= z_limit\n",
+ " \n",
+ " \n",
+ " # don't change this pad_factor, only 2 is supported\n",
+ " pad_factor = 2\n",
+ " L_new = pad_factor * L\n",
+ " N_new = pad_factor * N\n",
+ " # pad array\n",
+ " M = lbd * z * N / L**2 / 2\n",
+ " psi_p = zero_pad(psi)\n",
+ " \n",
+ " # helper varaibles\n",
+ " k = 2 * jnp.pi / lbd\n",
+ " df = 1 / L_new \n",
+ " Lf = N_new * df\n",
+ " \n",
+ " # freq space coordinates for padded array\n",
+ " f_y = jnp.fft.fftfreq(N_new, 1 / Lf, dtype=jnp.float32).reshape(1,1, N_new)\n",
+ " f_x = f_y.reshape(1, N_new, 1)\n",
+ " \n",
+ " # real space coordinates for padded array\n",
+ " y = jnp.fft.ifftshift(jnp.linspace(-L_new/2, L_new/2, N_new, endpoint=False).reshape(1, 1, N_new), axes=(-1))\n",
+ " x = y.reshape(1, N_new, 1)\n",
+ " \n",
+ " # bandlimit helper\n",
+ " cx = lbd * f_x \n",
+ " cy = lbd * f_y \n",
+ " tx = L_new / 2 / z + jnp.abs(lbd * f_x)\n",
+ " ty = L_new / 2 / z + jnp.abs(lbd * f_y)\n",
+ " \n",
+ " # bandlimit filter for precompensation, not smoothened!\n",
+ " W = (cx**2 * (1 + tx**2) / tx**2 + cy**2 <= 1) * (cy**2 * (1 + ty**2) / ty**2 + cx**2 <= 1)\n",
+ " \n",
+ " # calculate kernels\n",
+ " H_AS = jnp.sqrt(0j + 1 - jnp.abs(f_x * lbd)**2 - jnp.abs(f_y * lbd)**2)\n",
+ " H_Fr = 1 - jnp.abs(f_x * lbd)**2 / 2 - jnp.abs(f_y * lbd)**2 / 2\n",
+ " delta_H = W * jnp.exp(1j * k * z * (H_AS - H_Fr))\n",
+ "\n",
+ " # apply precompensation\n",
+ " psi_precomp = jnp.fft.ifft2(jnp.fft.fft2(jnp.fft.ifftshift(psi_p, axes=(-1, -2))) * delta_H)\n",
+ " # output coordinates\n",
+ " dq = lbd * z / L_new\n",
+ " Q = dq * N * pad_factor\n",
+ " \n",
+ " q_y = jnp.fft.ifftshift(jnp.linspace(-Q/2, Q/2, N_new, endpoint=False).reshape(1, 1, N_new), axes=(-1))\n",
+ " q_x = q_y.reshape(1, N_new, 1)\n",
+ " \n",
+ " H_1 = jnp.exp(1j * k / (2 * z) * (x**2 + y**2))\n",
+ "\n",
+ " if skip_final_phase:\n",
+ " psi_p_final = jnp.fft.fftshift(jnp.fft.fft2(H_1 * psi_precomp), axes=(-1,-2))\n",
+ " else:\n",
+ " H_2 = np.exp(1j * k * z) * jnp.exp(1j * k / (2 * z) * (q_x**2 + q_y**2))\n",
+ " psi_p_final = jnp.fft.fftshift(H_2 * jnp.fft.fft2(H_1 * psi_precomp), axes=(-1,-2))\n",
+ " \n",
+ " psi_final = zero_unpad(psi_p_final, psi.shape)\n",
+ " \n",
+ " return psi_final, Q / 2, delta_H"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "2024-05-25 13:52:45.935229: W external/xla/xla/service/gpu/nvptx_compiler.cc:760] The NVIDIA driver's CUDA version is 12.3 which is older than the ptxas CUDA version (12.5.40). Because the driver is older than the ptxas version, XLA is disabling parallel compilation, which may slow down compilation. You should update your NVIDIA driver or use the NVIDIA-provided CUDA forward compatibility packages.\n"
+ ]
+ }
+ ],
+ "source": [
+ "L_box = 128e-6\n",
+ "D_box = L_box / 16\n",
+ "N_box = 512;\n",
+ "lbd = 500e-9\n",
+ "y_box = jnp.linspace(-L_box/2, L_box/2, N_box, endpoint=False).reshape(1,1, N_box)\n",
+ "x_box = y_box.reshape(1, N_box, 1)\n",
+ "\n",
+ "U_box = ((x_box)**2 <= (D_box / 2)**2) * (y_box**2 <= (D_box / 2)**2) *\\\n",
+ " (jnp.exp(1j * 2 * jnp.pi / lbd * y_box * np.sin(20/ 360 * 2 * jnp.pi)))\n",
+ " \n",
+ "M_box = 8; \n",
+ "z_box = M_box / N_box / lbd * L_box**2 * 2 \n",
+ "U_prop, Q, delta_H = scalable_angular_spectrum(U_box, z_box, lbd, L_box)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAGiCAYAAAC/NyLhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAeuElEQVR4nO3df2xc5b3n8c8Yjyeusac2Bk8GBzAlpQUn3sWhaSxuHPLDFYvJIioVUtQbif4Tin3jhW6Fg3Ztur3Yiu6mBaUQAb38UO/WlTYxmy6QxhGJQxRFZO34ZmxQlN4Y4mQ9WNB4xgnOOLG/+0eXIya/YJzBzuN5v6SvhM95fHz8iOStcY5tn5mZAABwRNZ03wAAAKkgXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAAp0xruJ5//nmVlZVp1qxZqqys1LvvvjudtwMAcMC0heuPf/yjGhoa9NRTT+nAgQP6u7/7O91zzz06evTodN0SAMABvun6IbsLFy7UHXfcoRdeeME79t3vflf333+/WlpapuOWAAAOyJ6ODzo2Nqauri49+eSTScdramq0d+/e89YnEgklEgnv7YmJCf31r3/VNddcI5/P97XfLwAgvcxMIyMjCofDyspK7Yt/0xKuTz75ROPj4yopKUk6XlJSomg0et76lpYWPf3001N1ewCAKTIwMKDS0tKU3mdaH84499WSmV3wFVRjY6NisZg3/DsYAMwM+fn5Kb/PtLziKi4u1lVXXXXeq6uhoaHzXoVJUiAQUCAQmKrbAwBMkcn8c8+0vOLKyclRZWWlOjo6ko53dHSoqqpqOm4JAOCIaXnFJUmPP/64fvKTn2jBggVatGiRXnzxRR09elRr1qyZrlsCADhg2sL14IMP6tNPP9Uvf/lLDQ4Oqry8XG+99ZZuvPHG6bolAIADpu37uC5HPB5XMBic7tsAAFymWCymgoKClN6Hn1UIAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAAp6Qcrt27d+u+++5TOByWz+fTG2+8kXTezNTc3KxwOKzc3FwtWbJEfX19SWsSiYTq6+tVXFysvLw8rVy5UseOHbusTwQAkBlSDtepU6dUUVGhjRs3XvD8+vXrtWHDBm3cuFH79+9XKBTSihUrNDIy4q1paGhQe3u72tratGfPHp08eVK1tbUaHx+f/GcCAMgMdhkkWXt7u/f2xMSEhUIha21t9Y6dPn3agsGgbdq0yczMhoeHze/3W1tbm7fm+PHjlpWVZdu2bftKHzcWi5kkhmEYxvGJxWIptyet/8bV39+vaDSqmpoa71ggEFB1dbX27t0rSerq6tKZM2eS1oTDYZWXl3trzpVIJBSPx5MGAJCZ0hquaDQqSSopKUk6XlJS4p2LRqPKyclRYWHhRdecq6WlRcFg0Js5c+ak87YBAA75Wp4q9Pl8SW+b2XnHznWpNY2NjYrFYt4MDAyk7V4BAG5Ja7hCoZAknffKaWhoyHsVFgqFNDY2phMnTlx0zbkCgYAKCgqSBgCQmdIarrKyMoVCIXV0dHjHxsbG1NnZqaqqKklSZWWl/H5/0prBwUH19vZ6awAAuJjsVN/h5MmT+stf/uK93d/fr56eHhUVFemGG25QQ0ODnnnmGc2dO1dz587VM888o2984xv68Y9/LEkKBoP66U9/qieeeELXXHONioqK9POf/1zz5s3T8uXL0/eZAQBmplQfQ9y5c+cFH2lcvXq1mf3tkfimpiYLhUIWCARs8eLFFolEkq4xOjpqdXV1VlRUZLm5uVZbW2tHjx79yvfA4/AMwzAzYybzOLzPzEyOicfjCgaD030bAIDLFIvFUn5ugZ9VCABwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKekFK6Wlhbdeeedys/P13XXXaf7779fhw4dSlpjZmpublY4HFZubq6WLFmivr6+pDWJREL19fUqLi5WXl6eVq5cqWPHjl3+ZwMAmPFSCldnZ6cee+wx7du3Tx0dHTp79qxqamp06tQpb8369eu1YcMGbdy4Ufv371coFNKKFSs0MjLirWloaFB7e7va2tq0Z88enTx5UrW1tRofH0/fZwYAmJnsMgwNDZkk6+zsNDOziYkJC4VC1tra6q05ffq0BYNB27Rpk5mZDQ8Pm9/vt7a2Nm/N8ePHLSsry7Zt2/aVPm4sFjNJDMMwjOMTi8VSbs9l/RtXLBaTJBUVFUmS+vv7FY1GVVNT460JBAKqrq7W3r17JUldXV06c+ZM0ppwOKzy8nJvzbkSiYTi8XjSAAAy06TDZWZ6/PHHddddd6m8vFySFI1GJUklJSVJa0tKSrxz0WhUOTk5KiwsvOiac7W0tCgYDHozZ86cyd42AMBxkw5XXV2dDh48qD/84Q/nnfP5fElvm9l5x851qTWNjY2KxWLeDAwMTPa2AQCOm1S46uvrtXXrVu3cuVOlpaXe8VAoJEnnvXIaGhryXoWFQiGNjY3pxIkTF11zrkAgoIKCgqQBAGSmlMJlZqqrq9OWLVv0zjvvqKysLOl8WVmZQqGQOjo6vGNjY2Pq7OxUVVWVJKmyslJ+vz9pzeDgoHp7e701AABcVCpPcjz66KMWDAZt165dNjg46M1nn33mrWltbbVgMGhbtmyxSCRiq1atstmzZ1s8HvfWrFmzxkpLS23Hjh3W3d1tS5cutYqKCjt79uxXug+eKmQYhpkZM5mnClMK18U+8CuvvOKtmZiYsKamJguFQhYIBGzx4sUWiUSSrjM6Omp1dXVWVFRkubm5Vltba0ePHv3K90G4GIZhZsZMJly+/x8kp8TjcQWDwem+DQDAZYrFYik/t8DPKgQAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCnZ030DQCYpevE/66pwcdqvO/q/9ujkS39K+3WBKxHhAqbQv3/tP+rqvLK0X/ffPhpVrwgXMgPhAqbQP/xls27Ivibt1/2Xk++pN+1XBa5MhAuYQp/91zadDKf/j13iTyPSP6f9ssAViXABU+gfK5YoUPbNtF/343/9V0l7035d4EpEuIApdPLXjyhRcFPar5voe1WEC5mCcAFT6L+/+4FuuGo47df9H6eO6ddpvypwZSJcwBS67icbFC72p/26wT3D0ptpvyxwRSJcwBT6Tyu/qdyyWWm/7vHTY4QLGYNwAVPoo23/pOxrb0r7dU/1vCLpH9N+XeBKRLiAKfRfXvEpnO1L+3X/94hPr6b9qsCViXABU+jfzft7lV2d/j92kb/EpUjaLwtckQgXMIX+ucn0zfT/xCcdeEmECxmDcAFT6O2jv1G2zUn7dU9++kdJv0n7dYErEeECptDKtYUqzE7/T4f/P6e+oXfSflXgykS4gCn099f+k24OXJ326740dFjvxNN+WeCKRLiAKbTnN+/r8E3p/wbkQ6/GpGfSflngikS4gCm04dQjyh65Lu3XPZV4R9LmtF8XuBIRLmAKzfvl3crLTf9jhUf+78fqI1zIEIQLmELNn/1RZeNFab/uy2PvqS/tVwWuTIQLmEL/9g//U4mv4RdJfvynU/wiSWQMwgVMof92y39Qzk2Fab/uX9/rlrQ77dcFrkSEC5hCQxsKlHV1MO3XTRzOTfs1gSuVz8xsum8iVfF4XMFg+v/wAwCmViwWU0FBQUrvk/U13QsAAF8LwgUAcArhAgA4hXABAJxCuAAATiFcAACnEC4AgFMIFwDAKYQLAOAUwgUAcArhAgA4hXABAJxCuAAATiFcAACnEC4AgFMIFwDAKSmF64UXXtD8+fNVUFCggoICLVq0SG+//bZ33szU3NyscDis3NxcLVmyRH19fUnXSCQSqq+vV3FxsfLy8rRy5UodO3YsPZ8NAGDmsxRs3brV3nzzTTt06JAdOnTI1q1bZ36/33p7e83MrLW11fLz823z5s0WiUTswQcftNmzZ1s8HveusWbNGrv++uuto6PDuru77e6777aKigo7e/bsV76PWCxmkhiGYRjHJxaLpZIhMzNLKVwXUlhYaC+//LJNTExYKBSy1tZW79zp06ctGAzapk2bzMxseHjY/H6/tbW1eWuOHz9uWVlZtm3btq/8MQkXwzDMzJjJhGvS/8Y1Pj6utrY2nTp1SosWLVJ/f7+i0ahqamq8NYFAQNXV1dq7d68kqaurS2fOnElaEw6HVV5e7q25kEQioXg8njQAgMyUcrgikYiuvvpqBQIBrVmzRu3t7brtttsUjUYlSSUlJUnrS0pKvHPRaFQ5OTkqLCy86JoLaWlpUTAY9GbOnDmp3jYAYIZIOVy33nqrenp6tG/fPj366KNavXq13n//fe+8z+dLWm9m5x0715etaWxsVCwW82ZgYCDV2wYAzBAphysnJ0e33HKLFixYoJaWFlVUVOjZZ59VKBSSpPNeOQ0NDXmvwkKhkMbGxnTixImLrrmQQCDgPcn4+QAAMtNlfx+XmSmRSKisrEyhUEgdHR3eubGxMXV2dqqqqkqSVFlZKb/fn7RmcHBQvb293hoAAC4plSc5Ghsbbffu3dbf328HDx60devWWVZWlm3fvt3M/vY4fDAYtC1btlgkErFVq1Zd8HH40tJS27Fjh3V3d9vSpUt5HJ5hGCZD52t/HP6RRx6xG2+80XJycuzaa6+1ZcuWedEyM5uYmLCmpiYLhUIWCARs8eLFFolEkq4xOjpqdXV1VlRUZLm5uVZbW2tHjx5N6aYJF8MwzMyYyYTLZ2Ymx8TjcQWDwem+DQDAZYrFYik/t8DPKgQAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTCBcAwCmECwDgFMIFAHAK4QIAOIVwAQCcQrgAAE4hXAAApxAuAIBTLitcLS0t8vl8amho8I6ZmZqbmxUOh5Wbm6slS5aor68v6f0SiYTq6+tVXFysvLw8rVy5UseOHbucWwEAZIhJh2v//v168cUXNX/+/KTj69ev14YNG7Rx40bt379foVBIK1as0MjIiLemoaFB7e3tamtr0549e3Ty5EnV1tZqfHx88p8JACAz2CSMjIzY3LlzraOjw6qrq23t2rVmZjYxMWGhUMhaW1u9tadPn7ZgMGibNm0yM7Ph4WHz+/3W1tbmrTl+/LhlZWXZtm3bvtLHj8ViJolhGIZxfGKxWMoNmtQrrscee0z33nuvli9fnnS8v79f0WhUNTU13rFAIKDq6mrt3btXktTV1aUzZ84krQmHwyovL/fWnCuRSCgejycNACAzZaf6Dm1tberu7tb+/fvPOxeNRiVJJSUlScdLSkr00UcfeWtycnJUWFh43prP3/9cLS0tevrpp1O9VQDADJTSK66BgQGtXbtWv//97zVr1qyLrvP5fElvm9l5x851qTWNjY2KxWLeDAwMpHLbAIAZJKVwdXV1aWhoSJWVlcrOzlZ2drY6Ozv13HPPKTs723ulde4rp6GhIe9cKBTS2NiYTpw4cdE15woEAiooKEgaAEBmSilcy5YtUyQSUU9PjzcLFizQww8/rJ6eHt18880KhULq6Ojw3mdsbEydnZ2qqqqSJFVWVsrv9yetGRwcVG9vr7cGAICLSvlxjnN88alCM7PW1lYLBoO2ZcsWi0QitmrVKps9e7bF43FvzZo1a6y0tNR27Nhh3d3dtnTpUquoqLCzZ89+pY/JU4UMwzAzYybzVGHKD2d8mV/84hcaHR3Vz372M504cUILFy7U9u3blZ+f76359a9/rezsbP3oRz/S6Oioli1bpldffVVXXXVVum8HADDD+MzMpvsmUhWPxxUMBqf7NgAAlykWi6X83AI/qxAA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOIVwAAKcQLgCAUwgXAMAphAsA4BTCBQBwCuECADiFcAEAnEK4AABOSSlczc3N8vl8SRMKhbzzZqbm5maFw2Hl5uZqyZIl6uvrS7pGIpFQfX29iouLlZeXp5UrV+rYsWPp+WwAADNeyq+4br/9dg0ODnoTiUS8c+vXr9eGDRu0ceNG7d+/X6FQSCtWrNDIyIi3pqGhQe3t7Wpra9OePXt08uRJ1dbWanx8PD2fEQBgZrMUNDU1WUVFxQXPTUxMWCgUstbWVu/Y6dOnLRgM2qZNm8zMbHh42Px+v7W1tXlrjh8/bllZWbZt27avfB+xWMwkMQzDMI5PLBZLJUNmZpbyK67Dhw8rHA6rrKxMDz30kI4cOSJJ6u/vVzQaVU1Njbc2EAiourpae/fulSR1dXXpzJkzSWvC4bDKy8u9NReSSCQUj8eTBgCQmVIK18KFC/X666/rz3/+s1566SVFo1FVVVXp008/VTQalSSVlJQkvU9JSYl3LhqNKicnR4WFhRddcyEtLS0KBoPezJkzJ5XbBgDMICmF65577tEPf/hDzZs3T8uXL9ebb74pSXrttde8NT6fL+l9zOy8Y+f6sjWNjY2KxWLeDAwMpHLbAIAZ5LIeh8/Ly9O8efN0+PBh7+nCc185DQ0Nea/CQqGQxsbGdOLEiYuuuZBAIKCCgoKkAQBkpssKVyKR0AcffKDZs2errKxMoVBIHR0d3vmxsTF1dnaqqqpKklRZWSm/35+0ZnBwUL29vd4aAAAuKZUnOZ544gnbtWuXHTlyxPbt22e1tbWWn59vH374oZmZtba2WjAYtC1btlgkErFVq1bZ7NmzLR6Pe9dYs2aNlZaW2o4dO6y7u9uWLl1qFRUVdvbs2a98HzxVyDAMMzNmMk8VphSuBx980GbPnm1+v9/C4bA98MAD1tfX552fmJiwpqYmC4VCFggEbPHixRaJRJKuMTo6anV1dVZUVGS5ublWW1trR48eTemmCRfDMMzMmMmEy2dmJsfE43EFg8Hpvg0AwGWKxWIpP7fg5M8qdLC1AIALmMzf506G64s/QgoA4K7J/H3u5JcKJyYmdOjQId12220aGBjg8fgLiMfjmjNnDvtzEezPpbE/l8b+fLkv2yMz08jIiMLhsLKyUnsNlZ2um5xKWVlZuv766yWJ7+v6EuzPpbE/l8b+XBr78+UutUeTfVbByS8VAgAyF+ECADjF2XAFAgE1NTUpEAhM961ckdifS2N/Lo39uTT258t9nXvk5MMZAIDM5ewrLgBAZiJcAACnEC4AgFMIFwDAKU6G6/nnn1dZWZlmzZqlyspKvfvuu9N9S1Ni9+7duu+++xQOh+Xz+fTGG28knTczNTc3KxwOKzc3V0uWLFFfX1/SmkQiofr6ehUXFysvL08rV67UsWPHpvCz+Pq0tLTozjvvVH5+vq677jrdf//9OnToUNKaTN6jF154QfPnz/e+IXTRokV6++23vfOZvDcX0tLSIp/Pp4aGBu9YJu9Rc3OzfD5f0nz+C4SlKd6blH+e/DRra2szv99vL730kr3//vu2du1ay8vLs48++mi6b+1r99Zbb9lTTz1lmzdvNknW3t6edL61tdXy8/Nt8+bNFolEvF9Dc+7vQ7v++uuto6PDuru77e67707596FdqX7wgx/YK6+8Yr29vdbT02P33nuv3XDDDXby5ElvTSbv0datW+3NN9+0Q4cO2aFDh2zdunXm9/utt7fXzDJ7b8713nvv2U033WTz58+3tWvXesczeY+amprs9ttvt8HBQW+Ghoa881O5N86F63vf+56tWbMm6dh3vvMde/LJJ6fpjqbHueGamJiwUChkra2t3rHTp09bMBi0TZs2mZnZ8PCw+f1+a2tr89YcP37csrKybNu2bVN271NlaGjIJFlnZ6eZsUcXUlhYaC+//DJ78wUjIyM2d+5c6+josOrqai9cmb5HTU1NVlFRccFzU703Tn2pcGxsTF1dXaqpqUk6XlNTo717907TXV0Z+vv7FY1Gk/YmEAiourra25uuri6dOXMmaU04HFZ5efmM3L9YLCZJKioqksQefdH4+Lja2tp06tQpLVq0iL35gscee0z33nuvli9fnnScPZIOHz6scDissrIyPfTQQzpy5Iikqd8bp37I7ieffKLx8XGVlJQkHS8pKVE0Gp2mu7oyfP75X2hvPvroI29NTk6OCgsLz1sz0/bPzPT444/rrrvuUnl5uST2SJIikYgWLVqk06dP6+qrr1Z7e7tuu+027y+OTN4bSWpra1N3d7f2799/3rlM//9n4cKFev311/Xtb39bH3/8sX71q1+pqqpKfX19U743ToXrcz6fL+ltMzvvWKaazN7MxP2rq6vTwYMHtWfPnvPOZfIe3Xrrrerp6dHw8LA2b96s1atXq7Oz0zufyXszMDCgtWvXavv27Zo1a9ZF12XqHt1zzz3ef8+bN0+LFi3St771Lb322mv6/ve/L2nq9sapLxUWFxfrqquuOq/OQ0ND55U+03z+dM+l9iYUCmlsbEwnTpy46JqZoL6+Xlu3btXOnTtVWlrqHWePpJycHN1yyy1asGCBWlpaVFFRoWeffZa90d++lDU0NKTKykplZ2crOztbnZ2deu6555Sdne19jpm8R1+Ul5enefPm6fDhw1P+/49T4crJyVFlZaU6OjqSjnd0dKiqqmqa7urKUFZWplAolLQ3Y2Nj6uzs9PamsrJSfr8/ac3g4KB6e3tnxP6Zmerq6rRlyxa98847KisrSzrPHp3PzJRIJNgbScuWLVMkElFPT483CxYs0MMPP6yenh7dfPPNGb9HX5RIJPTBBx9o9uzZU///T0qPclwBPn8c/ne/+529//771tDQYHl5efbhhx9O96197UZGRuzAgQN24MABk2QbNmywAwcOeN8K0NraasFg0LZs2WKRSMRWrVp1wcdRS0tLbceOHdbd3W1Lly6dEY/qmpk9+uijFgwGbdeuXUmP7H722Wfemkzeo8bGRtu9e7f19/fbwYMHbd26dZaVlWXbt283s8zem4v54lOFZpm9R0888YTt2rXLjhw5Yvv27bPa2lrLz8/3/u6dyr1xLlxmZr/97W/txhtvtJycHLvjjju8x51nup07d5qk82b16tVm9rdHUpuamiwUClkgELDFixdbJBJJusbo6KjV1dVZUVGR5ebmWm1trR09enQaPpv0u9DeSLJXXnnFW5PJe/TII494f26uvfZaW7ZsmRcts8zem4s5N1yZvEeff1+W3++3cDhsDzzwgPX19Xnnp3Jv+LUmAACnOPVvXAAAEC4AgFMIFwDAKYQLAOAUwgUAcArhAgA4hXABAJxCuAAATiFcAACnEC4AgFMIFwDAKYQLAOCU/wetvJp/jAYa2gAAAABJRU5ErkJggg==",
+ "text/plain": [
+ "