From 20fd09ed198f57f97bc2cf0bda595c73ce602430 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Mon, 9 Jun 2025 17:54:52 -0700 Subject: [PATCH 01/25] propagate gradients through `isotropic_thin_3d` --- waveorder/models/isotropic_thin_3d.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/waveorder/models/isotropic_thin_3d.py b/waveorder/models/isotropic_thin_3d.py index 85a6a21a..4fd21ae4 100644 --- a/waveorder/models/isotropic_thin_3d.py +++ b/waveorder/models/isotropic_thin_3d.py @@ -1,6 +1,5 @@ from typing import Literal, Tuple -import numpy as np import torch from torch import Tensor @@ -38,7 +37,7 @@ def generate_test_phantom( def calculate_transfer_function( yx_shape: Tuple[int, int], yx_pixel_size: float, - z_position_list: list, + z_position_list: torch.Tensor, wavelength_illumination: float, index_of_refraction_media: float, numerical_aperture_illumination: float, @@ -50,7 +49,7 @@ def calculate_transfer_function( numerical_aperture_illumination, numerical_aperture_detection, ) - yx_factor = int(np.ceil(yx_pixel_size / transverse_nyquist)) + yx_factor = int(torch.ceil(yx_pixel_size / transverse_nyquist)) ( absorption_2d_to_3d_transfer_function, @@ -97,7 +96,7 @@ def calculate_transfer_function( def _calculate_wrap_unsafe_transfer_function( yx_shape: Tuple[int, int], yx_pixel_size: float, - z_position_list: list, + z_position_list: torch.Tensor, wavelength_illumination: float, index_of_refraction_media: float, numerical_aperture_illumination: float, @@ -114,7 +113,7 @@ def _calculate_wrap_unsafe_transfer_function( numerical_aperture_illumination = 0.9 * numerical_aperture_detection if invert_phase_contrast: - z_position_list = [-1 * x for x in z_position_list] + z_position_list *= -1 radial_frequencies = util.generate_radial_frequencies( yx_shape, yx_pixel_size ) @@ -133,7 +132,7 @@ def _calculate_wrap_unsafe_transfer_function( radial_frequencies, detection_pupil, wavelength_illumination / index_of_refraction_media, - torch.tensor(z_position_list), + z_position_list, ) zyx_shape = (len(z_position_list),) + tuple(yx_shape) From 6acc8b75667dd2e8018baef2eaa690694655b345 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 10 Jun 2025 12:10:51 -0700 Subject: [PATCH 02/25] pupil rolloff --- tests/test_optics.py | 18 +++++++++++++++--- waveorder/optics.py | 37 ++++++++++++++++++++----------------- 2 files changed, 35 insertions(+), 20 deletions(-) diff --git a/tests/test_optics.py b/tests/test_optics.py index d51b5bef..6090889f 100644 --- a/tests/test_optics.py +++ b/tests/test_optics.py @@ -8,11 +8,23 @@ def test_generate_pupil(): pupil = optics.generate_pupil(radial_frequencies, 0.5, 0.5) # Corners are in the pupil - assert pupil[0, 0] == 1 - assert pupil[-1, -1] == 1 + assert torch.isclose(pupil[0, 0], torch.tensor(1.0), rtol=1e-3) + assert torch.isclose(pupil[-1, -1], torch.tensor(1.0), rtol=1e-3) # Center is outside the pupil - assert pupil[5, 5] == 0 + assert pupil[5, 5] < 1e-3 + + +def test_generate_pupil_cutoff(): + """ + Test generate_pupil at the cutoff frequency. + """ + frr = torch.tensor([[0.5, 1.0, 1.5]]) + NA = 1.0 + lamb_in = 1.0 + pupil = optics.generate_pupil(frr, NA, lamb_in) + # At cutoff, sigmoid should be ~0.5 + assert torch.isclose(pupil[0, 1], torch.tensor(0.5), atol=1e-3) def test_generate_propagation_kernel(): diff --git a/waveorder/optics.py b/waveorder/optics.py index 46204517..c28cf344 100644 --- a/waveorder/optics.py +++ b/waveorder/optics.py @@ -118,34 +118,37 @@ def analyzer_output(Ein, alpha, beta): return Eout -def generate_pupil(frr, NA, lamb_in): +def generate_pupil( + frr: torch.Tensor, NA: float, lamb_in: float, slope: float = 4.0 +) -> torch.Tensor: """ - - compute pupil function given spatial frequency, NA, wavelength. + Generate a soft-edged pupil function using a sigmoid roll-off. The default + sigmoid softens the edge within ~1 voxel of the boundary. Parameters ---------- - frr : torch.tensor - radial frequency coordinate in units of inverse length + frr : torch.Tensor + radial frequency coordinates (units: 1/length) - NA : float - numerical aperture of the pupil function (normalized by the refractive index of the immersion media) + NA : float + numerical aperture (unitless) lamb_in : float - wavelength of the light in free space - in units of length (inverse of frr's units) + wavelength of light (same length units as frr^-1) + + slope : float, optional + steepness of the sigmoid roll-off (default 4.0 keeps ~90% of + sigmoid within a single voxel) Returns ------- - Pupil : numpy.ndarray - pupil function with the specified parameters with the size of (Ny, Nx) - + pupil : torch.Tensor + pupil function, pupil.shape == frr.shape, values in [0, 1] """ - - Pupil = torch.zeros(frr.shape) - Pupil[frr < NA / lamb_in] = 1 - - return Pupil + pixel_slope = slope / torch.abs(frr[0, 1] - frr[0, 0]) + cutoff = NA / lamb_in + pupil = torch.sigmoid(pixel_slope * (cutoff - frr)) + return pupil def gen_sector_Pupil(fxx, fyy, NA, lamb_in, sector_angle, rotation_angle): From 2c9e5b8202bd8c6bb14795b42062971a2fcf74fc Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 10 Jun 2025 12:31:19 -0700 Subject: [PATCH 03/25] fix bug --- waveorder/models/isotropic_thin_3d.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/waveorder/models/isotropic_thin_3d.py b/waveorder/models/isotropic_thin_3d.py index 4fd21ae4..d36487ee 100644 --- a/waveorder/models/isotropic_thin_3d.py +++ b/waveorder/models/isotropic_thin_3d.py @@ -2,6 +2,7 @@ import torch from torch import Tensor +import numpy as np from waveorder import optics, sampling, util from waveorder.filter import apply_filter_bank @@ -49,7 +50,7 @@ def calculate_transfer_function( numerical_aperture_illumination, numerical_aperture_detection, ) - yx_factor = int(torch.ceil(yx_pixel_size / transverse_nyquist)) + yx_factor = int(np.ceil(yx_pixel_size / transverse_nyquist)) ( absorption_2d_to_3d_transfer_function, From aa89bbd4063a4f8a6ee9b2567d4de178d4dd966c Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 10 Jun 2025 12:32:16 -0700 Subject: [PATCH 04/25] change z_position_list to a tensor everywhere --- docs/examples/models/isotropic_thin_3d.py | 3 ++- tests/cli_tests/test_compute_tf.py | 4 ++-- waveorder/cli/compute_transfer_function.py | 14 ++++++++------ waveorder/models/isotropic_thin_3d.py | 2 +- 4 files changed, 13 insertions(+), 10 deletions(-) diff --git a/docs/examples/models/isotropic_thin_3d.py b/docs/examples/models/isotropic_thin_3d.py index 577c4064..3a560933 100644 --- a/docs/examples/models/isotropic_thin_3d.py +++ b/docs/examples/models/isotropic_thin_3d.py @@ -5,6 +5,7 @@ import napari import numpy as np +import torch from waveorder.models import isotropic_thin_3d @@ -27,7 +28,7 @@ ] ) transfer_function_arguments = { - "z_position_list": (np.arange(z_shape) - z_shape // 2) * z_pixel_size, + "z_position_list": (torch.arange(z_shape) - z_shape // 2) * z_pixel_size, "numerical_aperture_illumination": 0.9, "numerical_aperture_detection": 1.2, } diff --git a/tests/cli_tests/test_compute_tf.py b/tests/cli_tests/test_compute_tf.py index a31d5b9e..ee1b155c 100644 --- a/tests/cli_tests/test_compute_tf.py +++ b/tests/cli_tests/test_compute_tf.py @@ -1,5 +1,5 @@ -import numpy as np import pytest +import torch from click.testing import CliRunner from waveorder.cli import settings @@ -22,7 +22,7 @@ ) def test_position_list_from_shape_scale_offset(shape, scale, offset, expected): result = _position_list_from_shape_scale_offset(shape, scale, offset) - np.testing.assert_allclose(result, expected) + torch.testing.assert_close(result, torch.tensor(expected)) def test_compute_transfer(tmp_path, example_plate): diff --git a/waveorder/cli/compute_transfer_function.py b/waveorder/cli/compute_transfer_function.py index 24497005..86c04fb2 100644 --- a/waveorder/cli/compute_transfer_function.py +++ b/waveorder/cli/compute_transfer_function.py @@ -1,7 +1,7 @@ from pathlib import Path import click -import numpy as np +import torch from iohub.ngff import Position, open_ome_zarr from waveorder import focus @@ -23,18 +23,20 @@ def _position_list_from_shape_scale_offset( shape: int, scale: float, offset: float -) -> list: +) -> torch.Tensor: """ - Generates a list of positions based on the given array shape, pixel size (scale), and offset. + Generates a 1D tensor of positions based on the given array shape, pixel size (scale), and offset. Examples -------- >>> _position_list_from_shape_scale_offset(5, 1.0, 0.0) - [2.0, 1.0, 0.0, -1.0, -2.0] + tensor([ 2., 1., 0., -1., -2.]) >>> _position_list_from_shape_scale_offset(4, 0.5, 1.0) - [1.5, 1.0, 0.5, 0.0] + tensor([1.5, 1.0, 0.5, 0.0]) """ - return list((-np.arange(shape) + (shape // 2) + offset) * scale) + return ( + -torch.arange(shape, dtype=torch.float32) + (shape // 2) + offset + ) * scale def generate_and_save_birefringence_transfer_function(settings, dataset): diff --git a/waveorder/models/isotropic_thin_3d.py b/waveorder/models/isotropic_thin_3d.py index d36487ee..f64bb2aa 100644 --- a/waveorder/models/isotropic_thin_3d.py +++ b/waveorder/models/isotropic_thin_3d.py @@ -1,8 +1,8 @@ from typing import Literal, Tuple +import numpy as np import torch from torch import Tensor -import numpy as np from waveorder import optics, sampling, util from waveorder.filter import apply_filter_bank From be06cbdb7865284c235d8fe3a03157c0bea82d58 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 10 Jun 2025 13:23:06 -0700 Subject: [PATCH 05/25] fix tf test --- tests/models/test_isotropic_thin_3d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/models/test_isotropic_thin_3d.py b/tests/models/test_isotropic_thin_3d.py index 9bd49bf7..09f66675 100644 --- a/tests/models/test_isotropic_thin_3d.py +++ b/tests/models/test_isotropic_thin_3d.py @@ -9,7 +9,7 @@ def test_calculate_transfer_function(invert_phase_contrast): Hu, Hp = isotropic_thin_3d.calculate_transfer_function( yx_shape=(100, 101), yx_pixel_size=6.5 / 40, - z_position_list=[-1, 0, 1], + z_position_list=torch.tensor([-1, 0, 1]), wavelength_illumination=0.5, index_of_refraction_media=1.0, numerical_aperture_illumination=0.4, From 9f9e1cd0e7e73534afbd5ab81954213f02c592f4 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 10 Jun 2025 13:48:56 -0700 Subject: [PATCH 06/25] loosen optics tests to accomodate rolloff --- tests/test_optics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_optics.py b/tests/test_optics.py index 6090889f..590b2dee 100644 --- a/tests/test_optics.py +++ b/tests/test_optics.py @@ -39,7 +39,7 @@ def test_generate_propagation_kernel(): assert propagation_kernel.shape == (3, 10, 10) assert propagation_kernel[1, 0, 0] == 1 - assert propagation_kernel[1, 5, 5] == 0 + assert torch.abs(propagation_kernel[1, 5, 5]) < 1e-3 def test_gen_Greens_function_z(): @@ -53,7 +53,7 @@ def test_gen_Greens_function_z(): ) assert G.shape == (3, 10, 10) - assert G[1, 5, 5] == 0 + assert torch.abs(G[1, 5, 5]) < 1e-3 def test_WOTF_2D(): From 3520b38dc0205b71fbbb13e5310ae753f14d50db Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 10 Jun 2025 13:54:42 -0700 Subject: [PATCH 07/25] soften wotf optics test for rolloff --- tests/test_optics.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/test_optics.py b/tests/test_optics.py index 590b2dee..19f6522a 100644 --- a/tests/test_optics.py +++ b/tests/test_optics.py @@ -69,7 +69,9 @@ def test_WOTF_2D(): ) # Absorption DC term - assert absorption_transfer_function[0, 0] == 2 + assert torch.isclose( + absorption_transfer_function[0, 0], torch.tensor(2.0), rtol=1e-3 + ) # No phase contrast for an in-focus slice assert torch.all(torch.real(phase_transfer_function) == 0) From d6c8bceedcdb7cecdc32b270784a0c00dd36ad4d Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 10 Jun 2025 14:11:47 -0700 Subject: [PATCH 08/25] fix sqrt(-1) -> nan bug --- waveorder/optics.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/waveorder/optics.py b/waveorder/optics.py index c28cf344..4d83ecbc 100644 --- a/waveorder/optics.py +++ b/waveorder/optics.py @@ -464,7 +464,9 @@ def generate_greens_function_z( oblique_factor = ( (1 - wavelength_illumination**2 * radial_frequencies**2) - * pupil_support + * pupil_support.type( + torch.complex64 + ) # complex to avoid sqrt(-1) -> nan ) ** (1 / 2) / wavelength_illumination if axially_even: From fd7c79e6fc888e44597b88e136a33f82ea62dfaa Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 10 Jun 2025 14:20:52 -0700 Subject: [PATCH 09/25] fix optics wotf test --- tests/test_optics.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/test_optics.py b/tests/test_optics.py index 19f6522a..d5a00df8 100644 --- a/tests/test_optics.py +++ b/tests/test_optics.py @@ -70,8 +70,11 @@ def test_WOTF_2D(): # Absorption DC term assert torch.isclose( - absorption_transfer_function[0, 0], torch.tensor(2.0), rtol=1e-3 + torch.real(absorption_transfer_function[0, 0]), + torch.tensor(2.0), + rtol=1e-3, ) + assert torch.abs(torch.imag(absorption_transfer_function[0, 0])) < 1e-3 # No phase contrast for an in-focus slice assert torch.all(torch.real(phase_transfer_function) == 0) From 93acbbafb00fe35dfd684437d99697862e1a4602 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 17 Jun 2025 16:11:23 -0700 Subject: [PATCH 10/25] enable tilted pupils --- waveorder/models/isotropic_thin_3d.py | 44 +++++++++++------ waveorder/optics.py | 70 +++++++++++++++++++++++++++ 2 files changed, 99 insertions(+), 15 deletions(-) diff --git a/waveorder/models/isotropic_thin_3d.py b/waveorder/models/isotropic_thin_3d.py index f64bb2aa..69f783ae 100644 --- a/waveorder/models/isotropic_thin_3d.py +++ b/waveorder/models/isotropic_thin_3d.py @@ -1,6 +1,5 @@ from typing import Literal, Tuple -import numpy as np import torch from torch import Tensor @@ -44,13 +43,17 @@ def calculate_transfer_function( numerical_aperture_illumination: float, numerical_aperture_detection: float, invert_phase_contrast: bool = False, + tilt_angle_zenith: float = 0.0, + tilt_angle_azimuth: float = 0.0, ) -> Tuple[Tensor, Tensor]: transverse_nyquist = sampling.transverse_nyquist( wavelength_illumination, numerical_aperture_illumination, numerical_aperture_detection, ) - yx_factor = int(np.ceil(yx_pixel_size / transverse_nyquist)) + yx_factor = int( + torch.ceil(torch.tensor(yx_pixel_size / transverse_nyquist)) + ) ( absorption_2d_to_3d_transfer_function, @@ -67,6 +70,8 @@ def calculate_transfer_function( numerical_aperture_illumination, numerical_aperture_detection, invert_phase_contrast=invert_phase_contrast, + tilt_angle_zenith=tilt_angle_zenith, + tilt_angle_azimuth=tilt_angle_azimuth, ) absorption_2d_to_3d_transfer_function_out = torch.zeros( @@ -103,6 +108,8 @@ def _calculate_wrap_unsafe_transfer_function( numerical_aperture_illumination: float, numerical_aperture_detection: float, invert_phase_contrast: bool = False, + tilt_angle_zenith: float = 0.0, + tilt_angle_azimuth: float = 0.0, ) -> Tuple[Tensor, Tensor]: if numerical_aperture_illumination >= numerical_aperture_detection: print( @@ -114,15 +121,21 @@ def _calculate_wrap_unsafe_transfer_function( numerical_aperture_illumination = 0.9 * numerical_aperture_detection if invert_phase_contrast: - z_position_list *= -1 - radial_frequencies = util.generate_radial_frequencies( - yx_shape, yx_pixel_size - ) + z_positions = z_position_list * -1 + else: + z_positions = z_position_list.clone() - illumination_pupil = optics.generate_pupil( - radial_frequencies, + fyy, fxx = util.generate_frequencies(yx_shape, yx_pixel_size) + radial_frequencies = torch.sqrt(fyy**2 + fxx**2) + + illumination_pupil = optics.generate_tilted_pupil( + fxx, + fyy, numerical_aperture_illumination, wavelength_illumination, + index_of_refraction_media, + tilt_angle_zenith, + tilt_angle_azimuth, ) detection_pupil = optics.generate_pupil( radial_frequencies, @@ -133,17 +146,17 @@ def _calculate_wrap_unsafe_transfer_function( radial_frequencies, detection_pupil, wavelength_illumination / index_of_refraction_media, - z_position_list, + z_positions, ) - zyx_shape = (len(z_position_list),) + tuple(yx_shape) + zyx_shape = (len(z_positions),) + tuple(yx_shape) absorption_2d_to_3d_transfer_function = torch.zeros( zyx_shape, dtype=torch.complex64 ) phase_2d_to_3d_transfer_function = torch.zeros( zyx_shape, dtype=torch.complex64 ) - for z in range(len(z_position_list)): + for z in range(len(z_positions)): ( absorption_2d_to_3d_transfer_function[z], phase_2d_to_3d_transfer_function[z], @@ -257,17 +270,18 @@ def apply_transfer_function( # simulate absorbing object yx_absorption_hat = torch.fft.fftn(yx_absorption) - zyx_absorption_data_hat = yx_absorption_hat[None, ...] * torch.real( - absorption_2d_to_3d_transfer_function + zyx_absorption_data_hat = ( + yx_absorption_hat[None, ...] * absorption_2d_to_3d_transfer_function ) + zyx_absorption_data = torch.real( torch.fft.ifftn(zyx_absorption_data_hat, dim=(1, 2)) ) # simulate phase object yx_phase_hat = torch.fft.fftn(yx_phase) - zyx_phase_data_hat = yx_phase_hat[None, ...] * torch.real( - phase_2d_to_3d_transfer_function + zyx_phase_data_hat = ( + yx_phase_hat[None, ...] * phase_2d_to_3d_transfer_function ) zyx_phase_data = torch.real( torch.fft.ifftn(zyx_phase_data_hat, dim=(1, 2)) diff --git a/waveorder/optics.py b/waveorder/optics.py index 4d83ecbc..a5f046e8 100644 --- a/waveorder/optics.py +++ b/waveorder/optics.py @@ -1,4 +1,5 @@ import itertools +import math import numpy as np import torch @@ -151,6 +152,75 @@ def generate_pupil( return pupil +def generate_tilted_pupil( + fxx: torch.Tensor, + fyy: torch.Tensor, + NA: float, + lamb_in: float, + n: float = 1.0, + tilt_angle_zenith: float = 0.0, + tilt_angle_azimuth: float = 0.0, + slope: float = 4.0, +) -> torch.Tensor: + """ + Generate a soft-edged 2-D pupil that may be tilted on the Ewald sphere. + + Parameters + ---------- + fxx, fyy : torch.Tensor + Cartesian frequency grids (units: 1/length) with identical shape. + NA : float + Numerical aperture of the objective (must satisfy NA ≤ n). + lamb_in : float + Illumination wavelength (units: length) + n : float, optional + Refractive index of the immersion medium (default 1.0). + tilt_angle_zenith : float, optional + Polar angle θ (radians) between the pupil axis and +z (0 = untilted). + tilt_angle_azimuth : float, optional + Azimuth φ (radians) of the tilt in the xy-plane (0 = +x). + slope : float, optional + Controls sigmoid roll-off (≈ 90 % change in one pixel when slope=4). + + Returns + ------- + pupil : torch.Tensor + 2-D soft mask with values in [0, 1] and shape == fx.shape. + """ + if NA > n: + raise ValueError("NA must be ≤ n (otherwise angle would be complex).") + + # constants + K = n / lamb_in # Ewald-sphere radius + cos_alpha_max = math.sqrt(1 - (NA / n) ** 2) + + # sampling metrics + # Assume fxx, fyy are on a regular grid ⇒ pixel spacing in fx direction: + df = torch.abs(fxx[0, 1] - fxx[0, 0]) + pixel_slope = slope / df + + # sphere coordinates + fz_sq = K**2 - fxx**2 - fyy**2 + inside_sphere = fz_sq >= 0 + # clamp to avoid negative round-off, but keep gradients + fz = torch.sqrt(torch.clamp(fz_sq, min=0.0)) + + # tilt unit vector + sx = torch.sin(tilt_angle_zenith) * torch.cos(tilt_angle_azimuth) + sy = torch.sin(tilt_angle_zenith) * torch.sin(tilt_angle_azimuth) + sz = torch.cos(tilt_angle_zenith) + + # dot-product test + dot = fxx * sx + fyy * sy + fz * sz # v · s + threshold = K * cos_alpha_max + pupil_soft = torch.sigmoid(pixel_slope * (dot - threshold)) + + # mask outside sphere + pupil = pupil_soft * inside_sphere.to(fxx.dtype) + + return pupil + + def gen_sector_Pupil(fxx, fyy, NA, lamb_in, sector_angle, rotation_angle): """ From 4e22a57ae569afc5c9c9bcff7612c3de39126680 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Thu, 26 Jun 2025 16:55:37 -0700 Subject: [PATCH 11/25] wip phase-only recon --- waveorder/models/isotropic_thin_3d.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/waveorder/models/isotropic_thin_3d.py b/waveorder/models/isotropic_thin_3d.py index 69f783ae..09a9c732 100644 --- a/waveorder/models/isotropic_thin_3d.py +++ b/waveorder/models/isotropic_thin_3d.py @@ -201,11 +201,21 @@ def calculate_singular_system( ), dim=0, ) + # XXX TEMP PHASE ONLY + sfYX_transfer_function = phase_2d_to_3d_transfer_function.unsqueeze(0) + U = torch.ones_like(sfYX_transfer_function) + S = torch.linalg.norm(sfYX_transfer_function, dim=1) + Vh = sfYX_transfer_function / S.unsqueeze(0) + return U, S, Vh + # XXX END TEMP + YXsf_transfer_function = sfYX_transfer_function.permute(2, 3, 0, 1) Up, Sp, Vhp = torch.linalg.svd(YXsf_transfer_function, full_matrices=False) + U = Up.permute(2, 3, 0, 1) S = Sp.permute(2, 0, 1) Vh = Vhp.permute(2, 3, 0, 1) + return U, S, Vh @@ -349,10 +359,12 @@ def apply_inverse_transfer_function( U, S, Vh = singular_system S_reg = S / (S**2 + regularization_strength) sfyx_inverse_filter = torch.einsum( - "sj...,j...,jf...->fs...", U, S_reg, Vh + "sj...,j...,jf...->fs...", U, S_reg, Vh.conj() ) - absorption_yx, phase_yx = apply_filter_bank(sfyx_inverse_filter, zyx) + # absorption_yx, phase_yx = apply_filter_bank(sfyx_inverse_filter, zyx) + phase_yx = apply_filter_bank(sfyx_inverse_filter, zyx)[0] + absorption_yx = torch.zeros_like(phase_yx) # ADMM deconvolution with anisotropic TV regularization elif reconstruction_algorithm == "TV": From 1867a0229bccf76c668445a35b5963dbf8a0e3d8 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 1 Jul 2025 10:54:46 -0700 Subject: [PATCH 12/25] clean up phase-only comments --- waveorder/models/isotropic_thin_3d.py | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/waveorder/models/isotropic_thin_3d.py b/waveorder/models/isotropic_thin_3d.py index 09a9c732..5d93c325 100644 --- a/waveorder/models/isotropic_thin_3d.py +++ b/waveorder/models/isotropic_thin_3d.py @@ -201,22 +201,26 @@ def calculate_singular_system( ), dim=0, ) - # XXX TEMP PHASE ONLY + + # phase only reconstruction sfYX_transfer_function = phase_2d_to_3d_transfer_function.unsqueeze(0) U = torch.ones_like(sfYX_transfer_function) S = torch.linalg.norm(sfYX_transfer_function, dim=1) Vh = sfYX_transfer_function / S.unsqueeze(0) return U, S, Vh - # XXX END TEMP - YXsf_transfer_function = sfYX_transfer_function.permute(2, 3, 0, 1) - Up, Sp, Vhp = torch.linalg.svd(YXsf_transfer_function, full_matrices=False) + # Absorption and phase reconstruction + # Gradients do not work with complex-valued SVD, so only phase reconstructions + # are supported for now. - U = Up.permute(2, 3, 0, 1) - S = Sp.permute(2, 0, 1) - Vh = Vhp.permute(2, 3, 0, 1) + # YXsf_transfer_function = sfYX_transfer_function.permute(2, 3, 0, 1) + # Up, Sp, Vhp = torch.linalg.svd(YXsf_transfer_function, full_matrices=False) - return U, S, Vh + # U = Up.permute(2, 3, 0, 1) + # S = Sp.permute(2, 0, 1) + # Vh = Vhp.permute(2, 3, 0, 1) + + # return U, S, Vh def visualize_transfer_function( @@ -362,10 +366,12 @@ def apply_inverse_transfer_function( "sj...,j...,jf...->fs...", U, S_reg, Vh.conj() ) + # Phase only reconstruction # absorption_yx, phase_yx = apply_filter_bank(sfyx_inverse_filter, zyx) - phase_yx = apply_filter_bank(sfyx_inverse_filter, zyx)[0] absorption_yx = torch.zeros_like(phase_yx) + phase_yx = apply_filter_bank(sfyx_inverse_filter, zyx)[0] + # ADMM deconvolution with anisotropic TV regularization elif reconstruction_algorithm == "TV": raise NotImplementedError From e38af9e8e509be61c02fe0628cff76b16e1bc006 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 1 Jul 2025 13:58:10 -0700 Subject: [PATCH 13/25] bugfix --- waveorder/models/isotropic_thin_3d.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/waveorder/models/isotropic_thin_3d.py b/waveorder/models/isotropic_thin_3d.py index 5d93c325..a22804c5 100644 --- a/waveorder/models/isotropic_thin_3d.py +++ b/waveorder/models/isotropic_thin_3d.py @@ -368,9 +368,8 @@ def apply_inverse_transfer_function( # Phase only reconstruction # absorption_yx, phase_yx = apply_filter_bank(sfyx_inverse_filter, zyx) - absorption_yx = torch.zeros_like(phase_yx) - phase_yx = apply_filter_bank(sfyx_inverse_filter, zyx)[0] + absorption_yx = torch.zeros_like(phase_yx) # ADMM deconvolution with anisotropic TV regularization elif reconstruction_algorithm == "TV": From 6a7d4e2134e341b151d99c2dd688d30f55e73fca Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 1 Jul 2025 16:55:24 -0700 Subject: [PATCH 14/25] reattach gradients --- waveorder/optics.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/waveorder/optics.py b/waveorder/optics.py index a5f046e8..845a1d8f 100644 --- a/waveorder/optics.py +++ b/waveorder/optics.py @@ -1,5 +1,4 @@ import itertools -import math import numpy as np import torch @@ -192,7 +191,7 @@ def generate_tilted_pupil( # constants K = n / lamb_in # Ewald-sphere radius - cos_alpha_max = math.sqrt(1 - (NA / n) ** 2) + cos_alpha_max = torch.sqrt(1 - (NA / n) ** 2) # sampling metrics # Assume fxx, fyy are on a regular grid ⇒ pixel spacing in fx direction: From f5420fdf13434422dcc3f896cd35845b38ce4984 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 1 Jul 2025 17:02:23 -0700 Subject: [PATCH 15/25] better derivative through NA_ill > NA_det condition --- waveorder/models/isotropic_thin_3d.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/waveorder/models/isotropic_thin_3d.py b/waveorder/models/isotropic_thin_3d.py index a22804c5..4dd71073 100644 --- a/waveorder/models/isotropic_thin_3d.py +++ b/waveorder/models/isotropic_thin_3d.py @@ -118,7 +118,11 @@ def _calculate_wrap_unsafe_transfer_function( "numerical_aperture_illumination to 0.9 * " "numerical_aperture_detection to avoid singularities." ) - numerical_aperture_illumination = 0.9 * numerical_aperture_detection + numerical_aperture_illumination = torch.where( + numerical_aperture_illumination >= numerical_aperture_detection, + numerical_aperture_detection, + numerical_aperture_illumination, + ) if invert_phase_contrast: z_positions = z_position_list * -1 @@ -137,6 +141,7 @@ def _calculate_wrap_unsafe_transfer_function( tilt_angle_zenith, tilt_angle_azimuth, ) + detection_pupil = optics.generate_pupil( radial_frequencies, numerical_aperture_detection, From de78f61a6bee6bcc84f338ca06dd7009e001aa81 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 1 Jul 2025 17:52:34 -0700 Subject: [PATCH 16/25] quieter --- waveorder/models/isotropic_thin_3d.py | 1 - 1 file changed, 1 deletion(-) diff --git a/waveorder/models/isotropic_thin_3d.py b/waveorder/models/isotropic_thin_3d.py index 4dd71073..f8ff96c0 100644 --- a/waveorder/models/isotropic_thin_3d.py +++ b/waveorder/models/isotropic_thin_3d.py @@ -364,7 +364,6 @@ def apply_inverse_transfer_function( # TODO Consider refactoring with vectorial transfer function SVD if reconstruction_algorithm == "Tikhonov": - print("Computing inverse filter") U, S, Vh = singular_system S_reg = S / (S**2 + regularization_strength) sfyx_inverse_filter = torch.einsum( From f7cf345d3277fc713674f5c5f5438507940455b5 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Tue, 1 Jul 2025 18:13:38 -0700 Subject: [PATCH 17/25] add optimize_phase_recon.py script --- docs/examples/visuals/optimize_phase_recon.py | 271 ++++++++++++++++++ 1 file changed, 271 insertions(+) create mode 100644 docs/examples/visuals/optimize_phase_recon.py diff --git a/docs/examples/visuals/optimize_phase_recon.py b/docs/examples/visuals/optimize_phase_recon.py new file mode 100644 index 00000000..a649d064 --- /dev/null +++ b/docs/examples/visuals/optimize_phase_recon.py @@ -0,0 +1,271 @@ +from datetime import datetime + +import napari +import numpy as np +import torch +from biahub.cli.utils import model_to_yaml +from biahub.settings import StitchSettings +from iohub import open_ome_zarr +from iohub.ngff import TransformationMeta +from torch.utils.tensorboard import SummaryWriter + +from waveorder import optics, util +from waveorder.models import isotropic_thin_3d + + +# === Core Functions === +def run_reconstruction( + zyx_tile: torch.Tensor, recon_args: dict +) -> torch.Tensor: + + # Prepare transfer function arguments + tf_args = recon_args.copy() + Z, _, _ = zyx_tile.shape + tf_args["z_position_list"] = ( + torch.arange(Z) - (Z // 2) + recon_args["z_offset"] + ) * recon_args["z_scale"] + tf_args.pop("z_offset") + tf_args.pop("z_scale") + + # Core reconstruction calls + tf_abs, tf_phase = isotropic_thin_3d.calculate_transfer_function(**tf_args) + system = isotropic_thin_3d.calculate_singular_system(tf_abs, tf_phase) + _, yx_phase_recon = isotropic_thin_3d.apply_inverse_transfer_function( + zyx_tile, system, regularization_strength=1e-2 + ) + return yx_phase_recon + + +def compute_midband_power( + yx_array: torch.Tensor, + NA_det: float, + lambda_ill: float, + pixel_size: float, + band: tuple[float, float] = (0.125, 0.25), +) -> torch.Tensor: + _, _, fxx, fyy = util.gen_coordinate(yx_array.shape, pixel_size) + frr = torch.tensor(np.sqrt(fxx**2 + fyy**2)) + xy_abs_fft = torch.abs(torch.fft.fftn(yx_array)) + cutoff = 2 * NA_det / lambda_ill + mask = torch.logical_and(frr > cutoff * band[0], frr < cutoff * band[1]) + return torch.sum(xy_abs_fft[mask]) + + +def extract_tiles( + zyx_data: np.ndarray, num_tiles: tuple[int, int], overlap_pct: float +) -> tuple[dict[str, np.ndarray], dict[str, tuple[int, int, int]]]: + Z, Y, X = zyx_data.shape + tile_height = int( + np.ceil(Y / (num_tiles[0] - (num_tiles[0] - 1) * overlap_pct)) + ) + tile_width = int( + np.ceil(X / (num_tiles[1] - (num_tiles[1] - 1) * overlap_pct)) + ) + stride_y = int(tile_height * (1 - overlap_pct)) + stride_x = int(tile_width * (1 - overlap_pct)) + + tiles = {} + translations = {} + for yi in range(num_tiles[0]): + for xi in range(num_tiles[1]): + y0, x0 = yi * stride_y, xi * stride_x + y1, x1 = min(y0 + tile_height, Y), min(x0 + tile_width, X) + tile_name = f"0/0/{yi:03d}{xi:03d}" + tiles[tile_name] = zyx_data[:, y0:y1, x0:x1] + translations[tile_name] = (0, y0, x0) + return tiles, translations + + +def log_optimization_progress( + step: int, + optimization_params: dict[str, torch.nn.Parameter], + loss: torch.Tensor, + tb_writer: SummaryWriter, + recon_args: dict, + yx_recon: torch.Tensor, +) -> None: + # Print progress + print(f"Step {step + 1}/{NUM_ITERATIONS}") + for name, param in optimization_params.items(): + print(f"\t{name} = {param.item():.4f}") + print(f"\tLoss: {loss.item():.2e}\n") + + # Log metrics and images + tb_writer.add_scalar("Loss", loss.item(), step) + for name, param in optimization_params.items(): + tb_writer.add_scalar(name, param.item(), step) + + yx_pixel_factor = 2 + fyy, fxx = util.generate_frequencies( + [yx_pixel_factor * x for x in recon_args["yx_shape"]], + recon_args["yx_pixel_size"] / yx_pixel_factor, + ) + pupil = optics.generate_tilted_pupil( + fxx, + fyy, + recon_args["numerical_aperture_illumination"], + recon_args["wavelength_illumination"], + recon_args["index_of_refraction_media"], + recon_args["tilt_angle_zenith"], + recon_args["tilt_angle_azimuth"], + ) + tb_writer.add_image( + "Illumination Pupil", + torch.fft.fftshift(pupil).detach().numpy()[None], + step, + ) + tb_writer.add_image( + "Reconstructed Phase", yx_recon.detach().numpy()[None], step + ) + + +def prepare_optimizer( + optimizable_params: dict[str, tuple[bool, float, float]], +) -> tuple[dict[str, torch.nn.Parameter], torch.optim.Optimizer]: + optimization_params: dict[str, torch.nn.Parameter] = {} + optimizer_config = [] + for name, (enabled, initial, lr) in optimizable_params.items(): + if enabled: + param = torch.nn.Parameter( + torch.tensor([initial], device="cpu"), requires_grad=True + ) + optimization_params[name] = param + optimizer_config.append({"params": [param], "lr": lr}) + + optimizer = torch.optim.Adam(optimizer_config) + return optimization_params, optimizer + + +def optimize_tile( + zyx_tile: torch.Tensor, + recon_args: dict, + optimizable_params: dict[str, tuple[bool, float, float]], + tb_writer: SummaryWriter, + num_iterations: int = 10, +) -> torch.Tensor: + optimization_params, optimizer = prepare_optimizer(optimizable_params) + + for step in range(num_iterations): + + # Update params + for name, param in optimization_params.items(): + recon_args[name] = param + + # Run reconstruction and compute loss + yx_recon = run_reconstruction(zyx_tile, recon_args) + loss = -compute_midband_power( + yx_recon, + NA_det=0.15, + lambda_ill=recon_args["wavelength_illumination"], + pixel_size=recon_args["yx_pixel_size"], + band=(0.1, 0.2), + ) + + # Update optimizer + loss.backward() + optimizer.step() + optimizer.zero_grad() + + log_optimization_progress( + step, optimization_params, loss, tb_writer, recon_args, yx_recon + ) + + return yx_recon.detach() + + +# === Configuration === +# INPUTS +INPUT_PATH = "/hpc/projects/intracellular_dashboard/ops/ops0031_20250424/0-convert/live_imaging/tracking_symlink.zarr" +INPUT_FOV = "A/1/001007" +SUBTILES = ["0/0/001001"] # or "all" + +# OUTPUTS +OUTPUT_PATH = "./optimized_recon.zarr" +OUTPUT_CHANNEL_NAME = "recon" + +# TILING +STITCH_CONFIG_PATH = "./stitch_config.yaml" +NUM_TILES = (6, 6) +OVERLAP_FRACTION = 0.2 + +# OPTIMIZATION +NUM_ITERATIONS = 10 +LOGS_DIR = "./runs" +FIXED_PARAMS = { + "wavelength_illumination": 0.450, + "index_of_refraction_media": 1.0, + "invert_phase_contrast": True, +} +OPTIMIZABLE_PARAMS = { # (optimize?, initial_value, learning_rate) + "z_offset": (True, 0.0, 0.01), + "numerical_aperture_detection": (True, 0.15, 0.001), + "numerical_aperture_illumination": (True, 0.1, 0.001), + "tilt_angle_zenith": (True, 0.1, 0.005), + "tilt_angle_azimuth": (True, 260 * np.pi / 180, 0.001), +} + +# === Main Execution === +input_store = open_ome_zarr(INPUT_PATH) +zyx_data = input_store[INPUT_FOV].data[0][0] +_, _, z_scale, y_scale, x_scale = input_store[INPUT_FOV].scale + +output_store = open_ome_zarr( + OUTPUT_PATH, layout="hcs", mode="w", channel_names=[OUTPUT_CHANNEL_NAME] +) +tiles, translations = extract_tiles(zyx_data, NUM_TILES, OVERLAP_FRACTION) +model_to_yaml( + StitchSettings(total_translation=translations), STITCH_CONFIG_PATH +) + +if SUBTILES == "all": + selected_keys = tiles.keys() +else: + selected_keys = SUBTILES + +for key in selected_keys: + zyx_tile = torch.tensor(tiles[key], dtype=torch.float32, device="cpu") + + print(f"Processing tile {key}") + timestamp = datetime.now().strftime("%d%H%M") + log_dir = f"{LOGS_DIR}/tile_{key.replace('/', '_')}_{timestamp}" + tb_writer = SummaryWriter(log_dir=log_dir) + + # Prepare reconstruction arguments + recon_args = FIXED_PARAMS + for name, value in OPTIMIZABLE_PARAMS.items(): + recon_args[name] = torch.tensor( + [value[1]], dtype=torch.float32, device="cpu" + ) + recon_args["yx_shape"] = zyx_tile.shape[1:] + recon_args["yx_pixel_size"] = y_scale + recon_args["z_scale"] = z_scale + + initial_recon = run_reconstruction(zyx_tile, recon_args) + optimized_recon = optimize_tile( + zyx_tile, + recon_args, + OPTIMIZABLE_PARAMS, + tb_writer, + num_iterations=NUM_ITERATIONS, + ) + tb_writer.close() + + # Write to napari viewer + scale = [z_scale, y_scale, x_scale] + viewer = napari.Viewer() + viewer.add_image( + initial_recon.numpy()[None], name=f"initial-{key}", scale=scale + ) + viewer.add_image( + optimized_recon.numpy()[None], name=f"optimized-{key}", scale=scale + ) + viewer.add_image(zyx_tile, name=f"tile-{key}", scale=scale) + + # Write to output store + pos = output_store.create_position(*key.split("/")) + pos.create_image( + "0", + optimized_recon[None, None, None].numpy(), + transform=[TransformationMeta(type="scale", scale=[1, 1] + scale)], + ) + input("Press Enter to continue...") From 481b4a7ff96fd3305295ce06bb864859587b3061 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Wed, 2 Jul 2025 15:40:07 -0700 Subject: [PATCH 18/25] propagate gradients through fluorescence recon --- .../models/isotropic_fluorescent_thick_3d.py | 24 +++++++++++-------- waveorder/sampling.py | 2 +- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/waveorder/models/isotropic_fluorescent_thick_3d.py b/waveorder/models/isotropic_fluorescent_thick_3d.py index cf67c862..6ea19a61 100644 --- a/waveorder/models/isotropic_fluorescent_thick_3d.py +++ b/waveorder/models/isotropic_fluorescent_thick_3d.py @@ -1,6 +1,5 @@ from typing import Literal -import numpy as np import torch from torch import Tensor @@ -31,6 +30,7 @@ def calculate_transfer_function( z_padding: int, index_of_refraction_media: float, numerical_aperture_detection: float, + detection_phase_zernike_vector: Tensor = torch.tensor([0.0]), ) -> Tensor: transverse_nyquist = sampling.transverse_nyquist( wavelength_emission, @@ -43,8 +43,8 @@ def calculate_transfer_function( index_of_refraction_media, ) - yx_factor = int(np.ceil(yx_pixel_size / transverse_nyquist)) - z_factor = int(np.ceil(z_pixel_size / axial_nyquist)) + yx_factor = int(torch.ceil(yx_pixel_size / transverse_nyquist)) + z_factor = int(torch.ceil(z_pixel_size / axial_nyquist)) optical_transfer_function = _calculate_wrap_unsafe_transfer_function( ( @@ -58,6 +58,7 @@ def calculate_transfer_function( z_padding, index_of_refraction_media, numerical_aperture_detection, + detection_phase_zernike_vector=detection_phase_zernike_vector, ) zyx_out_shape = (zyx_shape[0] + 2 * z_padding,) + zyx_shape[1:] return sampling.nd_fourier_central_cuboid( @@ -73,20 +74,23 @@ def _calculate_wrap_unsafe_transfer_function( z_padding: int, index_of_refraction_media: float, numerical_aperture_detection: float, + detection_phase_zernike_vector: Tensor = torch.tensor([0.0]), ) -> Tensor: - radial_frequencies = util.generate_radial_frequencies( - zyx_shape[1:], yx_pixel_size - ) + fyy, fxx = util.generate_frequencies(zyx_shape[1:], yx_pixel_size) + radial_frequencies = torch.sqrt(fyy**2 + fxx**2) z_total = zyx_shape[0] + 2 * z_padding z_position_list = torch.fft.ifftshift( (torch.arange(z_total) - z_total // 2) * z_pixel_size ) - det_pupil = optics.generate_pupil( - radial_frequencies, + det_pupil = optics.generate_tilted_pupil( + fxx, + fyy, numerical_aperture_detection, wavelength_emission, + index_of_refraction_media, + phase_zernike_vector=detection_phase_zernike_vector, ) propagation_kernel = optics.generate_propagation_kernel( @@ -95,16 +99,16 @@ def _calculate_wrap_unsafe_transfer_function( wavelength_emission / index_of_refraction_media, z_position_list, ) - point_spread_function = ( torch.abs(torch.fft.ifft2(propagation_kernel, dim=(1, 2))) ** 2 ) optical_transfer_function = torch.fft.fftn( point_spread_function, dim=(0, 1, 2) ) - optical_transfer_function /= torch.max( + optical_transfer_function = optical_transfer_function / torch.max( torch.abs(optical_transfer_function) ) # normalize + # NB: this is a /= operation, but in-place operations do not propagate gradients return optical_transfer_function diff --git a/waveorder/sampling.py b/waveorder/sampling.py index aab179d3..e6aa6a3d 100644 --- a/waveorder/sampling.py +++ b/waveorder/sampling.py @@ -65,7 +65,7 @@ def axial_nyquist( """ n_on_lambda = index_of_refraction_media / wavelength_emission - cutoff_frequency = n_on_lambda - np.sqrt( + cutoff_frequency = n_on_lambda - torch.sqrt( n_on_lambda**2 - (numerical_aperture_detection / wavelength_emission) ** 2 ) From d400ffec5820bf3b96ce787d9176f3fe6de43968 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Wed, 2 Jul 2025 15:40:31 -0700 Subject: [PATCH 19/25] zernike coefficient optimization --- waveorder/optics.py | 28 +++++++++++++++++++++++----- waveorder/sampling.py | 1 - waveorder/zernike.py | 40 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 63 insertions(+), 6 deletions(-) create mode 100644 waveorder/zernike.py diff --git a/waveorder/optics.py b/waveorder/optics.py index 845a1d8f..16bd35d0 100644 --- a/waveorder/optics.py +++ b/waveorder/optics.py @@ -4,6 +4,8 @@ import torch from numpy.fft import fft2, fftn, fftshift, ifft2, ifftn, ifftshift +from waveorder import zernike + def Jones_sample(Ein, t, sa): """ @@ -154,12 +156,13 @@ def generate_pupil( def generate_tilted_pupil( fxx: torch.Tensor, fyy: torch.Tensor, - NA: float, + NA: torch.Tensor, lamb_in: float, n: float = 1.0, - tilt_angle_zenith: float = 0.0, - tilt_angle_azimuth: float = 0.0, + tilt_angle_zenith: torch.Tensor = torch.Tensor([0.0]), + tilt_angle_azimuth: torch.Tensor = torch.Tensor([0.0]), slope: float = 4.0, + phase_zernike_vector: torch.Tensor = None, ) -> torch.Tensor: """ Generate a soft-edged 2-D pupil that may be tilted on the Ewald sphere. @@ -168,7 +171,7 @@ def generate_tilted_pupil( ---------- fxx, fyy : torch.Tensor Cartesian frequency grids (units: 1/length) with identical shape. - NA : float + NA : torch.Tensor Numerical aperture of the objective (must satisfy NA ≤ n). lamb_in : float Illumination wavelength (units: length) @@ -180,6 +183,8 @@ def generate_tilted_pupil( Azimuth φ (radians) of the tilt in the xy-plane (0 = +x). slope : float, optional Controls sigmoid roll-off (≈ 90 % change in one pixel when slope=4). + phase_zernike_vector : torch.Tensor, optional + Zernike phase coefficients (radians), shape [N]. Returns ------- @@ -217,7 +222,20 @@ def generate_tilted_pupil( # mask outside sphere pupil = pupil_soft * inside_sphere.to(fxx.dtype) - return pupil + if phase_zernike_vector is None or phase_zernike_vector.numel() == 0: + return pupil + + norm = NA / lamb_in + rho_sq = ((fxx / norm) ** 2 + (fyy / norm) ** 2).clamp(min=1e-6) + rho = torch.sqrt(rho_sq).clamp(max=1.0) + theta = torch.atan2(fyy, fxx) + phase = torch.zeros_like(rho) + + for j, coeff in enumerate(phase_zernike_vector): + m, n = zernike.noll_to_zern(j + 1) + phase = phase + (coeff * zernike.zernike(m, n, rho, theta)) + + return pupil * torch.exp(1j * phase) def gen_sector_Pupil(fxx, fyy, NA, lamb_in, sector_angle, rotation_angle): diff --git a/waveorder/sampling.py b/waveorder/sampling.py index e6aa6a3d..950b5d0f 100644 --- a/waveorder/sampling.py +++ b/waveorder/sampling.py @@ -1,4 +1,3 @@ -import numpy as np import torch diff --git a/waveorder/zernike.py b/waveorder/zernike.py new file mode 100644 index 00000000..bed82663 --- /dev/null +++ b/waveorder/zernike.py @@ -0,0 +1,40 @@ +import torch + + +def noll_to_zern(j: int) -> tuple[int, int]: + n = 0 + j1 = j - 1 + while j1 > n: + n += 1 + j1 -= n + m = -n + 2 * j1 + return m, n + + +def factorial(n: int) -> int: + return 1 if n < 2 else n * factorial(n - 1) + + +def zernike_radial(m: int, n: int, rho: torch.Tensor) -> torch.Tensor: + R = torch.zeros_like(rho) + for k in range((n - abs(m)) // 2 + 1): + num = (-1.0) ** k * factorial(n - k) + denom = ( + factorial(k) + * factorial((n + abs(m)) // 2 - k) + * factorial((n - abs(m)) // 2 - k) + ) + R += num / denom * rho ** (n - 2 * k) + return R + + +def zernike( + m: int, n: int, rho: torch.Tensor, theta: torch.Tensor +) -> torch.Tensor: + R = zernike_radial(m, n, rho) + if m > 0: + return R * torch.cos(m * theta) + elif m < 0: + return R * torch.sin(-m * theta) + else: + return R From 83bb026173942fd9f0def0515b5aedd5e8992432 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Wed, 2 Jul 2025 16:54:00 -0700 Subject: [PATCH 20/25] comment out biahub dependency for now --- docs/examples/visuals/optimize_phase_recon.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/docs/examples/visuals/optimize_phase_recon.py b/docs/examples/visuals/optimize_phase_recon.py index a649d064..acbc4854 100644 --- a/docs/examples/visuals/optimize_phase_recon.py +++ b/docs/examples/visuals/optimize_phase_recon.py @@ -3,8 +3,10 @@ import napari import numpy as np import torch -from biahub.cli.utils import model_to_yaml -from biahub.settings import StitchSettings + +# Commenting biahub dependency for now +# from biahub.cli.utils import model_to_yaml +# from biahub.settings import StitchSettings from iohub import open_ome_zarr from iohub.ngff import TransformationMeta from torch.utils.tensorboard import SummaryWriter @@ -213,9 +215,10 @@ def optimize_tile( OUTPUT_PATH, layout="hcs", mode="w", channel_names=[OUTPUT_CHANNEL_NAME] ) tiles, translations = extract_tiles(zyx_data, NUM_TILES, OVERLAP_FRACTION) -model_to_yaml( - StitchSettings(total_translation=translations), STITCH_CONFIG_PATH -) +# Commenting biahub dependency for now +# model_to_yaml( +# StitchSettings(total_translation=translations), STITCH_CONFIG_PATH +# ) if SUBTILES == "all": selected_keys = tiles.keys() From 6cfb49d46f8ad79de2c8da166cfa1fdd246b4141 Mon Sep 17 00:00:00 2001 From: Shalin Mehta Date: Wed, 2 Jul 2025 22:16:42 -0700 Subject: [PATCH 21/25] add tensborboard to dependencies --- pyproject.toml | 1 + ...ts.out.tfevents.1751519743.login-02.125407.0 | Bin 0 -> 59371 bytes 2 files changed, 1 insertion(+) create mode 100644 runs/tile_0_0_001001_022215/events.out.tfevents.1751519743.login-02.125407.0 diff --git a/pyproject.toml b/pyproject.toml index a01928e9..060b1b32 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -68,6 +68,7 @@ all = [ "napari-ome-zarr>=0.3.2", # drag and drop convenience "pycromanager==0.27.2", "jupyter", + "tensorboard", ] dev = [ "black==25.1.0", diff --git a/runs/tile_0_0_001001_022215/events.out.tfevents.1751519743.login-02.125407.0 b/runs/tile_0_0_001001_022215/events.out.tfevents.1751519743.login-02.125407.0 new file mode 100644 index 0000000000000000000000000000000000000000..2a3b423b0172dca948534ed47a3bf5bf3629e33a GIT binary patch literal 59371 zcmeFZcQl;c_cuHu(NhqFD2asVL>WYlkRnzsY|{_M{#$9>3uuA`Lux%BIQ zzC=WGl5caJwez-hc9+odcXao?C2c6dAm;1n?&INY=V9yZAmQVOLfLu;Nd$O1`#O3{ zIG)^PgK&0pwE52`E~5QHm;O6$al}J@`gsZH^OJE!&N47EoYwL1@nH+5{xld9%9h0a z-~0Nz(s}5ANd5OdW`^?&bipDHRW*65K*cB-kZKNY(uYk^S}2BFkEBc zcK1U$dOO?Oy4l!zI(qy1c{|!TIQlx;`#O8LvvuZksd{hJM9KX3d8V=}CiDgW`*|S- zK?c76m*<_`-270^?zaC95%qph#dqCPMc}_fthPb7>1Y1;5L^tG7?^#X-F$6q-H|{H z!H(|EzD{ha3==8_-9JMF|ND@S{#OV+i}ZizeVKuif#v^r$Trv+3-=+7kU{6~U?UJQ9&3@i`+kI7vz@bh$bVDGHbVDGLdn zJ?1=+YV4wk+SNK&8-=e*`G>0Hqd7-tH=*#li(rJ-qHXqnCiXwS|1U`3T(8R<2c%7@ zatwQ%kHnJNdeb3H_Z6tEXZ1hU*j#O)>eE!94j&sA+FgL7T)9{$B0mV^Pz^yj8e5J# zlJJes8e?9K9hEgpvVV^Jz@I~X`(rnhw!fM` zA6{kSDCw_7OYNIKnS5?xIh~H?I&SM_kdUk3jea6i{X3$ig84mB&HN+@US8g9mtZZMD!_eIEYHIrF#)CdZBgyFW2IT z?>p#_KP1G`23xf1nq}Gb^!29__eZs1dbm#aRT#0t6Uu=L|NA?c&zSFPr^&*y@0trt z)%9`#Ki7|ndb#<2@xHnk26*~g{*YX7sfr|ScRDDR`-SAT8WX)qRnYz&(iOms0aw|sxSh^j{4f;C}q+zdUH0h^HK zc+Vda_WMU+=quWvH=ko0-JqsYYj0^XgVV9j9pU3BcN1oDmdU-%x$VBc&U$W-LU3V2 zzx^u5DyB9*e>-iRGswyPa$~wq1l_Pbz@JL#o_xNzu)wbho(j%ZL~C`qxu?LDGYv5f z4bJ0k?zvCV+wlF&uM>NVQ%xTGUg=}xG_)!(!ePJC5=-5hn^|raSCbCy!cegfwn{g0`=tZZMw> zyxPsl-8f?pmq>5Ws*8(@!{KnN_H4uoRYX%1JbiG<4jw{*%NG&kZprIsF#0k&H^;J{ z!Fu}mq^d1_GWS((jj4Z^@RnP*tRsw8MsnMdfd92g zQ%5mFqvD;6WQ0*`9hZHWf#$8Pt+9L_>=P-WykKu9rVK4qEErB%i@!G^XL~+o|C5Po z5E~MSn{}kOt*&^C#$?Xu;>Pux9!!1`2gRM^cfwgC>YmuhSgw23bj-Z{_L#n z6#w2cQ@^lA(-Keg=g$tM`L0J{I`=?NaYhxPR4dk3*KfUL|h$)^m^Wqr(pWG`z0gx_K*4;1gxVsU(nnf#IQ#5&hZ8 zTu++kw#|v0JiSmfu-Lw`qR0M0L7bx075hgsV;PlYeRdmWH>L#Q2`X#+hI*l7RkQof z>#IbHhW|b4oC6;Q9{U1;Yw?2L*j7*M=-&IUn@p> zh(}Vc_N0W8x)rT82c?d4$sP+xVuMmy&jaVNPef?Z%=C2V-`_Xy)l1rRSXxk?_WP=sI61YE zjqe-DCR+7vngw8n*UYmL*yD+-kx%^j$=$3}28vN$eim8HY3cSPV6iZ9QRnYA$A3&e z)>%r$^<<81Tw=u7>F63(y~Va~O@!#yRm@~{dP@{lPppZQq$~LO`=4mf*247cZ0}l- zkB^b=#K9H?N1-@Aa9&2G-K@YWd|Qa!_yCaS_bs~N%GQ?1-C zwuIJ7uKHNFRjV00-S@q>Zj!@Z&WHG^27e=@eULPs&n@25oVNSpho`sq7Wf?;Y_9F} zvm5A_b!nd>y1 ztvX^jj&aou>*KX{>2#4f-%i10EVKlyg(Xi^GTaa{U&*62O)sav&fKYzY6`k|&Gjs5w@A=In| zn_zJdJO$nwYotoeFwb;%L>JMV&MAZU=994_ZFU~>I$}Q!;#rR?dopplc;}nwvzUPv zZ$Z~%0W|{m=Rr{$t@_3!ld(?=Uhdo(9%2IZ{m9Bb5<-z<5UGAOr|SaQp)q{g3MzdY zCWRZTn~OMBaY5YN>lDdsxFJbk! zaBu)Eq3A9RMhWHhtmRjx*J?iw{_e5#&`={PKu;k`r7{Kn+0~~8PYo5+ojTg#1JeWInPleBZsG51NBQ_<%@rnEbTJpfT(Bi(*8r=H& z<<4#vy~ixX$2R=DuxnvVKO9dP;SUp^k~)tMH(7|`K$oQVtEvX^y7e&^uAs#+aZZXi zC_y1VzbwAF?_8$u=bANEBQZ5Ui2UJBkNk($XF$wD>T)XHyQ@kYF?JnGc?xGcMed7^ z5i#HJR#jFu-od&^m|M$BaS#{e(N->-o0}7feQf>6AFA}oALA_pVt0QTLgF4E9Lrv% z3>^LbnwWwgM04I8a9f)`eq$=difeSaX1-%{Ex$ogY2vnzChl^!Z}Il^{H|7DTNrXS zUYj6h;&`_~ECv21t0$au#ss0X8t1519qRd@Yk`SIBwv3(_DC*3c~)+u{TDi>-!n;v|;T-H@V`oF3iOf1P!v- zGWH$KK2t&-c+|Gf!NdjOxN*!E@%Lx`@!zuQEJ`lEB^57t20Kx2pn>_MK@`a$Xw}CI z!kKBk^b-7@2Ha7IGBrI@>e3jfhj6q-=pn$P_pW1I5=qSu5UEc|qHE~k1usR4b1=En z{TkV_s|xdeYo96>>G>;KCQzwzX&4`h4!Vv#{=2M(8>kB1oR+G2ybG;hF*hp{IQ28# zf^dkgL-hf%gz)Edt(ljYAdC&R^stX)M7Qxr+vQ@PfYVk659ZWM@7cbMhRRhclv$hL zj^#w>F7*ZXcH)pf^Z-MQtpKJO+^_!`Kd20B{8b&?n#v#-_c^!rI(q$~L}7!`@K%Pd zE>~P%@&~9s#U@QHedx;zpaJG&2llfBZ+Eb?S|}ES_a^8_Up8p!fggU|vZ=t}$q~N4 z1>0!5Ioo1m+!0wW>;qz(NRyVV+EKSAOmlF{argDC_(B)2-p-76sX;!82(P1h1_#H@ zS6j8!Md>p14$=yC;Ndo?426;N0Q)d-(o z;PZ+3e)SE%d8jUZXzG$6x9~iZ^`ivG{m4?lkIJ*O6e{(QWRFPb zqrarCu~>vUHtc7G(%ieN5XCW$azejBB;$(;NrG4&R|Es%+LTJZBXj2xzxmigeIRQ> zx+jg~n^w9jXk@Q3c2eP5zLRxJ))-l5;Jf-d0168keYw_6=WE&OG=17&KtKmoi-k+G z!~Xmd&OiO*V6!toADMzbH)*ndf(3js}=cHSV2o)!Lz)iF}^By?M*3l{l|8FKM zp}(N{C@2PS2=||2cJO3Y3#DTn=_q@kCvh)BO?XhV8!qF5SVI}n`KpkdW(WUa6X&7H zzh%mbiZ5$kI-q?fXNvDTqatKBLj{Yvm@g2&oH0I2NEzA+$L>Q2EX0az()06t#GxS{ zdF+gWJy_1%>dF!F{Q2{bG?bK-a6MQ{C>T2Q``h!w)6pldt;u$ELi<4|5JDoK+EN`! zL)!I33<-Avj^H|~q5}1efhN*$sQhu2G8+Gd4lM3e#WO@meUKv4m%V#n^Ti2vC`2?C ziq=gbJ#xY%CNIgi-#O*2(AwTE&Q}R{s-?ZR|0o;dc{*C=x4_IZgx@7~cPEnpG07pz z8_D>sdkN}9eIs*sxVBS}m-=O5tWaUYjtoIp;aoT8-!{rKM8ml&9}(CRjAD~l@AWL6 zHzN6LV+#aJMbZ85yKl_+huy;kQ@aXTfh07gMd!tCU@P3viB$ap zP_uq?bQE~j^L(m{a+s#3ra3+>)wn-!>t;paysoR@=^G$up`>OPApXy*=B1_CJl-FSU>QNMJ6N~;M{=^)a?%}6?tg=*EL&LK+@s&BAcvwc zUtvC22BcO)IF`KaRwoTPm+DEQR&0K+c{P{Tpm7LjmHBA{7sO27lJEC{@9?3GZ(p_Y zV$e78)0dVG4-b2Kv|!2yyXzL8hrE5TVt?zw2-%)AL*rJo7aXyZs@H&j#4lWrY@;i% z5Lc>m*p1UKJKQaPaCxSOshsbdv=RSa+bBZ<89bj%}V(X8)FOMw`ddUmK!BegTP^{b zPf5pkxq^+c1pzsyH_0u#bn;ZO`D+rIpnXS%55;?y@Hh4Rh?8k}K}!;X#>xzPr2o+nST>;niv#Ra_kPV0^Gy4Yd@QgAch5j)#gqq|#P!1!f$AG&Uk^ zCC+?J4MzWz-1)-bj3WSRS~&Y2gq;sb-~aq+QSK)Ywsrva7XY3EYjQJ5A+4-9XQ9v; zrvQ7eRZ$19;u^8lkw?aI5Bjmwa*UK(`(UH18`UrOU|B}&ywu%7fP!9SA$}pDk(jXH zG#p}N*mBZ)_X_v;*t0t?oo1F0L?Gz&fKviF0E|}FCTd^`aSQN{zSD9F7uUjUAl%QI z>)<)B;XAU$ZsWY8pW_|BKF^E7C09qiJF0kfYeGFy8Wp5wFiZb2_DeIh#gixyu48*1 zXUXTqJFrJNatcq45~lnhQHi;IvMqe$J_7?{vrP4v#xx6e_g3_3SPzv-SYXr1`rZa= zN#j%qW@Vg4gasx-TP^(P(LuQ06wgO3)rl|qE*T5(tRVCpi*I6qfpXNH>p5Nim58VQ zW21PLQYR^jjBI7z()4*gLjS~W+&TUl<4xKY0->_9vOM2U2{UH#G)(QrW2Jc}|GMA` zD+YtJNsMP+MU^|eiinZ4n_Ql>HamD*fO03ciUlrD?%dm!?<7KTFi9uS@Li#8d);)?+b&Nu57sah zgYk6lxF#hy&xH03$;LTi{Vl5)5Gd`-A0sn<_1+CRvlNS@*jtNDvk$cx@s`!v3mMQM z|IvcgnhBW$1PW726z}cqZo0X&Hf%!Nf`+VzUBaDW7eU&q8lBL@+ZNEcLHMf|C1vp5 zp0;z3-LFO6y!Y{Djf>5*V7Zna>@=xri_~4Mu*E~?(9P*KfhKGTaycr9-L6&|LApGf zETri#lt5gd?hu3_mAxZcr3=ys?4`Lcoum}NjE!adT3CtQ8r61_%wF#w7JU*tNSYHd zmUpNkxs~M-4jmp}C_CtYT>7k`FEi9F@xk_@2sxK@kC(ZxIfbg)q>_NZZdR7VJlGbV zRhiJZg*6$3lOy#nWEvA6SLWcCTJG=TbK0VwWs`1SUKXgx{yNqg^X`T{ufy9M$_&x1 zoQn3gLKNIbg;1fLe6Q;g-IQ75s;=r(#%&4vQ0fkY{X9z*E6yawmLM`*IPa`=-p7u-$)>=67grOT%q$<=#ayve%C8eC1eMa+=qGKS$yio_e1PFN1PtEYAHaS0w^Cy+!8v>D8FK#!zscQl@64+o zuOcPpnOsm-g{_BfocpoOtfj5RE4=ZZXWP;$maW&fS=jwYFN{3a>E^_{+{l5YV{`C! zjii6lUyi(8{=&3#3XHfaTPBhtp8r+$vk`T)%(rrFMi`u5e+{08Lg1C8Lyn{TU2De< zp1Z5Mo)0Nb?Ie)?Its4g5{VD`ods1Y^NBy;N$I7LU zv_L@pTb*A8)tmyZsXB%p8N|06YHifyjrgHta_10Uc_1z?|J53W*<7#37Uo@S^Enx7 zUyS4WR_27e2XwD^x2NEMSFA1hMI!3u!~yLkFwSueS9t|hjd^v}L*^E(O~RFk1s)8B z5oc0ro5y$oEW2N~sL`9f*orCeLvP*4+vgm$<~jc3;=k^uHyzyV4m))R#K|UdFK;CJ zkIfLx3z>qYu$8X^JAHOU>(1>|$hc&}qE4zUswU*Ltz{dd_W|rs}MF@3@JA*SJ5F zZ%oAuWc;{-$4z;Zn3G z4Yr~Y!Mn_E;>5K)SzM83LBaanqV=CMWxxZK$qs($5){1#eZ~~HIqDR#x7vZ)oQ0WJ zE0TMtPjrF1afy7#C%(`+GLb5aCu;OTA{UA1Hj4aRglSH@qiTO|FQfK8m-k z3f^jEq*yD^b#scr`;ooSXAL`HT1^{_Od1MN3C-j?N@F^tS0js1O842eP26b&2yY#B z6#Rmle_XI>t67V-8`HUlYjWMIjyYhzdo8EyxbjwdCE|!QJ&B06mz`Pl=rEnBl@gn0 zIx4fEseC5?=+)cyx$}%cxDa{;g>FvY`T_ip-a<2#>fd^V-05(~y{RW*K#SsH%%G1o zFRuqr#_dYvw#Ls$NnM#|YOEBwFq%~RD!^$?_nSn_iTg^aR=G0S+N|W$Tt?6IP$8-# zCK*rHH_tNiIJBGdy;2>?Z8wbm#J%9NFZVKeFga#i_S<%vUwd*JK+VKC`)9Iuu{(MA zEGC7E0I4$`o*RoAc@+d}t6e>smfeaV_nkTQ5IEK%o=H39%g~~qd9XcXxY9^?%Mjxk z1Q=rE++#-Ov%BFkVmg8$8`Z6Q%clUvRvXK#pJtPey!kG)2E_aCp3g^uPmjr4cX+2j zg=2|E@eWT^|DI?1P{rjy^QtDQ5U`Xge>YbU)$bg%7}LR|L++ECSl!JpB>`iY>7JaX zLk8N*?w3}B1Tb~q#OB`9CprpI0%Py}BO*G7d7mD-v}U2CaRHl)<1tt1k;6@Xfuh%@ zWoJL;dcBQqGM;(9UVYg^d?_VMiIxUrjaXqSVqA9|AGb!-E3doNzw52!T0cM2QS$=R zNI<4evGra9+8+>Yk=JUUWWQOl$L2M6S9NWd@}!UsDWw~C@80rTgSR{8Gh26>PT3~eb9=23VWAwiFdGh+o_Nye z*YArZ`lzM&=S*#;Z|A`%QyZAriff%gE#|D3o&)pomJSIF{C3XUCzn9zFr~lbHTT|0 z(tEnrYwjq}#hruV+_rpbGs}~v!JWCtLJ_y`7HsY}wI=iBcf%$X)cI)y3WOOMXa$-^m6nc<`a-GV)28=!7 zQ>5gkAdgeB2id#T|AzLt-(rw?ro2g6##7!r3lhS`Rsl3aS4N3}H(KIzlKEBY%R_S* zajTu&Zbci)pbjF?V+wPnNFJqmY%G$ax`-P`akG#u){S#)uq#ROr1^dTdY|eG| z{evWPhTO~;s&FuwDPMQd4G?<;(g7XrO+rIoh4-s+2S?+bLimn%!>Rgg>odGuw+|2Y zHo-`LqYqWj6bJxSMU|74XS9oXwp;eyw5gZ#@)+)wTda|PjXS$nZ?-7fngJLwrZx+h zbaiEk%3qTGBEW+#w0^&PWcg03mGM+Xeq+szQ}riCyvhkXTMQ;;iu}qnMdDuMI`{b> zNvmnMkDZVoq!HJMY2}!8K0{wY%bN-kBIvePvrK2$ z!T4;x9p@{33X@CFuuEU3VV2M}8_vIjzb8~^!xD$=h?%Bx??$ZaoAVG7YrTif|8PN0 zJZxC=p9W`kx4R@a+rO8wSjDqGoX_p>8I^;`pj^Z>y1C9m0K3Z5f0Rc%Tv$#2u;S`J z*`ukV$VbW;%xE!opr*f2MAhUWmG_-)wB00Lu3`XiUs(T9=XQr zt%KRuHq9}Dk*)$p6hoQvmdJ|V=6(WV5%w{4vveW5i})U_%ZI1VpRao2E4e;qN)nge zJTH6hEHRxLob9O&94=ZnGeZ~`d4|{`Gt`nxu)X)uR?^L;H}rM|dr-rkwJ>Jp`npTxmQ`@Jlxsg@U3BYDIfp@ z*RAnBzYJ7Ge!RydhwbP$D5-`I=4VMu{T_s;z|&Y8UL1;*7FtpnVn4@p$@3NWkCfR@ zk@}s0*=TDh9_OfpY*Gw?m{{!7#J3BqU$ zD6YfgZ7sU089Yvc8kUprLkYvU{@NwoascUbZ6VWsG1GLXqBnF9mF}9FaJv6F?4we2 zjZJ5-?xgsfe}sY*SD$>B|IHLisEHr6GyIV2&eyRqm~dWh-pTcLVr{JWch|jwW_7+p7v7Bnnf#WJZe+&?LF;VaO^L;Mpqnh2&9E_%W z*47LLhtVtQc}i_az~b)6?UYTCazyjeZKXF)1z)&EDO~KIi9dct-dZhfGE^8ucun*K z&u<_?4uWsJ_W_*2DzCY~V5^ zy(S2`Wt?jM5kvHoFV9f@PzIIr52Or98`bc?fqs5|K|%ZTAD)9W2NsJ3IU!JtFys6u zD95W)G3U*mqk&plC5lH|=*{YzN2qmn*RN>TwFdv;leUL0J?+>EZ6_a@p^ngf!){Vf zPY*~ZX)7oxD#B|fKxaW>TC(3cff61&O{D?@!}NDsx$N>5*0dJ}s-m&Lj`xr-=ij@D0*ZM-;SH-w#@R!&=hUs5dMv12F6GHQa`$DFQZX_4IAL$H6hpZo+F z4`t@pVE6&#O_#eAFcGkDXihr&vTr1T;8qXb#QYlNw8Yji_Pw_}E9~g6BKnl^japWXGT{r&xk8XHLoi6U?I9w2zENDj|!AM94I@VbKZ zy{5wEEH$C3&vEWv2>|JP7F(4&FqF=wcG^`RC{`{E#&{J%N<-|cvA52XvCxoIK_8+gY zhMrI6&-jVQ$!K#inBaT~>@M(lE!$4qe&4D1t+6GD zl8-{%*&%PFjO3J9VwwFbN1w%K23qwSwh&MP1%*x(f`2lC4cA8W1di^v0H$x$Ph0hb z2JC}W*wWv>fA6r25{sI%){r-|d+V`#sh?)6mB+;NZC~Xfatb-m5{^{|m*r+=W_tJu zy0Hl_*%VJu z2B}qhD0WI{KrDk`ckrsQOmP6|nF(g|KvA(ll)5A0INaagKP=2D2HSkr+2?0fpobMP zphRg*45T9Mvq7;zYRaus|LkwSlU(73TOpYKPVnLSR1>fhOErrPH*~&HEHWn4KW&Mf z3M9^2Kw%C;AV7_vQBi>A3{5G|m|R2TrJt1^9E|+br>t*BKv*c_)|g5!;uf*RBgD}m zN9=B9X%7Bqe>p$K!DPMQYF5`gFbYRwQEEGVGDwSy6j3~(^4nGB+)@L`dxhG1usbpd z+5HCuw8QR@+VS7*A+<-)I6ugYgCP3ftZ#;|Db2@00w}jOj`bfBVm!yKxUZmjvU2cS zD_JuRyZhVA5xfm+5LwxIA+WoZjZk@?pJHqNP}~sy$cYPZhdC#@K*KhlF8im}u7BcHzoj1!F*bC)N^K`*B6(>kjo(~8UsoxIhsY~lL$KUH937T}w0suD50Yv|)29VW z9!7HV+j#7lOZts`Sc~ha*wNAB!^C<5^5^Bb4bD1j^?+CcRXG~zpGw;5MQm|$bAn2Z zq0+t>zxfLmbBBiMy;s@r*An~RzkJWbH;&^6AL#nTL9WXvXoF%9`aWXx>mt_W2(DT0 zLd^O@x_$5xz!CoF&sq&2?z5aeXLbLrjdWKi9$?Cp@_|&M%mW$_NkiwlB?YmaM}I1{ z5w!t#Tj$)4#6v;{-^YQhJOxG|2v^VL7{*~Oq%8W8H*VaZJPC!CQXIxNvL7vFrH^Lu zXC3ZL*sRZTHX{%foZKp4&OOr$k#&9=m|phr(>m=`??N9Q$#&>yhJ@^&PqR1jyJch{ zBRcTC2D}oe(}_;`GpqaQq&m=Sg7BcG#ERQ9BEp7u8GdZRhzPo#KL?P1g-JnpOo9G;Sh~O!p-vk*4kHS}@zzJwO8A%qw<_zPzeip2I=$Yq))X8OO z5kNZ>^afUXsLG)Zybz>}E1|>}Y}o?&>pLr_JV3(CG>wDSEl*-L$50G&Fgu>6VUc`* zfvIBZ6>N;g26OWCq^Dkszu+)|q?ws~?B#UqBzb$Qhpi8eI8XbptPQz zt5AwyAiwU}A0$5|{qjvAO|R&bVD{&ZkD68yM_|ddC8Y7xlZ|r+8`nbkyFhUFULD8V z^(UV8^!GQ}BFXBO*mC(P2-Htt4W#1KMh}an2fz|`0b;hg-pxY-6v>q*I_(jY@LzI4 z2p}%Dh~3~=mSuO${Oi2`6yU|y6(|-Jk9awI{pmn;hTqQCvq1{K!{Z1 z+eY0mcwRSWnoUqE^r#>pK|db`@pf;z+(xAzUR|&7P!%IusLN0GsRs`WXR5r8$*Og{ zs|YEm&PsnuS{};xkzU)#2B38S;%y6y4G)ufw{732LDwU>M+P}G<7Aex{`AX&x~i-r zOUAZN&5uD>i@_25$RK{8>bvaa=h^o}W%8+ptuVq_U9AA$c|l!i#W_fDPd4J$^xT|{ zrR5?E8F(E_dhr)dR;?fFnoSVGnLTJun6&75?Jlu=C*Crd2Z)otjc_~#8}k-EJ0iRI z!X)RU4=`6nSsB0q(PSv&YEcZ2J!xrW`nW9>4n%K(0JSmDQBit4lzZ~-5CGS+->uau z5ppwf{6$ulmZ6{#@AJ^0*h$|Cl$T~Mq1{-LkEPxu7OLLogr->*&I7=t;04_YOn8%`a|Nvf!-(`={OV6 z8o*4fcnZ&R7|%VD%%TD#z==r8&)|0ALxm}`3(V`gZ@9@4fF<;BVu!UeuzzD?>{Svw zlE%f;b!DpYbo%={i5#7v^6J1oRyTCwx&Jm_#FI@BI?}k8ZSBfm+5uBKh7$I_De{6g z6>UnXP@;#2N6p)}gAs>6=I;m>rb{5RUYKY*DHfqVu@GHG<5ylVIQxJWE$1s3#Gg!N z92u5XlFliI{bs+^eC&k|a(7f!kY{xQS&L|qeBH!D0yL*V^faJA+H^u-zPr0yTU!eY zMZSuaETA_mGp;d9abcyv$XtHSILl4@YM^#aXj3AIcvzDH58{@=gFW1I7}t}=T$i~5 z=6VqWQeA0iy2Tgh5Ae4Oi;19<>ZEtfrqjqobddO!KM`nOdSUDR7}<>kh%y48}@n z0hyC1+6i0lD{HZWO#KR=jx4j9>y-8^j=bh1{iLn}nc8rqm_ z*|h170=hSd$-lVvC0v&uoz4L*-J^s>+MmZZKP6%F1b@r|Evwhfeu~EG1*EGt zXkPi_C2-QpxHDR&4_rZ<)0ia?()!KZAV1{ahn9aeC~g_*jY*Hmx8QsB@~JAaN`IF8 z*yebx)c*YWtRaS3vu~WBk-Dc4MQtE}ZHyi_8L5ajEWWw-^HIHmeWALJrYL4SFE=;+ z#kk%*S6Amz9ql;R@Z^_{qmmAFx~q>>mGfS~#e4CRrl6Neot+0i{V^xbkpb~u|FTyU z|6@G~UD!{dEu#m<$VBbMmUr*og{Gy#VHy_hzDsG1QrD_0@*B>#t)j1PW2ExvE^Z$m71-|NRTe?ec{IUS4K!TlHZ7n-8T=G$d^pooL-(>*mVm zlY)<229D_pr66_NVDWX^ue#TjK*6+e@M;SU+(8i~&4h5j68gnyDi>ygX)LLy+#;o= zz?+QG0b)2Fyt^1~Dq(Tfah9hzd&`F1-sE|Psj0$C$M@Uxpe>hO+X=-c&TkI8e*jSB zfju+MTXjmHsNtbflLXIXZ8UVE^;I8k9j8oE>#I|^p?U+GGqf;QX$ZG|sc0&0Bt;#5 z(uTU;%v`v*INbby8g_bE#opurP%Itq&&Zagd*gOaZVQYqz|L0 zL!UQd(3E%`6KHs)3&7un(dkw^c98ti5eNDs@R2qk>7MNL&}kI!{iHF8{dToT8H_)t zC*d;k!R3)g_D^J7o~vVHN!>@!Fhwe0Cmv74ffK<=!gv;$pJ9u5lNs(%K&HL zVcbv;t%RNTJ`{V7ZBrsA!0O+9srdc58m(nsaVrxh7I#8pb3T#PO*-U}broU4Xetko z&bwAGSJWUZ9+#Oe>cYqe_&J}Hl&OYAdA_bL;XC{cwgRKLchj`L$P11#9Va;0asHr* zmQgTgWUf{&_ZP-B)K{#KdK*NNC=Tq01b>^89^KzotODCVcGHIZYq zuKjAs&03~F6~H(cD`#Vg&WQ68^A#1jT%64%|ZNP^~Wq|73ko5JUucvt(y$kn1P++%6ZOWa{JXV1lCOwUQO~A$(rivxUadZSS8x_J6S0|{{rVF54;=G}TC&817w4HW#|Zi; zt?J;moDb|NHfBV9M%{tt3*+%N;A@cTpBxx3*qzp2J{r|Pm~Uuod@8=!9|@8c@3HHb zu>Sjh@_$T!7D%FcuJ!LBlYf}?!J!AyFYVqZGZQ0KCRN%8tp6_@ zp!uIy{*MIyM*{y}NMLjFnCe-`&xIs()7Hs&547|&-)cO5%5%70{rT2^-cV_w!cPCz zTK7MP7yJDfdi)qz4IS-0+BS&*U1vVpA7ylG&r-@Y1&#T&?r1hXLL2Xf5v&KsFX1LL zGP@>A2B2IloUAPN=L@@%av$~bY788;MR0HCR<`$oT4y6m`)FJLN3vz$v~to2q7Vl zul)ahNqymiY2jO{#H+SS$7mA`)d0^gyg8VMwl;*A@~C4AQ2#tgQpX`+-(t0iDFFEjJ^%4iwLq=$rrv^F=> zkHp^AP^6~HRrC=Q6zuKoMfAk7fzvg{YwhCZ^-?l2t}`t=`4LBo2z=(pj~^RRQ@l$4 zRmh<{*YPUnnHE1VJqeo*W#>af2JF&L-l}&(rfqBo#F`509Ly1~V<~^N_ub2S=P-TrK|V z=HVDA2V=-GF}FBPHSkBl@TSLHhq$IsLDXBTR{zC*W^0p?6oH51Rd5+8segtyb8~Vo zwwHjfz-GiYjmg@t1dM!@{c0+x zx7L_9xGa2p>0%;&1P;IK9~|&0c&b0>U1{z?Sy?@PUlH{6xpxTkH0U*Dqy|#B!xB(;90pwZ)|+L+v>*@cai_ z6iaVk-{x%lRB!KH_ikaMQg9$tH(S;jtcYA+B#=HxDtVZ(I#DzCkL?}Zaltc|mX>r^ zuibt1ZM1AC`Pw~TFb{>pK$j?b8w<3o!F2*u6UWOR6*~;1DAO!uW@jJmPI>$d?~#Il z^?e(!o|>$+58CuxPZxSru%>U2c=H}X3<-#ve-Rcb8-vX-QG95z5mJ_W@ z2YXPd+65*ily%WGWqNuV2xuH^T(+F+vzQnZxn5XS)_1lo$iDa0Bd^|cPjLD8jJ~_` z8i+_aXe;_)ALt0|$j*`CHr0rVzQD@Lx{lP}&uMabR@#Dt^2E!yQ5@yD(b^qbRR_2> zF2dzRaEPxvnqDM)tu18lEsZXKj3cm#BK}suf`$w4jG^pSKfg))1OQLFbQEmpNzUZ^ z;uX}rK< zWv*&&*dkb2Q-K_^FK7s>crV_+d-s_LE-xJ{e`j?JH$82+I6M9c0K#k6C;UVr5uh&2 zbx>7Rb$omrsQe>U3ec_u9wnEcz*_6};OHc<(LhV|Fe@*r?0RCeva-BeO9~1i4R3*C zngEo2ES|*bjQ4&Jg0Q@oJTQ$UtmMv&3@pN;5)CcV24cr0F>%t zO3TV_T)+Ohw6r5=yL9+xS7v4=DEi#HcP~TeQ8x}3*LlP`_32np_=X3cZF8vD4mn><9?w!9o z+p$0-hu|!gvLI9gsL0i5Pg*;0h7;^Uv{MQA&e_>{75xuzMxdr_T5fQ=x{*uP>B>hj zDF*~%9PBTD>qG2i76sr-M@w}t%eh=EslOdXlPajovOF$SJXZ3MW5X!gGyihhCs3fy z5P`+|Nf``H*hKK~^5!KadCqso01v>Uj`;#9pe>)NXVkquvMLqnFr zG@!sg!M=b7c&5ZE>K9f059Q5hPwzxWM<)Tv?W~LaDgy>Wc> z=n<^aqydYy=Q;b?clERII5T82?!}8Q683pANIiFVi8Y0SWMH$b8lHNr?}@~>X0>*| z6X+q3D!*^TPp0HX8!*oy5aZ=9Umgcj2HA{^MwzeOb=nW~gFv1EJ4R+_XS)DUN1?*z zyU&aVK_IfmV9PXwcYgf}{q^k;;0k8TRFL6eU_UzaS3zYHRL3HKDg!TJcJCC#Kp$9b zPEHOUkB?>50DpX6Uw@b;^rSU|N`2FJoh&CKlX^izsASrRAH0Jx{stYS-~h<4OXXv- zfwsN9y>Q7~*dEWCF$Cgruz6JoD8aAzI!9cc-WfW}*4>?$+bb&&U1ik@m6S~MyAFZq0;R9?rVK`l zB7qp&L-(cE?nB62GKg0#*N&drnstV&Ts@B95|8Y_$b0Yx1iSVn@Kj+V%ilfWQOze^ zkN(!9J`NeQ{lfC~kF4lXY?fpUu`C7NRi7(hu{X=&*>ZiOU2|6!}?9)ew8zApS$kbjBVpb48=O8`X>bK-v3RWCLHh;%wKpF8?*VzY$uCPwCmawk(Ob7} z)tx)@Ob{>(eCEw6BM2nI<#T!0`s**wj*x;rz=a$-$=Gvef@z`xp6D4Ew@Gzym!HcR zd8!w(y>xbsu0a0B2<>;_W9yeZvIfk%#K?BAzgiJ4qowhg57cW4bExd^JY#y~T>=AE z-NQuvryAg2)wRf4CGW)-(`#Rw@AXp&ArXkXbT1PVB^UlV6FIDn7JEQP0%WYv917Zf z#z)lEJoO)10rs%$jZRDZc#GvyD?qhI;2IyK3If%rfmN~toRYWh47X}>x6)o&_C~n6 z;%>1(B7;GYaGS0$?d0@s+NYtiQ$K9!Ai>L^G?B=w!XN#LS0(#C(T!InIL%QA0@+Cy zHd-pn59g1bB$1dkbf75m=Oz!E@i)N!nc3Nw*4Nts8Rw9+4pr@=hZK~V*E@9qyMaKy zJ>Bd$aPLpL+7ri7^QgFYwmk&i2XL()w9WjCghGZxOiWB-X(ViU2<7wazEYvn@KD~E zm>A(dqpM?;`zS_;y1K2cZR{>ls05VGs;Eyv$^~!SSOy;NcO5+>cxY~puK~oCnn1Od z@mbe|9YCCS09J1H{1jw3p`^qS#A#+_5b0lb&(7^sDGZ+S7u+1t~fY zNOg^kF#v%lgWYXvc z@0uz|10`Q;|IY89-xxh@It#h`VQLEWu$PB9$8>fc52W6}$U3uJxL{yl;JLj-0A6R4 zguOSj0)bqV_FMf7$a@wBGxyG`*mvbqP*9LdJCA8BRis`)nqI!@Az(P*1`vSy+S&21 zX+j{`w}O9t1A!6%L8QHXao1GbHZL7S-RjAcC*XgfSTwo??X1H7^-B;4*G;(|9|x-U zkw1=xAD;lu&iZuo_6a<};h?Z}5*A*XNH_)g;dJzOKRG!$?UTI6tb?QD1BrY2V9IMw z5Qs>n@f((sU*Ml1yj#n%d1pYa?so3nx#$a^)V5-?w73YYOG^7Rq##+&bsXSv@o+d? zT^+J!34yF;O4)YT+V@?Dy{dEt{JI>t1u$Z;gqj)vyy8*d!5BFt#H?CA&(3y$@CQ)g zDppRq>8I{(ug@#nPjV5=;BehM<$b_4&T`4p(a{;SUv?d<@B?1$HkRrV$HoB>XesS^ zUi_^%>sa8E|4506y)p!l1si_}dG!edBEQyHp3t)j0)k%*NC*h9s<25fS#Ruq8?gt* zoN`V_PY=Pv@ctc-O7Q2Akz_xH^#8@)dxk}|McaY}C|L|70Wl#;2FV#jf|4^7Q534k zU?DjlK#(X%kensQB1jfgBqLdJkRVY35s>s|oqO-Q-QVl?b^m(5x_)>(2*s|w*P3(8 zF~^$Y#n#jzEUJNlfxt`4YEsulFy@ClORTwoN+FzN#o~^c+cM!@6%!Jw1^k6o&5VI6 zMAGATh$djQQK(dGLMU_V9$fz`s`%HhGfZjq@em#R03qjE zH4%jUzlhY`VUd73i_F--8ESZ_8g!!((wAy#YKobNZJ>OtY(T{hLTOqNNr-3!zP?p%Oad<{IYPt3wBSS&nr)v4*=f{ z@>Cf|N!*7ZKP$peu?;AkF-}gbxwmfJf*<4>ml%&@*Mqp6@im^)us!OA8fzSOOH&b= z<$a>J&ML9k$!Z%RfO>pEvdx}pgh?BNChnxh`10vOk=S?D$yUGBqLFmxX_|074(o zQXi(K8e#e23)~vtfWxTV>Xi7=!%vyBva^w2&LLz3ZNEtqKgVFMSpE&FIfevCJp9@q zCn}89mpQizXbw!#H#oe3lQNQ!BR1p|6le?aiz&VD?t4CUb9?_{7=DQ&KDS}b#!`0z zaAVcAwdglj5pRdKF}f^q8z$di`a?I*Jt{OKVh$^_*WK3_sS_iktgK9G@y5`)x3_nj zpB`~D@HNz;@NryWK|od=jwV3km@3T0(8mp!rO<|kBX~CY!LXSNJ1~D1K9h~oMWoJ0 z{X>z_1o4>aFbVWUmt*{{+$|Ki|}ymd=G_SR(b@DX;a%R$Td(LKvh?8P>2 zKbGX9j>JG~$3fyL0E@@UWtn03WEk_zjoI0j8?AHGGUwD{Tu06aM??VN%8D`hr`{=y{g>K+%Mlco z6_M=(*U9K|UwjI=T#ODOf*R`a~nvmFAu2$!Hc=R@Q&+v(B^iLBLgtV`KjM3GXe>BIe!Cz#%4;zH%;bKnv`{LIgXbk^GinQYp{@FShCdYFj~yriH7BrWHPG;mwn|4r zpl4?@zg=~mmUe&Ib@uYD{tA3c+x5{BC_y-l`heaF*waY>@-SM{D(!*vW)?J*22_$b z%y)XSZyMzLF+s&H@a4+|lzh|OqjJpXs|@0p%U8_3HTnQnd=@Q~jNa@8{-m7Z{F=sl zGd~LU`V^@8M!UK~omIUT;D%3T0%2>BN~WD_StKnt^c(;gYgP(Ckx1u95=^xZ_q=mh zva=dfGP72BEtWX=`D2~MPj+ z<)Qc=jHIRgVp{GpGrx*I#}C=zDH=TnJd6yGV}wu^;BUF3wenPRV}T@3a=u1zXWs_Z zpC<8du%`FReBF?U)YmWXv+#rH1dJv>lt1^MwW}+F`~Z0mzk2??-=dczt-C@{kYsD?uS4Ur-GJ zGXH$>9O46knp*#IZ}JT>7v;ZI^y}*I#t^+@jyd)*d61KM!-)_Xx3#ujslv~&>lVMA z()%_s0i;@3GATmZJ9ZF&HhW72p{WR-^p^`j#*yEU6lY?}?!dgf9E2yGpke=mlG-$a zeHc$#2--!^+;-y<>tpp5>+c!7^O zFkdqd$}dhg79L@^926rSIp`p#RsxP<2LT4p_1yin4)`x7PT^2LzI_gAtgKiZs&u|Y zPzA&za9g;72D%V|njU|DTgCFtX_%XvL$Uq)?=v7?M5&(PpgUtL0-$|eLxYwg$Rt^E z12yY`8y)7I3#+Xa(!GU*TF@1P35|%bVEmMU3saXgAus@y2#GLJU9mq!WEXUKa#`WvA+0sWm(R@p3bBHw{5x#A&ze-@p>*xypj&-#4s9~37f(BAqqS9UsE z@hn*Q0`)}N?#>PkhrT)>zvl+D>7ew?GI;sGP+A=_*_}Wr>w}yxRdh`EWMyTYOKmpt zs?N8GO}WwE(*sr62X$ftBUrR09-Mr07?{&Sep4`WJy z#19o}*cGxOZ(gS&wX}48W)b+^CN?%5F0puN+y_t#eaDHkzVo$|xLS&EtEh|_D=RC& z+s_|AX3GSTLc0xK|Ix%?fMSeLHgKA``cHRTTiZdB0Q49|b^2svWDb(J5j>)>g($1y zpKBLbn0oKU{keeWhk5C>p`3!3A3HVV0dRCne=v-eTlNt?;-&HrV?HFnF50J+h$kCV z)Rn;MR0537cn5(0d%id<`yUq-t?DpKP{Sg=#hQi0!cg_4qaX}u*m1$NYzyG#f_eBg6onJyN_rKjJh?bfNb9syd3Ql|Oj z&m_e7vwiDTV=R;}d|l)S(R;$e!mO+X85wg>`xf&$GAoY)gdx5Kv6h^aG-lU=@9tfo z9qt#)1W4oHmI9(FHYJ6Fl>!dKqTb5di-d)T-wDjg&Gm+2X(uWb)-xaZ1$-c^$$?Mt z9UL0U4f}|L#mof|l`{=$fGXNj1D*T2n$*diyp+o_F%X)V39oNXH`D@hgihYwKyu)^ z7RWE=ml@XDS02e0x3Pq70~s#9vbW26auJcyX{f8hf|{_TN{FzL`@!ASZ$s-U32 zR>ZIl!~BK^p3n-uIT~$R|KjmqHNXNZkuly}Xg~=ctyDgin|EGeJQO}paMQ5wTd789 z^DlW>*-J7|7RtZ(g8stOcKjI~V9fuREdT~7cJjYIRE;5g`E$Lj1^9#JDS?a66#NAj zOSJ(mLIxop1Mp!nR)?dk0I3vb`};uLy6W;(`VM69=!P1A*z-vY#(cKl%I~Tqm8)h7 z#`0Zy2%VLsQtMb;kE2qR=U)P-S{G_G{XDHnAY`PZm@oTrV}vfg z+|Tg44P7pE?Y`g7Z(Q(0!x1#DDxIDrfv`*eD4jWnLf3-Tdf3{~&@dv;;ZdH4ZmH4g z=4P+x`4H7gW}`;8Xia)r+ApG)I|iz(w3dRx6PmCDL=@_ZzK2V+CcM}-=ku;esjHV) zwGyEYd%{Q*i^QT-qL(({BIwmxfD3&7{Q00(vV^B#)Zyn_Bxlgh97ReSL4{LmzWwFP z7i=U!Ahd#jPQP~l>A8J1DJdzL4p^u<<6_^hQiwtCWhxH%KnwdDwlyHwHoKGv2BCP*^q0lJ`P+$@#9ZwG$idBU1^(b5Rq`t1c$DerD`|0{ zied)=JugM^6zG;D>AynoYi7$|7&m}q)-FmLg5Fu|gbJn|ebotShaW<6Hs!BKL}{dL z-^sO{TkP!XljkJ($~uw8Ca8K~fR_;vIwSU7UK?`Dm0J-e>L%UfK4U%I-DG1K2&WRm z+PjdxAY0IK8srRmzEewAWtH&!!_hYi2^te?r2I)V_0fCJD=*3E+C+)3iU|lx)ERC7 z8xnM}ARD5+diNboO-+E(*+9Fs+j-cGm2xiJzADFH2ipjTG`fnOsdwiv-@RKIiI6XZ z20JX3*E0yc5(5FKgEcFj%8l#k73*sZxHLdZl<)$5h!jB&>)LA|A+eF|F9oC4@xR9W z7JTn8E2+2&b#bUqK{4Rau;i2`s0N_WqGrM~7gYJRl_Znal)^sE2Ek4l* zoBQ#jBWL;wv*INg;M0FW=Ut`)qyQ>V`=8#q>u2w=yBx8>&297^7wd9;a~H5GpF%l< zi;8U2)i1qyx4E(LwyG$AB^Mq%DiuFGb}H+|>9Kp9$vqqRYzzGlp#hyOA3=Nzs3qDK z=hsh8qY35vy1Eo|y*a46$D!dNm*x*$cuhrGE~AiJVK4BmQ{a)?%`_taMz=tz0C1G^ z_%^i-oyGR25Xd)alHLZ3@_=Fl?5Ex`95b`9JUZp-!?kA%ETEwMIg}S=SqG~nRN*^; zUopg2*5l;R1{4rIBs%lYjniHv2BnTxv_r_NSBX4k|2%OI?}LzbM@oq+1!!>) zSuiC_w`Ip}jbnpPeW`A2Eim{~)fa0DzH(F{csBTKaE?#q9FPVHn==qypk{jI%%6#) z26*#e&jWIzs;axKxBULS@B;z_A4Eb<%ytJVaJTs+~B8Ffemz0#Cdy|}@vKmDOelC%Ngp(VKLZiuy z6}4H>G6~u~?rv_wg*qU|K^6xb=Rv&#U8f>wsYCK1xkX1uSHF`V7Z`YMiOhtXP@gdM zsg2DftPWUWkIEiFj-<(q%-8Hj^`)@BAioX>?Qx8eYpIE->sAVvk2s)OHsmj;u%0nq z`mOf#=>h}>vJ&_m##b_}bh?m^j^)IU-5#(?KpCBR+{1)ZBu!LV4Mw|cPJ9*|@r9hQ z#N|9%Cw=xP{l~%ntS_?^XmG+PrQQKf0vi}I9Lv!&8+ za%1D6bDXeP2m$bRASg^Qwg2E4bm(~8Mzj0_4ZhqLK%K&^j2$d4=75NIMo!NkGT9E2 z8Yy5}A&ERcvmGlm}C-pl}`E#Huct3MwPxSyRi2~z6uS}<)`Bk&@wOH9nn_`r1o z#RG+S>>wQh{*wd0@~cq>@b@i+mi)MYw5Y$^`w(&}32qQb>r2xy{NC zBtEpPy;QD=St9<~-P^NZe4v1`xV|g+3F|dRWzSGOYYB z?E%C9yyy_6x=3+DGJMdsF-g#_%NU*g2nXlGZ-&qdSOY8Jq;%Nuv-ugLf+4h3zxDE) zKtuuZJ}pM(t{tqo;@yALe>ai#q>9nbCPf$k2G!;;A7sv{Hr$~K zenth*MOzgzfWjw~6^eucxW3x<8+GpyN)yd@znd9)66lYkEH6QcX){$4Ehp-<91aL% zy&Ul&BsZ6No!HuPw)q6EnI%^t@|yepmNBf~$msF;d^Y zeKTtDNo5Fn#2j!jHz&uj zMau2hO$RI;Sh^ZGa)F?=~6}FnrC0ml|ZNBvJBJ) zppD4xN;`mXg;XH?`bb9RpAAw*T3S7vY;u5VLFPIfKLUIdTGI)WDAg6AH(;|!6p$wF z-gW{)1_zl(z&giLRl>d05HFRZ zGAijW)N6r|cGs5+?mb|X^b!^n)UChq&4do?w4w?MuiLtMooZ0h1A&CPs@H8tneA)s{2(*bPkfZJj5=3g@&1bN(qwo%HvT#!N)G#~!rGq?xXE`Mep1NqDd}Hev$S80 zwfOo*WAERqy}^H-;Pn^sLD(SYFMI{wXBV;8cIFhsB9JM`X7BLx_mNx-T^(1{Cuw!2#kuYpv-KHCoUnp;v{yqm|G2IjI z5?g}`1DFCRDYij0W=>(A~oW{@w)nub3cIaYy*+U4~ug?~P1{=P4)Nmc6oi z7vdb7EG1h_X;$0My64|t3X+r>)!oj!@Rc3lCP*!cwaT>t8Xp>lvG^R!6FCcghlNT? z^)mSsqQ=6X<9N4tqglACeW z_tC7qA}emZMs5{Iao1F}lOn};J)o43OXGo3O$Pui^4@vI+vtGH8t}Mu#y?HeI&h=d zthxfWcQRCTJmy+*Kb~+)Ntu?^7DSr)3p>vF5&D-IH?*^rM|%PU+y4lFqy?8kBj_u7 z9}-toRm9e&6JET|(RqoQni;d60$nZDB>uZ!e}Xa$P{DjN`r_r}q8ZKgK)oES{U6+8mqNo^@eIxu%l zJ^It*KJGfctX9xdS=Qv>7QX2FXu#OfY-Ep8?wJCaFL?9jhoPZ2>WUzE$R6*?2_)W? zCbk+Wr2g4SdF9F{l78B{hksr*-mKb(nk0f2X<~AH{Xjkscq~YCTFp?R8p(Ep1nf;> zB5Owl-*frv!rp)Re|FweetS~})2y53yu#R_= zy#R3H%f;EPLpm~0#eztsuqQTN=nZJ9l+-xc*~3V&U+#HQXc)`n;Clnt91gn8wC=*Q$Gd7XO#UX0&d$?yPWP7fWiCXoJQ&b28r~Vv zGbM>uiF*fkwq*Pht`MePQ9bexh7HN?+67f;RaMpE$smsoZ8HAxOFZxAyUQvJ3=Cb* zL_mK}>njNA@-pNE$Ps^CM*%Riw)locM-i6fA27bjK2cUls#~*&zyccB#v$4n zM-V$=fCnZXOD&9Z`tdbWTyS%3O_6oPHp35Y0GTK@Hg5@+g&yRYhD1e05vH<&tmT$0 z+gESxqb9yo`uh!#DgrMlt^&4({>&@@ZV(Le+l+?qzx=Z`ZP?&q2NZ@~y(;#dPWq|P ztk=ix9vCyvXKAe1K>@XrPVHhsI4Jh>yff@ctnIT)XQ*tfUM&<`UeR6%Q8?IM z01sG-+q+jc37+1Wa!KC%Td_nkYnzXe#)F#HIWnjQg40;`$` z>MZ0s#6uaOTY+D{euaAja(C8?$2j&6i_%cc&dvsnDU@T>p!QQ-1)5Pi9Vc~IAFxv! zWa@F<__n<g^f_9>vSs8@WOY$4cUW;r<1CrEm{nBF7(@=|Ktn zm3Ut>id16qAp=mjz_W_EY!FjW_-rqn7dipiC@8$)RayB|v`mhSnAmM^jXIbA{(X(^ znfDOaq=gb!SwY0TdaB_p1~V1RQd9h4#{fxr@tUxDclGO!-ZA;z(4H+p->8`IP6LrU zfK|8Ee3|zb-a|O5Kl8g}zh4e@J~TxGUqg+#xVZSKy?tlq9)L6qoO`Ew1{i}reICTX zR?seG$%VxM2Z9O4u0{WqHU?#JZLMo&!|WaIH_jjc0!Evw5iD(MD_sPwH;v3Je2tzfWfu9D}7@b zK%;cPms<*|beTuEer=8I>cRr<_|At9($y5l1vZdg2Jzn80U>T-vGo1>_l*r@;_*gz z?9R7`!4$Btk3udp0c8x`6NeE)5V~{k2tX5q?RPRn^z(#-xE6JSqa@`8GY}op1NVqF^DwhEM@Y9m9k^gV`|JL+nK(G3n{P@bRE#i@YWR3k%a8P8TB; z6db&=xCq)W3-~u+dKO&pbGUiZ;Hc8;TMr?eU?3#1eaJS910Esu#SPtj^~??*lS*@v z!?UYzp;DA>rOXZGLc_vgQixo1qE@#BT??d0=D=h4#Zwo7+JWOf$G~|)g|5b>W^0Rs zSo<5}PEh{Lr}?*f+YpnG&SsOxNJF zHeD&j6*!myK3<0i_TmMu67TKl zA-)Bj<(yc68?TZyUnV4sK`en=!_9`xW7UH$6o{0ueWs0Ww_KChyvM^|y}G=BY%eV> zh4;`~|4RVl61G831q+7x>`b&b0X|`q_mOV0Z3^fqfnJ7rD?4XLG>ciW5UvE+*x0zZ zZ~{?!eo$=}fF3561F7kxYH4Xbzio9>M8qAQ9^5E^`A*KxdpkRzb<}{5 z0lHH@pi~qK*ch@4^YBAyxwkeqOHN2XT@c@$Yb7W_8Mi|h4vsT?dDaVUGu`LS{@<6KQCXtEId_au)8MV@d9=@ES@2h zb2^T3&(3!QVay%6ziO4~$N#d|pcmtCPu%#ET>&T0;CKlwZ2~A%>g($jZ#XT=m#-*; zC=hraK=-&~K5)3fz+eIhl!LAh+!3^v5Fg)vU(hCMHuKLIZU`5Ig~!Pmaj9+kp%&!5 z@#+;Im8R>(&3cP&8twZ`t^EZ5@Xv@) zY5Vw`LgDqk>DJc74H@EEGfG@4tcGJd!^7H&tU@6H(g`Z7As`a2jLe{}Zl+OA71il> z63&QH&ZkH}Pb-5(^#X7+fhPbr0ah*)BbBP>*gAW9UVs8>DXzUnF(7atj_<1}q&N`# z;iPeI?y0?>g)tLGO&%j4)`O%3Q4a)d>)$GvCSQNJ)J69r5vo=oFGDZw!19HNT^>mw zdJcE8GSF&(8$`84|7VU>{_jt8o)Q^|#$3#pxxy|BK0~xbS4xY($;@_5b)r#Q4+|UJ^g$|M{K&cbiiRrfhL4h7BLy|NQPd zbh;ZP_5c0d|L;bpQtucfZRjjD1^)9LMn4-bll=JacliHgb?Q~sA#qy#T<6|@zDdJ; zg*?d>mH&Cs{~NQ@s7Lb>%U8_}L0SFp>&%fxkwoaff9!u}cj_x2ExrHyjG6I2?M{n= z)Hz7}+WG(MZbmthhCBX08lLVs5dW9q>6*j;dSv*2F+3Hj_^%_wlLic(m0VzWdRXilQkdZ$NSGu_5c6p|Ii)yFjMl! z8}VVlpk(IP8BcQW@@Cq#CsBcXbJQDeKbl+{a|wTXg#qVKGsjHxbacx1FStw@{W2gn zucUVj3FLdCH~m7#8+(#y+>9;m)gN(&3f&a{q>ohBM}t}m8#&!(X8jrOAuOqcrn&L#|~lgUhq zMywNvrdLw<5}IhPzc~9QIdUf#kASkzk2BHsR-p5lFgod?bxDwPO~y&0m7$TL8+kcx z8=aSt7^*VpT#ZSYt1l1ja6JwFIB-6rL?hqr)hBo{k)j!5x3 zBpvT-QnQ?e&giiXL%XNyL<>whG6C-cWh8X4kX?#r}3c#G9M($&<~Zn=#C zJM62WiVDTn^%u>QsRw{1_$s4+7~Arbp?c0aHDqK-Br;2Xp($5g19^| zv4E7VfU~z+68Q@Taf1^E6K(1fI7?J8)?um;h?nE(a(Y7MzwuQ5G_kNQG2wib65g_W z??lb=FnN0Q<&vcRvNd#clOlasnUwEoqx&$bn%G(+G0b3MN`0aogEvdFdkts?T^he zlrrSX0X@|VCIcOrFmo6Sc(SYEZgjr7<~wT^-cM=afQKXPl{*4?S^3YL`H&Rma^igz zFmBY22+J8)LgPTlqsL?8j2k{=;ID#xdNZ#I=70W$p=Rp}I_lkw{admIhj~Rz#U}V4 zu-{h5DF%&NIUFJeRO8?S-X2?71gLZv1+Iz=-G%OFb;_foow~iYMtWl1)#^$`ng$Jw z2Gj#tfi#C4&_P|w(5>3+6OXJhD1i}O6mHc9o#+o_Fb4={S~}ed!-(jrE4ZT;o(sh< z#l3B)<<5JjTJ!9#4QEbOkimpe12(5fJS3jU;4wWjmmL&5lJvG}H&Z>s%mZDeiuR3_I`(%uYaO^fKs zg)NMajEq3*vYwSwEoX_EKqOu6mI}!QH<+MH@ey+ao66Z;!7kpEYhUlCOZfAZsRGzxe~iHNV0q zU^%IawHRm(R|c?ND+u!pk*}Utq}saDj@PLE5rfFkw``Q;j#Tt~B&1&Ur@Vwd)lTvV z2m^30qaJvA4Xpp1xuYD2@vs3n^;o+90+H6Krl zP=me`b*^9CV063 z>vP+{y#{l2X|}7(-Xadc1JnL5b5`+AN}`0~sWteD{=VIOXECiWD8Sfb6>N<>k}>(S zLt9mG_bl$(FxdCfD;SoIfq&cOVMW=7mKG^Ez24AZL#v}b*|yZHkR^~Z+pj4DrlvbT zYjEzBmN30Mg)LD>%;M9WZ|4POu0{PnP^!w}#xI(OG^TEhx}U8kCT>`0Y@7#glP%v= z8^7gGG!N{s<8VfH(aXitGg-@^+TFqVMowxt2%*LQt%{EJU~;s!B-Y+c`XUdGByoSb-X zG1eZnJMRZ~J7C7`nbyFD)RScKy3GsY*_8mr^j3@8=fahax!*V<+2LwHf}sIUyvZ=WHSctiVp*oZ_LuT#90&_jt@$~9IS-3vtM+;r zD+*&rn;v>Kt`v>18Ke{?>ZI(@foT}eKXvB1y{8HOd73=k6~$f=+h^8!;dPSBMjT?1 z%A*G053D9nQu5ar261k3CIn1OiKQE+EX{^;R#U z8#aD6WmFU@)L{F1;N<&&zN%GGUXj>$n{(ALM^xhq*|iPDxKt>si#bZnmPlH=u6c>C?V!7e_b%W4$1N_M!H|So&!^M_dw! zJJ7$PVSc(08Q95Kw9lQT13d|14e(UmsH4W7kp-L$d2t-wOf0T^CJUofVWEvU6eh`a z!!k~RLDRuzr?^3@8()lgyFV&}@9g=12lWPGub%y^;T^Lgz}k9DJnALdjV?+1BkCj@ zmJYMfVf?4AXv1aB`gHdP@`qs~_dEV|n!Wzw4UGQmD&lo&lv>s z8oGn3Y9<-T&S!M;FM9`~I{nH|pe8T^p+pLl>M+8!#U&xR z-d1-OPj}ju#uYsIQ?F}+Qr17$4a>Gu-d)I1??i$qWQsnP)=Qz|dFU{~CHZaYRd} zsv|8F*C0~Pbgr$wcH2n0Y{S=`0*Qb*#rX}gS?_z>XA=%jEhpxx)#s#k)|kWLz?#+7 zGOOL;SVc>ij%<+f^N`oS2?jG>sHTqH8*1Q-JLxG*Lf*kjogOfOTPF&&j=G0^%++jFmc}PXHRQ=4mY|)TW_f*nXa=K<~i>p@?ng* z{c=uRozyhV%G#=mNZt&0?nxhNAa>ND3+*@sx0w<=S71WAD@&P65@!X+Y$~fYifp>) zajWfYAr0C~>TmLhnyvu$7YIC2-3w7pK~>-FKgVOe`hlqW+D(y&E9piWxU^*Od%Y zk0U^`2iNot^$N_aTj2e~jq^l!ZpOh;ZrD5$}k4n-#E>|2K7 zY)GHlpZ|aXBi@(Qo+rg6DbxFlr^~US=oq?d%rGw76igv6^tue=yfc&QobBG33%)nr z8NX}i05H2j!A^FaBK?)6qcYq;MzU`Uw2@g`rU_%s zhzvPJKCo;Bmn*$;9&23z{z-G1=^jScq_`f64j7YZD+n{ZoyhRF=r0cx9`cFV;0E~3 z#l>+@cNW)IyOW^hW(1h-@d;p`@h0>uhD6;63yFU-zL{Ovoye4wD)8?J@rHT}EEw8* zqGSIZC)v}5!VdI;J+6I}YRtI*`ZAU+5sCtxZg3$D3mJWha}$+>nfCF3dMc5LKe#V| zQ)(9>7Gh#kHh8S^*B^3(qI2s2KhPBM7A=ee<4ALRf?s9H-<JbO}G;TxMNMlEegbCehlR#r*tA`&fDbLuA@v ztG0394+ocf{WyiaB-`Z@J9`5q*dCW8m_mb7MO4tRu;lthe_Y-_JjPY65ZZ0S+GU=T%9Y)w z8CW0`>|Kk&3q&dMp%=R!Y`SUZZq#!nKXE%LzM9gn*X}KCSoRP|{*W<<8+MW` zQpAd+#LZ)g(bEg#X(pF`U4a{C$k5kSg$#s=>=HIX%XlIo&m@aw3?`>~ymu>=x&gY; zSvS_zvF44BWKT_T_~x>ut=hKXAsx;xgHuV}+Vqv1F-Q$&tabNz0qdKM^<0VW*4$14 z(TTNUeTh(nxj>*LNRpvf-OSi%|=a9T`k4ZaCDx{3fN#|_JQMPwVYPBT8Frn?<=zb=)#J;oTWfgX^HR+@n zDTzg=GS4s`oE=K}9aDN*k0g=AC5Dv#xZa&NQOggFyNP^rH&@4W?Y?6l1Y!%*ee}~w zajO80+^Z2{_XqcQA_GB&!*4SkZ@KFX#ly@;Ij5GV4eOJyyQW{mS;HB=Q1G2ldiEIA zy-$#z=3-iA+2QtBsNIB9R7AqfD?Gvc; zF--Y!w@BpJKa;2#3E^`zdl=EGCWA9=ZLsGpy6dx$rO!E797^*%VrreNFv3OL^ZrKi z>o?vbQ?~sdx4tFz{p%hSugx4hHPl0Zc3MY_^0^LP>nbPub1Lu2KB0Rs+D~HV&{W^Q z+-yQt4!*H}C(@P+YgMsEm%EqBQ}+tELJ8hgWm2HdeQds^hW~NO4 z;mBIrVd%GPAU9)bZFWT7E;L+LO#0g}3<4@FV1JlAt~)to>w2G9BmZrY4!IlWwO_9Z z?=OpKrK&o8#^$H|$2Vm>E6qb*NO5j>1U2}UDNkSS^Ask!AqzC#^69te zcG-lNH6pBD?-zz1j94L^C!JJkjVUjw33n(Fz*-!xH_oxNiC22UrBPMeMOUYIe`5RB zQynKDdytg>dO0KOBnMPYf-zwq?Y(cOp<_ZRMlD<37y8sVV&bDNxKmZCpt74ppIqt7 zU@tY(IkT)2zVFs+i(&RDxI9y<6|=V>^_xVeOY%-vkbt8u|7T_VYllxKL);q}Kk&Z^Y? zHI81s=%OzrHs614?a%pxyX$F@-F|V}%qJdRq|WR7EzIa>@7Yk3?B^{4cNJ-_WX`goK)tb&?wWMw!iDiR6%M}vVzheE+9d$NMw{}biLfIZ ztDUc)ZL24{%k$VT1@gU(+eI49*Q;a|Qq5O>X|S41vSVuAY;#RtIk|ZK2dsv&`r+QV z+dDJooJ>U|2~JT({Z|etgL0D`^}%PTOb}HahkaN#jCy61Zk;!V)$4|bg zQ311+;UJP;*Kz;dv4&YDtrVR=fIrR+t``lk*KWHfF3P(g6m2)&(@6q#-jhy@CyNiN zg@{CM?|GJufzjzsX)sw2yh%2jAE(jc@U5$|``Rf=s6PFx;Bex-DNvl;f6&n`PUKRY z7~%fw10V%mGXI{Xi0JAYoI3_Oqc6?-xxVAHvTc zz5=EwjdN11!YTPUu=36R%_r9#-|3;}tueyS8df-Qu_6u4m9>o*;f`PB&*58TNy)%N z4;FaedM5WXTYZ;0;+*u4++up3*n*?)=G#uHTUQlE?4a(4gJ0yh>tx(Kue=@paG<@? z?V1$NKJDx_k3SOZT|4#Og9^<}T_F6O2}y5mjf4^%^Z2 z9e&HM>q4G>&-qLO4UD;0K`)md4sJR`Iu$UiJaQx29RKE%W!XRJeM5i;_KSgzK|Ytm zh;!-D*zC#c`Ieo>p8wW7vzxZ<7k_T6<2L8tiQBTy(R`nFTQ3|iIH*lr8xV}QYx9RZ zlSDICwq%VgeOPr!KJ(>W?Uh`H?um=khcBpyadylvUT`}1feQiskgID#nG+d1 zf$-DpllB$sCbCXt#Q-xvYonYz^fuMsDYn^v;CP`a~}s&+rNECHbnj+k~RRNq$2A8Vny#%-=;RV)}up=}KX zo!Z-QH`P+C9CEX5dVFoLJW+WNL zqT4r!bN6&|R&>s{!C*colchoL=E6|6p(S{^;f0;}}Z`kxO8y9%x~#zi4hqzF3R|kj##)2bQ*$@-LWDJbZfy~k?chTsoj{=;d;lcaJU2ZvquVX` zrE%hd+oMMP+G~@4_??7{NSnxXOrBiAd1;OW6{>0L^c9aE3%4IG)nq@Lt^^CWi&m41 zRykw%@Uis=`?h0?rlHK5=v$t$aQLa$+H?cqkT*HRanD)=7v7(ZdH?}_YLjXB8)|W8 zYzKO}Nuu}nwe$g>=QN3hqZ24+2!lmlwLm$gcZVLQ))*wUu?|~(Kv~`zu_NYn=ix}e0toO$)VqAfvv%y z=Ey0)jh^vvt=iqArOnc&eXq<2pEmw{_|SODOA&G(t-a1$7iE{zOM6RHwj&yg)3$ATB1MXm2zX9ESC9=m29*Q zYkS?bNT)jl^t-(USA3g`^iQiC!F7&vr1)L0*Au^1$YL5gw9pXy9fgs!x!1pgz2GJS zz)PNIcyjRN(VA`eMe@Tyj+LCWyVqFe$x%y+4{J-1rT>5nA=m(}-d>WCt4-jKHM@v2 z{CrqwG{Q9qC^~2a0Vbl5!@k{i{b!;NRYyZ(D3<(;uzDa89o6p3etWxCq8AFaKeo26 z!dl`15Cp!Cxd^U=qMvIRnbe5vP(NRM<|#zbOrn$IkcTU*sh$^!+BXR8n3!HFrvwP5 zy~+kAM;a|C_KgIubWboLkKaBQ+ifC ziL7HWF*eo#$G(RX+e)-dHzt$mcY|EJQPEWse>_rxr{C~Mbs*_nbIiAwM zRr6i)J8l0vji@!1h_jf(W3+sZKmk7Z33U4I)=WbU*yh>ZMTh z(5=8YJ^qcb&G=pnuZ>MDxWXO1(m_v`m(fG9p|5E?n1xOJ;+$hQ|%v+p3M=*ZT?`&Y|LrIH0e0HW#Op*8T<6tN6$opu2k6n0(LHO zJteeh6l*#JxMGQhXXl4nS_jwp@U0(eGktUFq3F+ejB8x&zMTD>$uL+MT4Kw6^-M<> z+7!S+RuyxYY51`tp)=o%!{MBP@k<4lOGI@WH+dCDe5Z14@9x z#SaSDsz)^a?0Lm9YVjq7aXznF+K*91BZXb)rKh0j!i{RX(g!0iUVy6gUyC)5CNAIm zrAdb(7umPNeXEF)^7Krg3lUs;`cn7a4pTSJVvP>RSO41;s+0VZjUJI+gmPA)-7-HVHhMvttQ-P zBK?3h(xC2#QOeH1bm{f7rKF~(CwKwP%alUBZOdGG;p&ST`n&=6_lMbx)?is?`!<0A z3_am;fLwE(P#!0z08^!gP-1G|@P+Zt&c`#qINES3UPUx^5zFhYd)uIn7A_}^qNXl! zv*NW2|NOI|j&%gkLzs2(eg0>xIKMg~pd49f`wQFJ;U~Gc;ZQ?i{Q1Sc6vAq^bMWO=*eu=+Hw^$4#iIAQl_ z>UIyZTg#WU>`oax0E9tF8L}0NRXjTsODwa8-i@?Y>hy z@yGV>lZtS0Yp^iFaP=}{Z*P%)ZF=(#?g~Euw-Q5f=)nNuSl{~#s}5R~@Y5?RE3|6B zBsvQ!QL(#HynV;o-j`Fw^gr5r({QZc_3isYhEPc;NkYmjQ<;hskD)Uf8rV=tHlFZNkZ9VIGb3bqId%Iu!*Ne4n>-Sr(;d`CuaUREh z9OwS&2#FpFOSpM0rr}9a>4mCPt%CO?xn1aYyxG(1o1WqQ=0QtW7v+Q1r(4@U#c!AN zF1GKGud_XtuILcaN#*WBH;tkDwdE1zQb)t>(My6HweZw1{9bn-`(3eq(qjFdV8~2KK}Mr2?>d{S-}G`<(_T&`)~4Nr!t9U zF#`|q!=<*fBu4L3Pb}A=N6jJ4J!zCaMgHwLLp@y91Fx84L-!7nLTBq~&I>+1CN&hf z%OAzRh&H!xOwgL?%KITs{aNYd7L{^_#0|CvU0J8;ccr~LOaLK=4e7!%CwyKWLL;=S z?9L)3#LmE0IP)q<&#Uvm7JUJ>D9~UP9vHs3)l6P57TE}l zYXq%kyZ-9AGdptgyk@7xM5*oOW5k}?2t=)*RY;5j$M3H%Z!)1L{Zq0ZoUnV?TsN(g z{=fiI6-8yrC4zpYD`|XBfisWN)Vk>&t5nB-->1K*rAKrj>-T}oMSwp+U4G7kzwn7} zo~fE^l6cS$k!*C->zt=IWK9ULF|r<9(VwbliC_QBD^<6GLLppidRj%{Gf@&hp_qBW zu=i3yuH3*t;fJg@>3_6I-Y1}K8(0_QZ%izs5Liz<8K#L1odqg*9-?N_lOnFgC zd2dU9HYp}*q#@(7Ot{sbwVmO&hHEx73H%4~89g;kadsYFq@eWm|tJ}84Oqu zA<=)ybBsJy@Ursgzy{<;hLC9k?VQM&W44n82A-Gd4(lA{sta+msh2HiYpFnZoBEjS zNXJ9qHq(XU39jDg#Z|}Nz>Fb|rqArT(=hYoOf@ zfgN;);`ujAGafzsn3$OTDQ|fYBMu+pGAv4XEB-;$nMNnWS|DoZR6%hPC0y+!A=-z=;4f+LH>aT`fn$P~adtC2qSaaM1}LUVe{?&D)tL8(j_ zhcYEoPL@Kl{b+Nc=Zn||a@w7uDD=j06UDZO8rtVX+^Ox+K)WvRbPs0I^&tmyIrqr9 zgiY`3FY_;almoFt;}@@t{QfNy|3;i^$za7HHA*qsXesI~qaSPi>RgWR%DDO3;(XQG zX#r`^1GhJ!FnW=1&jBK+^Yhk_HKU{q?tK_Gk4YVcZC!LK7>$NEYt1{qT{-FW%F$<; z-Nk1qsz|mVl<1||+jYBe5b|oDj^{-cfddMeRpp7xDx-IXpS69v!o!r^(Ad6YQ`;Mb zwtdE_f(#qOmA>mw(Y5;09o=l$_?)wM{OVd!x&8)K5ou596NSRP4O6j{ zj}7=U>w*fUg3F}pQ`lZ{SP0?!pT#XLOd0q-+OPN7alds;q5Vbnk!1fr+X3MupFgxe zMoUA#wGk4sJ(dHKRv5rRe1I$e{ z&iyEnWbD~uu_g!hk=;)Ynse2sE?jaeYIM}FVe+)i=&-r;h!WoHw+V^PJeY17cYmdJ ziDa6s{!H+}WpHKX`&E~5!>7yf+aH==CtsXy(l}~bc=^=?oF1o6?Jhv4AhdzpC|YMD zX=NBn;Uvd=Y1Wjwc^`tUXVHV^3>lX71C^&-Rs0J|(jr++z5@-qUToV@(_-0(qEh+$ zYmru$yuG73Z6b_0vK#W-CmI_YYqSeOIj=xCsCZF&rNb$pYm|Kn_vs84@=5qu)T4Q(o)dUr`3`47wUa7=UAQtVJIMcd}Z{q8V` zyt@me0yV4{3et)k?m}ba$|0HOzRs`}flrm+bPRx`9uYpKx*zRwTG}YU`-s0b^ zfX(V~Ot0pH%h4|mmfpckP6R(ODP>Ujw^Oaq_HF*yRy6*<^Ac#{8x&cqj$V61DH`PM zF(6BOdUhh-4~fnkRw8y6ksyeR*O?PX(CMZbrP6xN+rK3FyQlITpIaDCN^%__36e!c z0Ei2ZdgMN3qaDh*PYd$UQjxuL$eGg_}UUdSm@bAO|Yl#x@!<`?y%v<06V(EIyb zn>8WH<=G&r7}cx%^6A3vD5@ZQ#siVf9MR`XOF0)rkCH{inm|@ZKf~@~M55tp>yB2} z^@XwiWm}zlHc&Id^>G{k`;&ftUpEQF@}ivV?FM3{DoAoO_eQ^PbBO+ zn?T|F4a9lBz!u}kK6D=+^zvoX8Op^OK~MSmPAwBbon8xzv6Mplio}k4a!gTLqt=(t~j@uXvpOJGMu!KUsD<(>8f7#mt2a(&tu`B<%d^Z%6M;=&S zu3&!h_Ke-;R=$YF#wQ+n-5-%-nis!v} zMO;dKatN1L^2mvzY=WD`j4|l5+m9D&r5%Pr$-RwjTp+d9$XP! zk0!Wqgc?slcde$#7B7yk%1zOu_O@FuyzojarpU3nZHnQ`2r+Hx$e*^&zlZ1ts zoTVmFU(7UWNpo3!@5-Or{yES(Pk2@;Xpc>hHO8d0z!KhIBOv{6%B)D^#5&g>qY*6` zn+t`F!>lk#h>O$Y_1R;th6aTWr7GI??UqMlqE$^YxjN5`H`=d#rpvcaz1fu=5uUzg zw(qVn$5+vv$BINF-l?+}D0wg6HHoT5Ll5v!=^GK5hE8(h#ox(t7%hxEh8sMmZt*91RZUPObLB9vY9b-OVAnQ^gjMs|kc*<1a)u$3al{dw_R=I&^o ziO>(OagO2ehiy;p zQ>Y%H zJOhb`KAp>xFTucofZXzCuPb7G&lmN!bH@CIdSR*J;Yr;SrV-;OokZ)8%*#J|J8(M0 zI;AJvMK#(b<$;LR(EU%R$(A>RKbMJ0)Mllw3P=wez_4$1M0sp)AV)wU5j4yrSU5h2 z-rZ)TR#@CwTGYFcb){oZL56-uu>{Y~oqUgmdL>jMTlMQF=U>P-ayt#yJDCLA^@c^%+Zz~W z?i@)7sw!6sF*ImRjk?kiCl`jLEO((>eqdd;r$R`*GW&8A`<^?Epb^HT~-opTK zLs@}cyJmSo^ojE#4o+t0i>wPyN5@=z(0i&NppL@TYEs%@n3AV-a=+2uE%kn$$V7kP zP}VEB?eyjMgiVofVbT3rbBb{Vg}blNSpKq|^TVi2yWV#ORZ?3&>^~W~9Ry`dfjzmK zF!0ajbTlS@xFCJh>yl}Lftl&w?n}~N;)g?|ZF^Zw!asGdGn^CoBKT+e_o2s3vjj)f z8GRY!jD0dj$76;?95M@K)w#O6f972J{9ymMc02~l#WIR`v&fd4UcTJd%|{7;4Pa>; zGWf${Tia3>MhnU86=l8?(LNut9@pC%)JE5znfz^R9j#$WwytKtYq$I(i5h6fTw-wQ z73`?Z-=AP`DqYNdGHuE8%Sb|4o!=`4+iID{XN9BoyiQ}iHOrJu86)Nlwt9N~J{>m$ zawe5Ou5Pw4$JL9V{bFLfGn~IiFLhtdwRrxfsO4dg3D@~8k2UFn)#Ydw%NEkJ)0|dp z=Bvb)x3r)0H)%30_Y}64BKjt|y#z^AS^SDyw1#%Q_xZ@3DJpb3Aa8obmmZz#`PsbE znZ0?8^6JDcnvKZLmm{*pC|A?PvC#|6OeOXK|9@nx8na zh4aLF-C)k>I_>S{dH#~jL84lJLa#dV@@iO}5LnW#)D=7t-_P^zr}os#-nrPNp36P} z+Id`#pYJ*pHz&_Z7NR2UyQXGLPa+M_10otokVOCQ_^%(PW@bBhcgxBaVU~_CAAR<0 z1lPRxYDL4@Ha0G%QQ^Htr1X^*1ZOM3^n7{qTo@E&->+U0`^#-^AT@*tf&Q z$uw1y^LiPsqyu-R%$p1x8BC$r(zc%ZhYVHnmQqq$h5-@`#o&o#u(lP(n35O2y}`>U z`TR};838w*i*s{f34p8ouVAtWV%+PR8sHFgaJwB&C|+J(xaa(-we=U?3@jKFn=l^u zcyJlk51_;~#UC7>m{0;SBj~6uUK{u2jf^^RF$H(@VTo)yf5r;zAi^$X@g%UJC^KjK z?i&9xYl78C5&S)f+w3VRU0=Sy$-AMc$$6mqy3!@yG(1Hhjp7wOpX1(kR>8KHKo9K} z5Fji7#yRkS^>lR!gb|(VvhvL?;In3DW&Isf>v0#e zg-fB4v-1=z5x|QiGvpN%l;b)(z&rC*AZzg)y7B{nMEC##fff7kVJl2>@Xg-7dBY4h zzT1CIG&D6ub(7!2Q5)ZRa&S=EV@%&(f7BDkmVks4d)U|Bso$IQ|(}OmH;O)eVJ_;uo*6 zhyu6+;GX&4#YGN@GfG&sd-oo@+XUBSLwH5QeBQ*w1hAA=~Ds(<9 zLEk?UT`-F#to-4$8DiQAQwuP3h2$iO&1qKMBru)3rg5cECxmQ!*0c$ zvJ*gv|93~=;lKcEKg=smn^=m7i8aD?Y3;7cRW#(VH`?r;ki!I%t}S61!9k|Yz4E2^ z>>e1!t(+mOS#|V!*7IPLq;cfl0L-~kMq`u0fcgYWp}V^q7t}VZ!aVIZ@Bs}8$_a2F zRg0-{PAFdyvEC@q+@KK_8@+TDzGR=^dh#juIqZ6HH{|@sRjl4~mX@$+)ZZT;8cO)D zs0Lno;Zul*vw#1k51|6EA_1%6fUNB4TnnnIcPCXZT)aprc8pSxR+g658y+ws^F3ry zMQfFT7kDbOe5;^_pLG!h_=GU9F-T}nZ-FNb11}7brV#wG*kMA}P~T<9y92l>Y-IRG zgqrU|V7LmFm!V+?=tqQkn?rlmb$S32KC^|9l7@D4&W4CJyJ`T&7;ikAHc$@?W zGm}T*d@Y!wyz*NQ6q{5nExZ%}PHpm{rcKOWLvyi*e~;jiQA|y@sKQ?I&5DIz<9Nh~h=I>@Gx7L55BUi=B%E`RW+8g$-z{FlR0D@U&!F|XTdq88kn6u3B)GnZ4e7|kiUqDiZU2A;mFiv@CLLkOh|9r z2>w+I%x_BberhrB_WV#)yTjJ?!i8th41BvW*gIhtF?={s^5@nDE_hqlDQqM8*$Aef zq1zvv2Bzd)Y`MIRKq*A6(H(xcDsySV_#dwDkCy!48Py{UlG9s_Nu>An@88S8lRTx6 zf%mAwQBND2J@r6{L<54dW&Y1!KK}=>%tX|PB&eoNgm3-j*mMToH!z0(`sF)SFo9x* zJzKp3?GICISUfs#N|BubMmp_xIdW49Tc-Qj*Kvz$K0So~1vHL}S#F&7zyr$-NsNY3XiYQDruDq=Pfo_ZS-iKdx(X-3 ziON`bdEN04AD}<5dJ?BN+o}tq`fo0JmE%&KR2O7b!rcr<*U9Nt(MzJNxR{? ziz9yRJ63L1jh_l!c+?pl9^c(|vhvD8jrc}C3((-1RZt-e84ixc)*Mme)A%dNxL!S8eXk zCXg$@fEfoQ(56XlxKLcGkk`a;?i75DV3LX`okn?&gc2N>gBI1Kj#sbbO3KPI%+gD4 zvru{I;0HnMOE()4hpiL*DsetXKYpAW!mO_T!N(f3KLzf3PLsWD$9FRad zBn>nZS4=0L`3D5#HBpli#R6`w06(Oz{$x>vR3nn+HN>OAsUqK_&P2Ci$u99wn6$+v zStc*I{?JXrHyVi2Ruy@gTC6kzuXtk>$Du*O+3FI)SKect;KXoqLeBK+U!;=-A$01N zNN?nvsOXYlosVaP&z&4b>?hrBIKp-$Rm)vCqP&it-dD_4(R0d*dp$`-`d=9P$CwdJ zF|Vg!krLRF7Q&CaE>hy8(A=Ltt2u;`_`mHf?p!4*EPO|?5d?K6@@S_>O!vJiMpE^< z`4X7Me}G>m120Mk0sqG^@RMpn97;k?u|*YG-{8oICqCZ$_wSMC<}5jqNa}TsjrRyD z@0qhD)Fd&W8gb2D0r+j8c;Ld7TUclnmzH}(o~Gqi%~G>u>}FFhn=LAqL`C&iK0c3% zbXHzbarTP^_&lLTlE1!eCHeWmdORmN7Oadw4;7F&1=FozJ(W!<*h4^&#zsG7M%I`Z zBIzL@NHV`)Z&D84A+Kxl9+`+x3Q;D+F}dgBF9HJ@I)q3h_mZEQtBf0+%cw}V3QK>& z;`kI}j`|?0ENL63K%@u4Es)oKs^r0!(nIAI@C198tKyT zEIbEu18GQ!*|69;bZyz#rpggmy(~8zVtg_(zwv3?gfAiCB9FrNk;rWnw);uW{`ph8 zu>B0Gwq&fHI0X;gk@_ByMEHWbU}$M*NEYWw(jZk}7ZSNyp3!yd*5ztjg{5jnsR}iI z`0&-{DV*e@t3Jaf^YBe7HgTMiZLWU^NXf~XWk7|M99~6D$h|rsnX|iJ*%qjqQO`~y)&66H!kzCNH$34FNHs4yT#{CxqL7BPX!S^~q@1BJjhq&2~sVQwCgNms;u5-P)XK#XAjMB@P zUeX<|i~&~Rq|N)JrI%dT#nNKKCCuYJv@O)HU6u1)RRZrUzoflixQ?@J} z@1fo>Um1O2*Q%-}4{JjUZr0DR{~5>vv6kZ>i zKt-Rmg?d{8%+!U_sqtYfxphv%PgY-w7Yr#YZXHAHX!ghgRMs6C26sJ_Nc`=vg%=Rm zrQD%mFzDQ%ZwGhitahQ)#8`p7dkeH_Ey7Ah96r8N$dSU3<;?Y)H*O%?t5q<$PB>el zS?^NEInPQ@Y6(+^i{IGTSe-(?RL^-c0yZ=-FhKBzLL6X~i4%6wgrI!gzI}~U2->=Y zmdN(;_fKjXn{$SOS)_o(K-wN0o?65_#9X}#hKFHXN+Ce7VI$mdt?hx-VSJ1<_|W$} zD#mB8&WxPnqqDaAG54`t`?K3 z3EEo8hPLKr2`MSv#52h;D>%&7jj8TJdM13HP^bQ}N0G;p890@%Rj5qXJ1#&EaQK{? zq;fSq2ggs_@wzGvp`zkqL*uu%Ht{@<%q3ZDgVbS_K^i}5t}z?TbwGPDX+|;(snaHEBBh z-n~W~=uj~bGc;(f%!bspFT+`z-ycqzK_0p!Qi^^=0jSDA3y>1Ys8#hNq&qmNKwV6Z z_?m7@+P4un*IoT*vI2vHzcW(6!6`?sPTF-StPP4n>>|~^?P0o8UR23535%osRm$UX zc$uWK<=t@&(zZ>43*AT=tZnWhaQ3M(zU%4P@#5-rJ2f@64~{QTpWAYi1;=_6;1(<2 z)9M^7m;!U$?`vBLcphkkuLQ50e*8WpgHib847*li&h^bnVYdPTG$IEGO4-5)Pe+o< zHpUlUuvl}O-g1Q(X%di-ti7~_4}EzT7Z>wGXl@AEwPWwmjN2&qx|Zh3C^mbh99+qMo-`?nZ-|Hu*BwQyVsF+(qq=N zYG!G9L3%uhhE%iy@U*tJc9x#q$r6+O;2L#op(h30x=S{RhHARm--Fb0K#4>; z7hseS*p9Nw`j9<5n2++cK8cXEZ}WeqL(&mMEod}usiArt)gXAk(Fa_fUk=471l3ep zLBV7`m6Z6PEzKVp4F0FuczA%fMm z%_J3WIKrMH{(pZUjTIKktfZpOY=8vDSk|-FFPjf*WTeVcnxNRXAv0iF%vsyGA9)TG z{p3hIWRo(_pJ%SGJzI;BerpAdli&pE`ThC6$=G-=^(RJ!&&$Kg6O)sa5t;2M-3^tF zyz=t$@!lYiUZ*sc4_$lBur{u}+2Ds0cL}^uxyKGBG;HMdt=YRb= zEJj1JL~hLQW?KuHz?e`FBC@l^#u@&bHwMJr5V6K_f!$bFU;m}@YWXS+6=}MitafS=TcED+Dc@$MP0ZAOXJ=+-&jTZ5@^kn} z*sa^Q3pO&7601=)Gw_DBK~X?3%UlkCv2~I(QqUx*bx51`_K*(aQi$OK)7TYym0DD~ z#>e)t5e#DP!-Ys%I-jU9)r1$|98FuM&5BHXhK0jk0M+!1rI|_ zjOL65yYPDg9(jx<92uyaGV`{5)gYC#Zr!?-i^~*R5HzAs&rh4KZoS}3m-Z7e1njt( z1WPgM>mW`@N=P*K^;I5!xaMn64IyMbFG};IFwpP_!|ihEr|=`&Q$7lRl$5*pputA= zHvp&6XtYG@^dVqGR2R*u!;dO_<{UakNpDnJzI_`WPQ50!<)XZZ|0Cu^7uC+K_Z)zE zfiFdcW@Hj4ga!vrKfjOZ5@m(Tb}hTK>N(9H+4wPjdPDsqre#j8-DJ;{7(odSup+up4CJkvC(_4rN2&gG@e=YVRB?WnN>DhK_%nLuwk3@YHE`HkN-5HXKfMkNA z-a_{vlC$x?yP-EeCmo>=I!X|_m0H%+Y7xiQw~CV3F2OdKk*Ny_4doar;})zBYlo>61>g^ZMhag?`SNI1g45Neoud39dgG2ZDDVjQ%HR5l-R#VbbJu zs0HK5prh?OvN(-iRhHhvx{ff|MIuecLTMd88DRT&J87Q-st0l1WQG)XbeGyTL6H4p zz7V=v2=ri>G-;pUqdi)8?%yAA&bbS3xp0SO%4Mx^#4wyT&C9i+ovg-j7&0Gjed(n&Zy|_)viXn41v=6Bw)E(A6 zz6*%_vz;fHC@f>yv=xLz^@A(Y}@Q7-Qqz z&Byl*zVMM^w!5eoBB3zc^+mH$*dPNSXmyRs*t{}X4}3r~KM-g{jdN9oeGeWyxc~rWDvQJ0)(-Bl1_#sldc_a1JcHSMN?qr#>xW=e-B)t)_2>G-7cNZuuisKKjnf-(zhr==T2RKRTmNgQR9iq1 zN=X=;K}=$aM4+s5*=gL&tao@w^w6Prl;a)R35UcA&dS)e_<0EGWF?xUE^RbHA}PSS z-dw(mtfdtvWz9Iwmgw0g_7xpz02v!}-Y+vIO^}agOW&Un=`m}d=E<|BvL){$v6cyp$!{0G_IMW+D5saUsRNh z*#SiU5|5|&))lsLUeg1R*A6&!yW`hB^;G>=-yt?}R8`V_1=|wKk$v6UA)k+iK-Nvt z5B0!yK8vGc9ogk}t$aP4^|8mN-dLt2)(J$L?w|@jzPW=cnBDV>f6bbz7>z{@(#)sH z-!Di3N1mWTQlY^;BIy*i`BIJbSEV;2+()66` z&RKo}^Ft)u8l5GkYIk(OQFyCB%cjeCw+#jP$HjOh+y84X+yaTeWk=bqG=fQ8C?7@X$OZdEmgvB?AO%w}lyN7Z;a~pmkv}G4+$tgkobFw{CI@)rN3=SINb} zb&aHkG_yY?vnLpOb_w8-^WY>I7p?i$@dR6vJc`5&!soNwQ~InjG*x&`6JW&X2#YA3YoDRx=>MRiLTU7l`}g;oohP@3 zDHJ(k7pT4)xSO7Czxw@_+dRPsXa3!MGn7`@zHN8MX!T!F)8b&opoeQMF4#{#KTP6W zk?SSj-d&;dIKev^b?Nf2cRUOHZ)lR}6U-A-ys$qEn!`Usps#Wa{?x^ab)LK&40k=Xj2A{$u$NtwtnrWP9Wf z<3}wBX4GpZOVBR>#F_M?q5)1&RmfUOCvNq{8g8s^{g@u&5^gB>zWgJy)`jSGr%Ni3?z5ulO7GoYv5j=7@bfK+sOa7Xcq?DOxN#VgB5f&|R^ zox+Xl`7d6Kyf_>tEocVL7w=UG-oQ5Wsb)H*T4R@now+~SGd02H8Xw$2f zqABqrOl|tyhSBX`Qnk21zW#|J6#+N5Vpb`MJ9Ls?r+!+-sb21JT9$JVub`%;{wQj- zaPg@bY1_sq1JFoS+BtTszM;O_ws`;X(Hkk7ZmHVHb~q1-q6K78U06B2w!GkW?GedN zD>VWg|MQ(2>!L1>qf>FIDy?~|GL4`oW+SGMnesh6j?=F2JpFAh^^@f{wRYIw_-paA z;RlhVM6cGNsHbJ1CRN`?kswkPwlLMJ$b5ezsrJq5*OoY*{zIJcXXy!Ku#f@P_QNC` z^Sq|}ypSvirL?_b@OI%LF_r_>h4{(gu1^p+c&`0Di^Qf|7I5v6=cmiBfggBY43`I& zc7TN8bwp@xZqC5FY6+VZAgUC#C6sEr{6T`7!etCIJv6NcZ%M}^0qNJHAq9livF_oQ zk)B(5EcxtVcJV6V2O^rZ^nn+;SIGhQ4jwvhazQgG3eGDVJ_)RDLr8ClaTHmwE}F;%TnC%WhARU}X*N+p=}*-ndO@ z^jL9^Wa7G&P~REGZK+c5*x1t`KBNr4{6PG6!ufwZC3znH{XI+gb-yng2!>V1)Z~TH@*iDf!6? znLm7gt-JFDux7)8>S>5lAs+rY7%M8V{#fR5K1!Cj@fUDuJx58=kANBl-3J((WMmt} z*3B9wfoCl{)_gnXmw#YP)hish(`|bgghfTOP%&U^AX>)04K;qLbtlFk&@NhC`dNjQ zXU(mX10Wi^i*1Ifg#?-I>-LkHnru(53Sdkrvr+TG!aT%(yl&BtN_3=X=Olf8+x0!; zP|;b;+AgJfy`u!vmUdwGw<3V z5ODMCPjcYsCl0g+Od6;%a${F!>LvVH=6)X5>w1TCKOpuJY4u``u|e*u9|8MK^Goa` z>C3Ma39s}&|u1e-9~Q6uDg?qj?T`{oAwM(O+9bc*Tv{X z-A%xJNBiz7xB;vAwDpB-cyx5e(r-Jh@x2egaM)HW4s_OpUdlgM4I}WCx}JL@$){9{Q{L`dT#aXxv_8Gy+gJUn=|^=Ge7+Gm&xCn zD{6X>+5mpjPUwn|Fc+V^igIIdLZJ#BN@dJkpzPu0w@Klhs(|R(;{mc2a>Sp?!Nn;wfb%JxQ%ehR({R^Owc)$1uI;dNNe`T zAQp`83=AxGQdT&`W9Fw9HJYB@ zm^xc|;VbpH-B_EmM!>CGm`M!t_@jDkWvRwX{%UF)fk@uKH9R?~C71&8Zg;jO`N!Gq z47|`ghYv0w#Q)wBFM-hk+sa+;yVpn+sm3P{f7LB@we|?P{ad~eBFoD(mCLs1&$hi@ z$MCR+bnJ-HHK`T4l4eoJ9hjvH?%lhmaGId9m;agD ze;tz=`fXC4T4s@dhU{hRXNUo3m!10#~ha>Lq(`Wvbl|j~50e zsU5=%+SbCAWn#d(;%a5Q^~1JvJ7*L>SOmjrr%L>4uhFjI@dxViV5Yt~@eJdLQ=jt& zZ_qLcrl5KJ@mwmVGqu~bp^W96+FGGGidUm`GQYK$V||L@r;@loME7I&wySwa5JbGJ zXG6+JF!=p-$n+bw>tjp@l|EF9q1!yy243=$N0HhUa>?pvARCYVli8Zu)cN!Z65t4z zWn`l@++Aop*&j4QCP+)By%lV_s1PoF{dpHxDyj`sL4OJA>GDE>xb!{bblhT6Y5lAC zH{P*-a1>*E^Z?+6J&(}XSr&&FfN6+7IBQY;<77Z}61(hKGS3<-cTu&Rvm67|lh3vU zKZuNs#L$y@6V#i9QIXRgjnWWh4sPGMCBfy&NYsM9ovSGnqxi+2Fo2e=HS4c&VC`uJrbW9dC@yz zQ}-w)02~}LqN7UdCvlA*+3BX8N&~mvomy?U$$PS#74sw5oiWn38+nrv*kv6);(kC< z*mI1K=`O-aIc)_TBvR-eOIe&h9na0=Ntf^4$x>JQF_BS;(dbL@jX1W*! zB|qJUD5qQX3(sKN`+83Igj0wz7=7nFU}$G-{CGwJ?%7Cv!sm{Pw4k@D`W{K1J8ZIW zghCd{dRCj;Q|3_*HCo}WjOoblHSMj#=mmICXszuwCqlGAJTs_b%673G+Y_}*{2 zTrW}eea+U(>n{@K@B9#H>v3l!m4B-0`|ZydUT(FH*CW$!+#$F{rPmb)M$Ga!hTS3{J93{KeW+&@%iZC?nZVWCCd21B*uSG#zD0SP0W})U zM2Ke1h@85@S3{69yu)5?sITAR(}XDo?1e(M zy7u4^bo209@EVzTt&sGiSG!f_3`XCVZ$l>2|58b*mtIW%UpUjf1yl2Z#SaQpnHGFZp%)d%C!7X=$ncz`bOsZpa_ZOI-yz z2T>1~@<4bvA{5T^;lNt|LDl!jBN=#MPc$N_>)!4AaYm>9O2T{qB91|Z{zJ}@ zj@H&wi_-H$LVBqfPbY-Wn91lev@dX|vNxR_Y9AjP%lPh#kaY+t;Jx^qZ${eMA$v}o z^tEzjlypzzNOrkX{}b2AzQ4CFDM+GplJ32X)c3s9CKh#fR zZnCHK(X**O&Z_Y(o@>eK{2?O0#y*Ms6Fme|Z~}HnB8_|+!lO(idb6ST+^}}>n}HX8 zOIP9C8wS91u>>N8^lWgifq*RM=M$S-f*TKikAeb_*)OpVu!#5gdwF{U!a7{FlW-)e z*K~z&!oKmxA!JJkq_eD=o(kd}%a~aF3YlY??0)8EQ zzY!!qgDsAl4cI^kvD9%rGCRg1qMvnovJXT0Nn{a=kPu*gf>T(NL4!=@^Dlupd4oMhlMjBx;SjdkK1Of! zt+e1R{vW?*XYU3Do#iG=_sDDA17-92(rIp Date: Mon, 7 Jul 2025 16:29:22 -0700 Subject: [PATCH 22/25] default tensors --- docs/examples/models/isotropic_thin_3d.py | 4 ++-- waveorder/models/isotropic_thin_3d.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/examples/models/isotropic_thin_3d.py b/docs/examples/models/isotropic_thin_3d.py index 3a560933..8d9f690f 100644 --- a/docs/examples/models/isotropic_thin_3d.py +++ b/docs/examples/models/isotropic_thin_3d.py @@ -29,8 +29,8 @@ ) transfer_function_arguments = { "z_position_list": (torch.arange(z_shape) - z_shape // 2) * z_pixel_size, - "numerical_aperture_illumination": 0.9, - "numerical_aperture_detection": 1.2, + "numerical_aperture_illumination": torch.tensor([0.9]), + "numerical_aperture_detection": torch.tensor([1.2]), } # Create a disk phantom diff --git a/waveorder/models/isotropic_thin_3d.py b/waveorder/models/isotropic_thin_3d.py index f8ff96c0..e47832d1 100644 --- a/waveorder/models/isotropic_thin_3d.py +++ b/waveorder/models/isotropic_thin_3d.py @@ -43,8 +43,8 @@ def calculate_transfer_function( numerical_aperture_illumination: float, numerical_aperture_detection: float, invert_phase_contrast: bool = False, - tilt_angle_zenith: float = 0.0, - tilt_angle_azimuth: float = 0.0, + tilt_angle_zenith: torch.Tensor = torch.tensor([0.0]), + tilt_angle_azimuth: torch.Tensor = torch.tensor([0.0]), ) -> Tuple[Tensor, Tensor]: transverse_nyquist = sampling.transverse_nyquist( wavelength_illumination, From 1dd0353977d95c35775ca96244f412f8abb4282c Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Fri, 1 Aug 2025 16:25:14 -0700 Subject: [PATCH 23/25] make 3d phase reconstruction compatible with tilted illumination --- waveorder/models/phase_thick_3d.py | 29 +++++++++++++++++++---------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/waveorder/models/phase_thick_3d.py b/waveorder/models/phase_thick_3d.py index 5a2e2547..07521c84 100644 --- a/waveorder/models/phase_thick_3d.py +++ b/waveorder/models/phase_thick_3d.py @@ -43,6 +43,8 @@ def calculate_transfer_function( numerical_aperture_illumination: float, numerical_aperture_detection: float, invert_phase_contrast: bool = False, + tilt_angle_zenith: torch.Tensor = torch.tensor([0.0]), + tilt_angle_azimuth: torch.Tensor = torch.tensor([0.0]), ) -> tuple[np.ndarray, np.ndarray]: transverse_nyquist = sampling.transverse_nyquist( wavelength_illumination, @@ -50,13 +52,13 @@ def calculate_transfer_function( numerical_aperture_detection, ) axial_nyquist = sampling.axial_nyquist( - wavelength_illumination, - numerical_aperture_detection, - index_of_refraction_media, + torch.tensor(wavelength_illumination), + torch.tensor(numerical_aperture_detection), + torch.tensor(index_of_refraction_media), ) yx_factor = int(np.ceil(yx_pixel_size / transverse_nyquist)) - z_factor = int(np.ceil(z_pixel_size / axial_nyquist)) + z_factor = int(np.ceil(z_pixel_size / axial_nyquist.item())) ( real_potential_transfer_function, @@ -75,6 +77,8 @@ def calculate_transfer_function( numerical_aperture_illumination, numerical_aperture_detection, invert_phase_contrast=invert_phase_contrast, + tilt_angle_zenith=tilt_angle_zenith, + tilt_angle_azimuth=tilt_angle_azimuth, ) zyx_out_shape = (zyx_shape[0] + 2 * z_padding,) + zyx_shape[1:] @@ -98,10 +102,11 @@ def _calculate_wrap_unsafe_transfer_function( numerical_aperture_illumination: float, numerical_aperture_detection: float, invert_phase_contrast: bool = False, + tilt_angle_zenith: torch.Tensor = torch.tensor([0.0]), + tilt_angle_azimuth: torch.Tensor = torch.tensor([0.0]), ) -> tuple[np.ndarray, np.ndarray]: - radial_frequencies = util.generate_radial_frequencies( - zyx_shape[1:], yx_pixel_size - ) + fyy, fxx = util.generate_frequencies(zyx_shape[1:], yx_pixel_size) + radial_frequencies = torch.sqrt(fyy**2 + fxx**2) z_total = zyx_shape[0] + 2 * z_padding z_position_list = torch.fft.ifftshift( (torch.arange(z_total) - z_total // 2) * z_pixel_size @@ -109,10 +114,14 @@ def _calculate_wrap_unsafe_transfer_function( if invert_phase_contrast: z_position_list = torch.flip(z_position_list, dims=(0,)) - ill_pupil = optics.generate_pupil( - radial_frequencies, - numerical_aperture_illumination, + ill_pupil = optics.generate_tilted_pupil( + fxx, + fyy, + torch.tensor(numerical_aperture_illumination), wavelength_illumination, + index_of_refraction_media, + tilt_angle_zenith, + tilt_angle_azimuth, ) det_pupil = optics.generate_pupil( radial_frequencies, From da77b2b0c71496eb83adb6f25fcbb54603a4c709 Mon Sep 17 00:00:00 2001 From: Talon Chandler Date: Mon, 4 Aug 2025 14:17:43 -0700 Subject: [PATCH 24/25] phase compatibility fix --- waveorder/models/phase_thick_3d.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/waveorder/models/phase_thick_3d.py b/waveorder/models/phase_thick_3d.py index 07521c84..1fd5d1f7 100644 --- a/waveorder/models/phase_thick_3d.py +++ b/waveorder/models/phase_thick_3d.py @@ -53,12 +53,13 @@ def calculate_transfer_function( ) axial_nyquist = sampling.axial_nyquist( torch.tensor(wavelength_illumination), - torch.tensor(numerical_aperture_detection), + numerical_aperture_detection, torch.tensor(index_of_refraction_media), ) - - yx_factor = int(np.ceil(yx_pixel_size / transverse_nyquist)) - z_factor = int(np.ceil(z_pixel_size / axial_nyquist.item())) + yx_factor = int( + torch.ceil(torch.tensor(yx_pixel_size / transverse_nyquist)) + ) + z_factor = int(torch.ceil(torch.tensor(z_pixel_size / axial_nyquist))) ( real_potential_transfer_function, @@ -153,7 +154,6 @@ def _calculate_wrap_unsafe_transfer_function( greens_function_z, z_pixel_size, ) - return real_potential_transfer_function, imag_potential_transfer_function From 67180f6cfeea7f3876d2b7c766edf05e353187af Mon Sep 17 00:00:00 2001 From: Sricharan Reddy Varra Date: Tue, 21 Oct 2025 14:35:36 -0700 Subject: [PATCH 25/25] added runs/ to gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 6906a010..09400a80 100644 --- a/.gitignore +++ b/.gitignore @@ -155,3 +155,4 @@ waveorder/_version.py # example data /examples/data_temp/* /logs/* +runs/* \ No newline at end of file