From b2a826b1972a9caa3393a4afd40bd3d596dcfc78 Mon Sep 17 00:00:00 2001 From: Maddie Cain Date: Wed, 27 Mar 2019 01:34:40 -0400 Subject: [PATCH] add projectors, propagators, and initializers --- CDTools/tools/__init__.py | 6 + CDTools/tools/cmath.py | 30 +- CDTools/tools/image_processing.py | 3 + CDTools/tools/initializers.py | 29 ++ CDTools/tools/projectors.py | 63 ++++ CDTools/tools/propagators.py | 161 +++++++++ build/lib/CDTools/__init__.py | 3 + build/lib/CDTools/tools/__init__.py | 4 + build/lib/CDTools/tools/cmath.py | 243 ++++++++++++++ build/lib/CDTools/tools/data.py | 344 ++++++++++++++++++++ build/lib/CDTools/tools/image_processing.py | 3 + build/lib/CDTools/tools/initializers.py | 27 ++ build/lib/CDTools/tools/losses.py | 123 +++++++ build/lib/CDTools/tools/projectors.py | 52 +++ build/lib/CDTools/tools/propagators.py | 168 ++++++++++ dist/CDTools-0.0.1-py2.7.egg | Bin 0 -> 26496 bytes tests/tools/test_initializers.py | 24 ++ tests/tools/test_projectors.py | 21 ++ tests/tools/test_propagators.py | 48 +++ 19 files changed, 1336 insertions(+), 16 deletions(-) create mode 100644 CDTools/tools/image_processing.py create mode 100644 CDTools/tools/initializers.py create mode 100644 CDTools/tools/projectors.py create mode 100644 CDTools/tools/propagators.py create mode 100644 build/lib/CDTools/__init__.py create mode 100644 build/lib/CDTools/tools/__init__.py create mode 100644 build/lib/CDTools/tools/cmath.py create mode 100644 build/lib/CDTools/tools/data.py create mode 100644 build/lib/CDTools/tools/image_processing.py create mode 100644 build/lib/CDTools/tools/initializers.py create mode 100644 build/lib/CDTools/tools/losses.py create mode 100644 build/lib/CDTools/tools/projectors.py create mode 100644 build/lib/CDTools/tools/propagators.py create mode 100644 dist/CDTools-0.0.1-py2.7.egg create mode 100644 tests/tools/test_initializers.py create mode 100644 tests/tools/test_projectors.py create mode 100644 tests/tools/test_propagators.py diff --git a/CDTools/tools/__init__.py b/CDTools/tools/__init__.py index 47a8cd2..6ea3699 100644 --- a/CDTools/tools/__init__.py +++ b/CDTools/tools/__init__.py @@ -2,3 +2,9 @@ from __future__ import division, print_function, absolute_import from CDTools.tools import cmath from CDTools.tools import losses +from CDTools.tools import data +from CDTools.tools import image_processing +from CDTools.tools import initializers +from CDTools.tools import losses +from CDTools.tools import projectors +from CDTools.tools import propagators diff --git a/CDTools/tools/cmath.py b/CDTools/tools/cmath.py index 34ab6a4..2335033 100644 --- a/CDTools/tools/cmath.py +++ b/CDTools/tools/cmath.py @@ -98,7 +98,7 @@ def cabs(a): def cphase(a): - """Returns the complex conjugate of a complex torch tensor + """Returns the phase of a complex torch tensor Pytorch uses tensors with a final dimension of 2 to represent complex numbers. This calculates the elementwise complex phase @@ -112,8 +112,8 @@ def cphase(a): """ return t.atan2(a[...,1],a[...,0]) - - + + def cconj(a): """Returns the complex conjugate of a complex torch tensor @@ -147,7 +147,7 @@ def cmult(a,b): torch.Tensor : A tensor storing the elementwise product """ - + real = a[...,0] * b[...,0] - a[...,1] * b[...,1] imag = a[...,0] * b[...,1] + a[...,1] * b[...,0] return t.stack((real,imag),dim=-1) @@ -171,7 +171,7 @@ def cdiv(a,b): return cmult(a, cconj(b)) / t.unsqueeze(cabssq(b),-1) - + # # Not entirely sure if these belong here, but heck with it. # We just need the ability to do fftshifts @@ -180,23 +180,23 @@ def cdiv(a,b): def fftshift(array,dims=None): """Drop-in torch replacement for scipy.fftpack.fftshift - + This maps a tensor, assumed to be the output of a fast Fourier transform, into a tensor whose zero-frequency element is at the center of the tensor instead of the start. It will by default shift every dimension in the tensor but the last (which is assumed to represent the complex number and be of dimension 2), but can shift any arbitrary set of dimensions. - + Args: array (torch.Tensor) : An array of data to be fftshifted dims (iterable) : A list of all dimensions to shift Returns: torch.Tensor : fftshifted tensor - + """ - + if dims is None: dims=list(range(array.dim()))[:-1] for dim in dims: @@ -211,7 +211,7 @@ def fftshift(array,dims=None): def ifftshift(array,dims=None): """Drop-in torch replacement for scipy.fftpack.iftshift - + This maps a tensor, assumed to be the shifted output of a fast Fourier transform, into a tensor whose zero-frequency element is back at the start of the tensor instead of the center. It is the @@ -219,25 +219,23 @@ def ifftshift(array,dims=None): every dimension in the tensor but the last (which is assumed to represent the complex number and be of dimension 2), but can shift any arbitrary set of dimensions. - + Args: array (torch.Tensor) : An array of data to be ifftshifted dims (iterable) : A list of all dimensions to shift Returns: torch.Tensor : ifftshifted tensor - + """ - + if dims is None: dims=list(range(array.dim()))[:-1] for dim in dims: length = array.size()[dim] cut_to = length // 2 cut_len = length - cut_to - + array = t.cat((array.narrow(dim,cut_to,cut_len), array.narrow(dim,0,cut_to)), dim) return array - - diff --git a/CDTools/tools/image_processing.py b/CDTools/tools/image_processing.py new file mode 100644 index 0000000..22c4a17 --- /dev/null +++ b/CDTools/tools/image_processing.py @@ -0,0 +1,3 @@ +from __future__ import division, print_function, absolute_import +import numpy as np +import torch as t diff --git a/CDTools/tools/initializers.py b/CDTools/tools/initializers.py new file mode 100644 index 0000000..dbe6361 --- /dev/null +++ b/CDTools/tools/initializers.py @@ -0,0 +1,29 @@ +from __future__ import division, print_function, absolute_import +import numpy as np +import torch as t + +all = ['gaussian'] + + +def gaussian(shape, amplitude, sigma, center = None): + """Returns an array with a centered gaussian + + Takes in the shape, amplitude, and standard deviation of a gaussian + and returns an array with values corresponding to a two-dimensional gaussian function + z = amplitude*exp(-(x-center[0])**2/sigma[0]**2+(y-center[1])**2/sigma[1]**2) + Note that [0, 0] is taken to be at the upper left corner of the array. + Default is centered at ((shape[0]-1)/2, (shape[1]-1)/2)) because x and y are zero-indexed. + + Args: + shape (array_like) : A 1x2 array-like object specifying the dimensions of the output array in the form (y shape, x shape) + amplitude (float or int): The amplitude the gaussian to simulate + sigma (array_like): A 1x2 array-like object specifying the x- and y- standard deviation of the gaussian in the form (y stdev, y stdev) + center (array_like) : Optional 1x2 array-like object specifying the location of the center of the gaussian (y center, x center) + + Returns: + numpy.array : The real-valued gaussian array + """ + if center is None: + center = ((shape[0]-1)/2, (shape[1]-1)/2) + y, x = np.mgrid[:shape[0], :shape[1]] + return amplitude*np.exp(-((x-center[1])/sigma[1])**2-((y-center[0])/sigma[0])**2) diff --git a/CDTools/tools/projectors.py b/CDTools/tools/projectors.py new file mode 100644 index 0000000..b2aa8cf --- /dev/null +++ b/CDTools/tools/projectors.py @@ -0,0 +1,63 @@ +from __future__ import division, print_function, absolute_import +from CDTools.tools.cmath import * +import torch as t + +__all__ = ['modulus', 'support'] + +def modulus(wavefront, intensities, mask = None): + """Implements the modulus constraint in torch + + This accepts a torch tensor representing the propagated simulated wavefront(s), + where the last dimension represents the real and imaginary components of + the propagated wavefield(s). It projects the modulus of the diffraction pattern + onto the modulus of the simulated wavefield. + + It assumes that the wavefront is stored in an array + [i,j] where i corresponds to the y-axis and j corresponds to the + x-axis, with the origin following the CS standard of being in the + upper right. + + Args: + wavefront (torch.Tensor) : The JxNxMx2 stack of complex propagated wavefronts + intensities (torch.Tensor): The measured diffraction pattern(s) stored as an JxNxM stack of real tensors + mask (torch.Tensor) : Mask for the intensities array with shape JxNxM, where bad detector pixels are set to 0 and usable pixels set to 1 + Returns: + torch.Tensor : The JxNxMx2 propagated wavefield with corrected intensities + """ + # Calculate amplitudes from intensities + amplitudes = intensities**.5 + # Normalize wavefront so the complex elements have modulus one + abs = cabs(wavefront) + wavefront[...,0]/=abs + wavefront[...,1]/=abs + if mask is None: + # Replace amplitude of wavefront with measured amplitude + wavefront[...,0]*=amplitudes + wavefront[...,1]*=amplitudes + return wavefront + else: + return wavefront[mask != 0] + + +def support(wavefront, support): + """Implements the support constraint in torch + + This accepts a torch tensor representing the propagated simulated wavefront(s), + where the last dimension represents the real and imaginary components of + the propagated wavefield(s). It projects the support of the imaged object + onto the simulated wavefront via a mask. + + It assumes that the wavefront is stored in an array + [i,j] where i corresponds to the y-axis and j corresponds to the + x-axis, with the origin following the CS standard of being in the + upper right. + + Args: + wavefront (torch.Tensor) : The JxNxMx2 stack of complex propagated wavefronts + mask (torch.Tensor) : Mask for the intensities array with shape JxNxM, where bad detector pixels are set to 0 and usable pixels set to 1 + Returns: + torch.Tensor : The JxNxMx2 wavefield with the mask applied + """ + wavefront[...,0] *= support + wavefront[...,1] *= support + return wavefront diff --git a/CDTools/tools/propagators.py b/CDTools/tools/propagators.py new file mode 100644 index 0000000..751f8d9 --- /dev/null +++ b/CDTools/tools/propagators.py @@ -0,0 +1,161 @@ +from __future__ import division, print_function, absolute_import +from CDTools.tools.cmath import * +import torch as t + +__all__ = ['far_field', 'near_field', 'inverse_far_field', 'inverse_near_field', 'get_exit_waves'] + + +def far_field(wavefront): + """Implements a far-field propagator in torch + + This accepts a torch tensor, where the last dimension + represents the real and imaginary components of the wavefield, + and returns the far-field propagated version of it assuming it matches the + detector dimensions. It assumes that the + propagation is purely far-field, without checking that the geometry + is consistent with that assumption. + + + It also assumes that the real space wavefront is stored in an array + [i,j] where i corresponds to the y-axis and j corresponds to the + x-axis, with the origin following the CS standard of being in the + upper right. The zero frequency component of the propagated wavefield is + shifted to the center of the array. + + Args: + wavefront (torch.Tensor) : The JxNxMx2 stack of complex wavefronts to be propagated + Returns: + torch.Tensor : The JxNxMx2 propagated wavefield + """ + + return fftshift(t.fft(wavefront, 2)) + +def inverse_far_field(wavefront): + """Implements the inverse of the far-field propagator in torch + + This accepts a torch tensor, where the last dimension + represents the real and imaginary components of the propagated wavefield, + and returns the un-propagated array. + + It assumes that the real space wavefront is stored in an array + [i,j] where i corresponds to the y-axis and j corresponds to the + x-axis, with the origin following the CS standard of being in the + upper right. The zero frequency component of the propagated wavefield is + assumed to be the center of the array. + + Args: + wavefront (torch.Tensor) : The JxNxMx2 stack of complex wavefronts propagated to the far-field + Returns: + torch.Tensor : The JxNxMx2 exit wavefield + """ + return t.ifft(ifftshift(wavefront), 2) + + +def generate_angular_spectrum_propagator(shape, spacing, wavelength, z): + """Generates an angular-spectrum based near-field propagator from experimental quantities + + This function generates an angular-spectrum based near field + propagator that will work on torch Tensors. The function is structured + this way - to generate the propagator first - because the + generation of the propagation mask is a bit expensive and if this + propagator is used in a reconstruction program, then it will be best + to calculate this mask once and close over it. + + Args: + shape (iterable) : The shape of the arrays to be propagated + spacing (iterable) : The pixel size in each dimension of the arrays to be propagated + wavelength (float) : The wavelength of light to simulate propagation of + z (float) : The distance to simulate propagation over + Returns: + torch.Tensor : A propagation term which accounts for the phase change that each plane wave will undergo on its journey to the prediction plane. + """ + + ki = fftpack.fftfreq(shape[0],spacing[0]) + kj = fftpack.fftfreq(shape[1],spacing[1]) + Ki, Kj = np.meshgrid(ki,kj) + propagator = np.exp(1j*np.sqrt((2*np.pi/wavelength)**2 + - Ki**2 - Kj**2) * z) + propagator = complex_to_float(propagator).astype(np.float32) + propagator = t.from_numpy(propagator).cuda() + + return propagator + + +def near_field(wavefront, angular_spectrum_propagator): + """This function accepts an 3d torch tensor, where the + last dimension represents the real and imaginary components of + the wavefield, and returns the near-field propagated version of it. + + + Args: + angular_spectrum_propagator (torch.Tensor) : The near field propagator + wavefront (torch.Tensor) : The JxNxMx2 stack of complex wavefronts to be propagated + Returns: + function : The wavefront propagated to the near field + """ + + return t.ifft(angular_spectrum_propagator * t.fft(wavefront,2), 2) + + + +def inverse_near_field(wavefront, angular_spectrum_propagator): + """This function accepts a 3d torch tensor, where the + last dimension represents the real and imaginary components of + the near-field propagated wavefield, and returns the exit wavefront via an inverse transformation. + + + Args: + angular_spectrum_propagator (torch.Tensor) : The pixel size in each dimension of the arrays to be propagated + wavefront (torch.Tensor) : The JxNxMx2 stack of complex wavefronts to be propagated + Returns: + function : A function to propagate a torch tensor. + """ + return t.ifft(t.fft(wavefront,2) * angular_spectrum_propagator**-1, 2) + + + + +def get_exit_waves(probe, object, translations): + """Returns a stack of exit waves accounting for subpixel shifts + + This function returns a collection of exit waves, with the first + dimension as the translation index and the final dimensions + corresponding to the detector. The exit waves are calculated by + shifting the object with each translation in turn, using linear + interpolation. + Args: + probe (torch.Tensor) : An MxM probe function for the exit waves + object (torch.Tensor) : The object function to be probed + translations (torch.Tensor) : The Nx2 array of translations to simulate + Returns: + torch.Tensor : An NxMxM tensor of the calculated exit waves + """ + + # Separate the translations into a part that chooses the window + # And a part that defines the windowing function + integer_translations = t.floor(translations) + subpixel_translations = translations - integer_translations + integer_translations = integer_translations.to(dtype=t.int32) + + selections = [] + for tr, sp in zip(integer_translations, + subpixel_translations): + + sel00 = object[tr[0]:tr[0]+probe.shape[0], + tr[1]:tr[1]+probe.shape[1]] + + sel01 = object[tr[0]:tr[0]+probe.shape[0], + tr[1]+1:tr[1]+1+probe.shape[1]] + + sel10 = object[tr[0]+1:tr[0]+1+probe.shape[0], + tr[1]:tr[1]+probe.shape[1]] + + sel11 = object[tr[0]+1:tr[0]+1+probe.shape[0], + tr[1]+1:tr[1]+1+probe.shape[1]] + + selections.append(sel00 * (1-sp[0])*(1-sp[1]) + \ + sel01 * (1-sp[0])*sp[1] + \ + sel10 * sp[0]*(1-sp[1]) + \ + sel11 * sp[0]*sp[1]) + + return t.stack([cmult(probe,selection) for selection in selections]) diff --git a/build/lib/CDTools/__init__.py b/build/lib/CDTools/__init__.py new file mode 100644 index 0000000..fbe00ae --- /dev/null +++ b/build/lib/CDTools/__init__.py @@ -0,0 +1,3 @@ +from __future__ import division, print_function, absolute_import + +from CDTools import tools diff --git a/build/lib/CDTools/tools/__init__.py b/build/lib/CDTools/tools/__init__.py new file mode 100644 index 0000000..47a8cd2 --- /dev/null +++ b/build/lib/CDTools/tools/__init__.py @@ -0,0 +1,4 @@ +from __future__ import division, print_function, absolute_import + +from CDTools.tools import cmath +from CDTools.tools import losses diff --git a/build/lib/CDTools/tools/cmath.py b/build/lib/CDTools/tools/cmath.py new file mode 100644 index 0000000..34ab6a4 --- /dev/null +++ b/build/lib/CDTools/tools/cmath.py @@ -0,0 +1,243 @@ +"""Contains basic functions for dealing with complex numbers in pytorch. + +Since pytorch doesn't have built-in support for complex numbers, but the +fast fourier transforms in pytorch assume a specific format for complex +arrays, this module uses that format to store complex numbers. It exposes +functions for converting between complex numpy arrays and torch tensors +stored in that format, as well as basic complex math operations implemented +on the torch tensors +""" +from __future__ import division, print_function, absolute_import +import numpy as np +import torch as t + + +__all__ = ['complex_to_torch','torch_to_complex','cabssq','cabs','cconj', + 'cmult', 'cdiv', 'cphase', 'fftshift', 'ifftshift'] + + +# +# These define the conversions to and from this format +# + +def complex_to_torch(x): + """Maps a complex numpy array to a torch tensor + + Pytorch uses tensors with a final dimension of 2 to represent + complex numbers. This maps a complex type numpy array to a torch + tensor following this convention + + Args: + x (array_like): A numpy array to convert + + Returns: + torch.Tensor : A torch tensor representation of that array + + """ + return t.from_numpy(np.stack((np.real(x),np.imag(x)),axis=-1)) + + +def torch_to_complex(x): + """Maps a torch tensor to the a complex numpy array + + Pytorch uses tensors with a final dimension of 2 to represent + complex numbers. This maps a torch tensor following that convention + to the appropriate numpy complex array + + Args: + x (torch.Tensor): A tensor to convert + + Returns: + np.array : A complex typed numpy array corresponding to the input + + """ + x = np.array(x) + x = x[...,0] + x[...,1] * 1j + return x + + +# +# And these define the basic operations on these arrays. Note that +# multiplication between a complex valued and real valued pytorch +# tensor will proceed as expected because of torch's broadcasting +# and thus doesn't need it's own function +# + +def cabssq(a): + """Returns the square of the absolute value of a complex torch tensor + + Pytorch uses tensors with a final dimension of 2 to represent + complex numbers. This calculates the elementwise absolute value + squared of any toch tensor following that standard. + + Args: + x (torch.Tensor): An input tensor + + Returns: + array_like : A tensor storing the elementwise absolute value squared + + """ + return a[...,0]**2 + a[...,1]**2 + + +def cabs(a): + """Returns the absolute value of a complex torch tensor + + Pytorch uses tensors with a final dimension of 2 to represent + complex numbers. This calculates the elementwise absolute value + of any torch tensor following that standard. + + Args: + x (torch.Tensor): An input tensor + + Returns: + array_like : A tensor storing the elementwise absolute value + + """ + return t.sqrt(cabssq(a)) + + +def cphase(a): + """Returns the complex conjugate of a complex torch tensor + + Pytorch uses tensors with a final dimension of 2 to represent + complex numbers. This calculates the elementwise complex phase + of any torch tensor following that standard. + + Args: + x (torch.Tensor): An input tensor + + Returns: + array_like : A tensor storing the elementwise phase + + """ + return t.atan2(a[...,1],a[...,0]) + + +def cconj(a): + """Returns the complex conjugate of a complex torch tensor + + Pytorch uses tensors with a final dimension of 2 to represent + complex numbers. This calculates the elementwise complex conjugate + of any torch tensor following that standard. + + Args: + x (torch.Tensor): An input tensor + + Returns: + array_like : A tensor storing the elementwise complex conjugate + + """ + return t.stack((a[...,0],-a[...,1]),dim=-1) + + + +def cmult(a,b): + """Returns the complex product of two torch tensors + + Pytorch uses tensors with a final dimension of 2 to represent + complex numbers. This calculates the elementwise product + of two torch tensors following that standard. + + Args: + a (torch.Tensor): An input tensor + b (torch.Tensor): A second input tensor + + Returns: + torch.Tensor : A tensor storing the elementwise product + + """ + + real = a[...,0] * b[...,0] - a[...,1] * b[...,1] + imag = a[...,0] * b[...,1] + a[...,1] * b[...,0] + return t.stack((real,imag),dim=-1) + + +def cdiv(a,b): + """Returns the complex quotient of two torch tensors + + Pytorch uses tensors with a final dimension of 2 to represent + complex numbers. This calculates the elementwise quotient + of two torch tensors following that standard. + + Args: + a (torch.Tensor): An input tensor + b (torch.Tensor): A second input tensor + + Returns: + torch.Tensor : A tensor storing the elementwise complex quotient + + """ + return cmult(a, cconj(b)) / t.unsqueeze(cabssq(b),-1) + + + +# +# Not entirely sure if these belong here, but heck with it. +# We just need the ability to do fftshifts +# + + +def fftshift(array,dims=None): + """Drop-in torch replacement for scipy.fftpack.fftshift + + This maps a tensor, assumed to be the output of a fast Fourier + transform, into a tensor whose zero-frequency element is at the + center of the tensor instead of the start. It will by default shift + every dimension in the tensor but the last (which is assumed to + represent the complex number and be of dimension 2), but can shift + any arbitrary set of dimensions. + + Args: + array (torch.Tensor) : An array of data to be fftshifted + dims (iterable) : A list of all dimensions to shift + + Returns: + torch.Tensor : fftshifted tensor + + """ + + if dims is None: + dims=list(range(array.dim()))[:-1] + for dim in dims: + length = array.size()[dim] + cut_to = (length + 1) // 2 + cut_len = length - cut_to + array = t.cat((array.narrow(dim,cut_to,cut_len), + array.narrow(dim,0,cut_to)), dim) + return array + + + +def ifftshift(array,dims=None): + """Drop-in torch replacement for scipy.fftpack.iftshift + + This maps a tensor, assumed to be the shifted output of a fast + Fourier transform, into a tensor whose zero-frequency element is + back at the start of the tensor instead of the center. It is the + inverse of the fftshift operator. It will by default shift + every dimension in the tensor but the last (which is assumed to + represent the complex number and be of dimension 2), but can shift + any arbitrary set of dimensions. + + Args: + array (torch.Tensor) : An array of data to be ifftshifted + dims (iterable) : A list of all dimensions to shift + + Returns: + torch.Tensor : ifftshifted tensor + + """ + + if dims is None: + dims=list(range(array.dim()))[:-1] + for dim in dims: + length = array.size()[dim] + cut_to = length // 2 + cut_len = length - cut_to + + array = t.cat((array.narrow(dim,cut_to,cut_len), + array.narrow(dim,0,cut_to)), dim) + return array + + diff --git a/build/lib/CDTools/tools/data.py b/build/lib/CDTools/tools/data.py new file mode 100644 index 0000000..4ef1796 --- /dev/null +++ b/build/lib/CDTools/tools/data.py @@ -0,0 +1,344 @@ +from __future__ import division, print_function, absolute_import + +import h5py +import numpy as np + +__all__ = ['get_entry_info', + 'get_sample_info', + 'get_wavelength', + 'get_detector_geometry', + 'get_mask', + 'get_data', + 'get_ptycho_translations'] + +# +# +# I will put here some thoughts about how to load data into this program. +# +# +# The reconstructions should have the ability to generate datasets. +# So you could write a reconstruction engine and then it would be +# able to simulate data directly in the engine for you to use as a +# reconstruction +# +# I don't even think there needs to be a loading tool for loading cxi files +# because there isn't really a better method beyond just loading the +# file into an h5py object. This file could host the simple cxi file +# browser, perhaps. But I think the reality is that we need individual +# loaders for each kind of experiment. Perhaps we could put some basic +# reuseable tools for inspecting cxi-type h5 files in this file. +# +# +# Then, there can be some more sophisticated tools that load data for +# specific use cases that are common - loading data for a 2D CDI experiment, +# loading data for a 2D Ptycho experiment, loading data for Bragg Ptycho in +# 3D, loading data for a 3D CDI experiment, etc. +# +# +# Perhaps one good way to package this is for the kind of data associated +# with any particular experiment to have it's own kind of dataset or view. +# So there would be a "2D Ptychography" data viewer, which would contain +# all the measured data that comes from a 2D ptychography experiment. +# The specialized functions would plop out these data viewers, and the +# reconstruction classes could be designed around a particular kind of +# viewer with the most general kind just requiring a generic data viewer. +# +# Data viewers could have simple tools like the ability to send themselves +# to the GPU, CPU, change the datatype, etc. I think the most generic thing +# is as a subclass of the torch Data objects, where they would for each slice +# return the index, a set of defining parameters (translation, angle, energy, +# whatever), and a diffraction pattern. They would also have a "setup" +# attribute, or "metadata", or whatever you'd want to call it, that contain +# the various fixed experimental parameters (energy, distance, etc.) +# +# And I think the cxi visualizer should really go into it's own script, +# because it's not a reuseable component. +# + + +# +# Functions to inspect the basic attributes of a cxi file represented as an +# h5 file object +# + +def get_entry_info(cxi_file): + """Returns a dictionary with the basic metadata from the cxi file's entry_1 attribute + + Args: + cxi_file (h5py.File) : a file object to be read + + Returns: + dict : A dictionary with basic metadata defined in the cxi file + + """ + e1 = cxi_file['entry_1'] + metadata_attrs = ['title', + 'experiment_identifier', + 'experiment_description', + 'program_name', + 'start_time', + 'end_time'] + metadata = {attr: str(e1[attr][()].decode()) for attr in metadata_attrs + if attr in e1} + return metadata + + +def get_sample_info(cxi_file): + """Returns a dictionary with the basic metadata from the cxi file's entry_1/sample_1 attribute + + Args: + cxi_file (h5py.File) : a file object to be read + + Returns: + dict : A dictionary with basic metadata from the sample defined in the cxi file + + """ + if 'entry_1/sample_1' not in cxi_file: + return None + + s1 = cxi_file['entry_1/sample_1'] + metadata_attrs = ['name','description','unit_cell_group'] + metadata = {attr: str(s1[attr][()].decode()) for attr in metadata_attrs + if attr in s1} + + float_attrs = ['concentration', + 'mass', + 'temperature', + 'thickness', + 'unit_cell_volume'] + for attr in float_attrs: + if attr in s1: + metadata[attr] = np.float32(s1[attr][()]) + + if 'unit_cell' in s1: + metadata['unit_cell'] = np.array(s1['unit_cell']).astype(np.float32) + + # TODO: Add my nonstandard "surface normal" attribute here + + # TODO: I should add the sample geometry as a valid metadata that can + # be copied over + + return metadata + + +def get_wavelength(cxi_file): + """Returns the wavelength of the source defined in the cxi file object, in m + + Args: + cxi_file (h5py.File) : a file object to be read + + Returns: + np.float32 : The wavelength of the source defined in the cxi file + """ + i1 = cxi_file['entry_1/instrument_1'] + if 'source_1/wavelength' in i1: + wavelength = np.float32(i1['source_1/wavelength']) + elif 'source_1/energy' in i1: + energy = np.float32(i1['source_1/energy']) + wavelength = 1.9864459e-25 / energy + else: + raise KeyError('Neither Wavelength or Energy Defined in provided .cxi File') + + return wavelength + + +def get_detector_geometry(cxi_file): + """Returns a standardized description of the detector geometry defined in the cxi file object + + It makes intelligent assumptions based on the definitions in the cxi + file definition. The standardized description of the geometry that it + outputs includes the sample to detector distance, the corner location + of the detector, and the basis vectors defining the detector. It can + only handle detectors defined as rectangular grids of pixels. + + The distance and corner_location values are technically overdetermining + the detector location, but for many experiments (particularly + transmission experiments), the distance is needed and the exact + corner location is not. If the corner location is not reported in + the cxi file, no attempt will be made to calculate it. + + Args: + cxi_file (h5py.File) : a file object to be read + + Returns: + distance (np.float32) : The sample to detector distance, in m + basis_vectors (np.array) : The basis vectors for the detector + corner_location (np.array) : The location of the (0,0) pixel in the detector + + """ + i1 = cxi_file['entry_1/instrument_1'] + d1 = i1['detector_1'] + + if 'detector_1/basis_vectors' in i1: + basis_vectors = np.array(d1['basis_vectors']) + else: + # This whole thing just to account for all the ways people can + # implicitly define the x or y pixel size for a detector. I've + # seen too many of these in the wild, unfortunately... + try: + x_pixel_size = np.float32(d1['x_pixel_size']) + except: + x_pixel_size = None + try: + y_pixel_size = np.float32(d1['y_pixel_size']) + except: + y_pixel_size = None + + if x_pixel_size is None and y_pixel_size is not None: + x_pixel_size = y_pixel_size + elif x_pixel_size is not None and y_pixel_size is None: + y_pixel_size = x_pixel_size + if x_pixel_size is None and y_pixel_size is None: + raise KeyError('Detector pixel size not defined in file.') + basis_vectors = np.array([[0,-y_pixel_size,0], + [-x_pixel_size,0,0]]).transpose() + + try: + distance = np.float32(d1['distance']) + except: + distance = None + try: + corner_position = np.array(d1['corner_position']) + except: + corner_position = None + + # Don't pretend to calculate corner position from distance if it's + # if it's not defined, but do calculate distance from corner position + # if distance is not defined. If neither is defined, then raise + # an error. + if distance is None and corner_position is not None: + detector_normal = np.cross(basis_vectors[:,0], + basis_vectors[:,1]) + detector_normal /= np.linalg.norm(detector_normal) + distance = np.linalg.norm(np.dot(corner_position, detector_normal)) + + if distance is None and corner_position is not None: + raise KeyError('Neither sample to detector distance or corner position is defined in file.') + + return distance, basis_vectors, corner_position + + +def get_mask(cxi_file): + """Returns the detector mask defined in the cxi file object + + This function converts from the format specified in the cxi file + definition to a simple on/off mask, where a value of 1 defines a + good pixel (on) and a value of 0 defines a bad pixel (off). + + If any bit is set in the mask at all, it will be defined as a bad + pixel, with the exception of pixels marked exactly as 0x00001000, + which is defined to mean that the pixel has signal above the + background. These pixels are treated as on pixels + + Args: + cxi_file (h5py.File) : a file object to be read + + Returns: + np.array : An array storing the mask from the cxi file + """ + + i1 = cxi_file['entry_1/instrument_1'] + if 'detector_1/mask' in i1: + mask = np.array(i1['detector_1/mask']).astype(np.uint32) + mask_on = np.equal(mask,np.uint32(0)) + mask_has_signal = np.equal(mask,np.uint32(0x00001000)) + return np.logical_or(mask_on,mask_has_signal).astype(np.uint8) + else: + return None + + +def get_data(cxi_file): + """Returns an array with the full stack of detector data defined in the cxi file object + + This function will make sure to check all the various places that it's + okay to store the data in, to ensure that it can find the data regardless + of whether the creator of the .cxi file has remembered to link the data + to all the required locations. + + It will return the data array in whatever shape it's defined in. + + It will also read out the axes attribute of the data into a list + of strings + + Args: + cxi_file (h5py.File) : a file object to be read + + Returns: + np.array : An array storing the data defined in the cxi file + list : A list of the axes defined in the axes attribute, if any + """ + # Possible locations for the data + # + # entry_1/detector_1/data + if 'entry_1/data_1/data' in cxi_file: + pull_from = 'entry_1/data_1/data' + elif 'entry_1/instrument_1/detector_1/data' in cxi_file: + pull_from = 'entry_1/instrument_1/detector_1/data' + else: + raise KeyError('Data is not defined within cxi file') + data = np.array(cxi_file[pull_from]).astype(np.float32) + if 'axes' in cxi_file[pull_from].attrs: + axes = str(cxi_file[pull_from].attrs['axes'].decode()).split(':') + axes = [axis.strip().lower() for axis in axes] + else: + axes = None + + return data, axes + + + +def get_ptycho_translations(cxi_file): + """Gets an array of x,y,z translations, if such an array has been defined in the file + + It applies two operations to the translations. First, it negates them, + because the CXI file format is designed to specify translations of the + samples and the CDTools code specifies translations of the optics. + Second, it transposes the array so that the first axis is translation + ID and the second axis is the (x,y,z) components of the translation + + Args: + cxi_file (h5py.File) : a file object to be read + + Returns: + np.array : An array storing the translations defined in the cxi file + list : A list of the axes defined in the axes attribute, if any + """ + + if 'entry_1/data_1/translation' in cxi_file: + pull_from = 'entry_1/data_1/translation' + elif 'entry_1/sample_1/geometry_1/translation' in cxi_file: + pull_from = 'entry_1/sample_1/geometry_1/translation' + elif 'entry_1/instrument_1/detector_1/translation' in cxi_file: + pull_from = 'entry_1/instrument_1/detector_1/translation' + else: + raise KeyError('Translations are not defined within cxi file') + + translations = -np.array(cxi_file[pull_from]).astype(np.float32).transpose() + return translations + + +# +# It might be useful to make some helper functions to help write cxi files +# + +# +# A function to place the skeleton of a cxi file down +# + +# +# A function to define the source attributes +# + +# +# A function to define the detector geometry +# + +# +# A function to save out a mask, converting it to the correct format +# + +# +# Perhaps a function to store the data and link it correctly? But this might +# have to change too much situation to situation +# + diff --git a/build/lib/CDTools/tools/image_processing.py b/build/lib/CDTools/tools/image_processing.py new file mode 100644 index 0000000..22c4a17 --- /dev/null +++ b/build/lib/CDTools/tools/image_processing.py @@ -0,0 +1,3 @@ +from __future__ import division, print_function, absolute_import +import numpy as np +import torch as t diff --git a/build/lib/CDTools/tools/initializers.py b/build/lib/CDTools/tools/initializers.py new file mode 100644 index 0000000..a7184c6 --- /dev/null +++ b/build/lib/CDTools/tools/initializers.py @@ -0,0 +1,27 @@ +from __future__ import division, print_function, absolute_import +import numpy as np +import torch as t + +all = ['gaussian'] + + +def gaussian(shape, amplitude, sigma, center = None): + """Returns an array with a centered gaussian + + Takes in the shape, amplitude, and standard deviation of a gaussian + and returns an array with values corresponding to a two-dimensional gaussian function + z = amplitude*exp(-(x-center[0])**2/sigma[0]**2+(y-center[1])**2/sigma[1]**2) + Note that [0, 0] is taken to be at the upper left corner of the array. + Default is centered at (shape[0]/2, shape[1]/2). + + Args: + shape (array_like) : A 1x2 array-like object specifying the dimensions of the output array + amplitude (float or int): The amplitude the gaussian to simulate + sigma (array_like): A 1x2 array-like object specifying the x- and y- standard deviation of the gaussian + center (array_like) : Optional 1x2 array-like object specifying the location of the center of the gaussian + + Returns: + torch.Tensor : The real-valued gaussian array + """ + x, y = np.meshgrid(shape) + return x diff --git a/build/lib/CDTools/tools/losses.py b/build/lib/CDTools/tools/losses.py new file mode 100644 index 0000000..8255ac8 --- /dev/null +++ b/build/lib/CDTools/tools/losses.py @@ -0,0 +1,123 @@ +"""Contains various loss functions to be used for optimization + +It exposes three losses, one returning the mean squared amplitude error, one +that returns the mean squared intensity error, and one that returns the +maximum likelihood metric for a system with Poisson statistics. + +""" +from __future__ import division, print_function, absolute_import +import torch as t + + +__all__ = ['amplitude_mse', 'intensity_mse', 'poisson_nll'] + + +def amplitude_mse(intensities, sim_intensities, mask=None): + """ Returns the mean squared error of a simulated dataset's amplitudes + + Calculates the mean squared error between a given set of + measured diffraction intensities and a simulated set. + + This function calculates the mean squared error between their + associated amplitudes. Because this is not well defined for negative + numbers, make sure that all the intensities are >0 before using this + loss. + + It can accept intensity and simulated intensity tensors of any shape + as long as their shapes match, and the provided mask array can be + broadcast correctly along them. + + This is empirically the most useful loss function + + Args: + intensities (torch.Tensor) : A tensor with measured detector values + sim_intensities (torch.Tensor) : A tensor of simulated detector intensities + mask (torch.Tensor) : A mask with ones for pixels to include and zeros for pixels to exclude + + Returns: + loss (torch.Tensor) : A single value for the summed mse + + """ + + # I know it would be more efficient if this function took in the + # amplitudes instead of the intensities, but I want to be consistent + # with all the errors working off of the same inputs + + if mask is None: + return t.sum((t.sqrt(sim_intensities) - + t.sqrt(intensities))**2) / intensities.view(-1).shape[0] + else: + masked_intensities = intensities.masked_select(mask) + return t.sum((t.sqrt(sim_intensities.masked_select(mask)) - + t.sqrt(masked_intensities))**2) / masked_intensities.shape[0] + + + +def intensity_mse(intensities, sim_intensities, mask=None): + """ Returns the mean squared error of a simulated dataset's intensities + + Calculates the summed mean squared error between a given set of + diffraction intensities - the measured set of detector intensities - + and a simulated set of diffraction intensities. This function + calculates the mean squared error between the intensities. + + It can accept intensity and simulated intensity tensors of any shape + as long as their shapes match, and the provided mask array can be + broadcast correctly along them. + + Args: + intensities (torch.Tensor) : A tensor with measured detector intensities. + sim_intensities (torch.Tensor) : A tensor of simulated detector intensities + mask (torch.Tensor) : A mask with ones for pixels to include and zeros for pixels to exclude + + Returns: + loss (torch.Tensor) : A single value for the summed mse + + """ + if mask is None: + return t.sum((sim_intensities - intensities)**2) \ + / intensities.view(-1).shape[0] + else: + masked_intensities = intensities.masked_select(mask) + return t.sum((sim_intensities.masked_select(mask) - + masked_intensities)**2) \ + / masked_intensities.shape[0] + + + +def poisson_nll(intensities, sim_intensities, mask=None): + """ Returns the Poisson negative log likelihood for a simulated dataset's intensities + + Calculates the overall Poisson maximum likelihood metric using + diffraction intensities - the measured set of detector intensities - + and a simulated set of intensities. This loss would be appropriate + for detectors in a single-photon counting mode, with their output + scaled to number of photons + + Note that this calculation ignores the log(intensities!) term in the + full expression for Poisson negative log likelihood. This term doesn't + change the calculated gradients so isn't worth taking the time to compute + + It can accept intensity and simulated intensity tensors of any shape + as long as their shapes match, and the provided mask array can be + broadcast correctly along them. + + Args: + intensities (torch.Tensor) : A tensor with measured detector intensities. + sim_intensities (torch.Tensor) : A tensor of simulated detector intensities + mask (torch.Tensor) : A mask with ones for pixels to include and zeros for pixels to exclude + + Returns: + loss (torch.Tensor) : A single value for the poisson ML metric + + """ + if mask is None: + return t.sum(sim_intensities - + intensities * t.log(sim_intensities)) \ + / intensities.view(-1).shape[0] + + else: + masked_intensities = intensities.masked_select(mask) + masked_sims = sim_intensities.masked_select(mask) + return t.sum(masked_sims - masked_intensities * + t.log(masked_sims)) / masked_intensities.shape[0] diff --git a/build/lib/CDTools/tools/projectors.py b/build/lib/CDTools/tools/projectors.py new file mode 100644 index 0000000..b772888 --- /dev/null +++ b/build/lib/CDTools/tools/projectors.py @@ -0,0 +1,52 @@ +from __future__ import division, print_function, absolute_import +from CDTools.tools.cmath import * +import torch as t + +__all__ = ['modulus', 'support'] + +def modulus(wavefront, intensities, mask = None): + """Implements the modulus constraint in torch + + This accepts a torch tensor representing the propagated simulated wavefront(s), + where the last dimension represents the real and imaginary components of + the propagated wavefield(s). It projects the modulus of the diffraction pattern + onto the modulus of the simulated wavefield. + + It assumes that the wavefront is stored in an array + [i,j] where i corresponds to the y-axis and j corresponds to the + x-axis, with the origin following the CS standard of being in the + upper right. + + Args: + wavefront (torch.Tensor) : The JxNxMx2 stack of complex propagated wavefronts + intensities (torch.Tensor): The measured diffraction pattern(s) stored as an JxNxM stack of real tensors + mask (torch.Tensor) : Mask for the intensities array with shape JxNxM, where bad detector pixels are set to 0 and usable pixels set to 1 + Returns: + torch.Tensor : The JxNxMx2 propagated wavefield with corrected intensities + """ + if mask is None: + return cmult(cphase(wavefront), intensities**.5) + else: + return cmult(cphase(wavefront), intensities.masked_select(mask)**.5) + + +def support(wavefront, mask): + """Implements the support constraint in torch + + This accepts a torch tensor representing the propagated simulated wavefront(s), + where the last dimension represents the real and imaginary components of + the propagated wavefield(s). It projects the support of the imaged object + onto the simulated wavefront via a mask. + + It assumes that the wavefront is stored in an array + [i,j] where i corresponds to the y-axis and j corresponds to the + x-axis, with the origin following the CS standard of being in the + upper right. + + Args: + wavefront (torch.Tensor) : The JxNxMx2 stack of complex propagated wavefronts + mask (torch.Tensor) : Mask for the intensities array with shape JxNxM, where bad detector pixels are set to 0 and usable pixels set to 1 + Returns: + torch.Tensor : The JxNxMx2 wavefield with the mask applied + """ + return wavefront.masked_select(mask) diff --git a/build/lib/CDTools/tools/propagators.py b/build/lib/CDTools/tools/propagators.py new file mode 100644 index 0000000..d977310 --- /dev/null +++ b/build/lib/CDTools/tools/propagators.py @@ -0,0 +1,168 @@ +from __future__ import division, print_function, absolute_import +from CDTools.tools.cmath import * +import torch as t + +__all__ = ['far_field', 'near_field', 'inverse_far_field', 'inverse_near_field', 'get_exit_waves'] + + +def far_field(wavefront, detector_shape, detector_center=None, + scaling=1): + """Implements a far-field propagator in torch + + This accepts a torch tensor, where the last dimension + represents the real and imaginary components of the wavefield, + and returns the far-field propagated version of it using the provided + geometrical information about the detector. It assumes that the + propagation is purely far-field, without checking that the geometry + is consistent with that assumption. Note that the pitch of the + real space array is assumed to be consistent with the detector geometry, + such that the pixel spacing on the detector corresponds to the full + size covered by the wavefield array. + + It also assumes that the real space wavefront is stored in an array + [i,j] where i corresponds to the y-axis and j corresponds to the + x-axis, with the origin following the CS standard of being in the + upper right. + + Args: + wavefront (torch.Tensor) : The JxNxMx2 stack of complex wavefronts to be propagated + detector_shape (array_like): The shape of the detector to simulate + detector_center (array_like): Optional, the pixel (i,j) coordinates of the intersection of the detector with the forward propagation direction. + scaling (int) : Default is 1, the downscaling to apply to the measured diffraction pattern + Returns: + torch.Tensor : The Jxdetector_shapex2 propagated wavefield + """ + + if detector_center is None: + # Default is the exact center. This is subtly different from the + # default for a shifted FFT, where for even-sized arrays, the + # zero frequency pixel is placed at (shape-1)//2, not (shape-1)/2 + detector_center = (np.array(detector_shape)-np.array([1,1]))/2 + + center = np.array(detector_center) + # Split the center into a pixel and subpixel shift + int_center = np.floor(center).astype(int) + subpixel_shift = (center - int_center) + + # A selection pulling out the final area from the simulated + # diffraction pattern + # To be used as arr[sel[0]:sel[1],sel[2]:sel[3]] + # Use the wavefront shape before downsampling + wf_shape = tuple(wavefront.shape) + sel = ((wf_shape[-3]-1)//2 - scaling*int_center[0], + (wf_shape[-3]-1)//2 - scaling*int_center[0] + + detector_shape[0]*scaling, + (wf_shape[-2]-1)//2 - scaling*int_center[1], + (wf_shape[-2]-1)//2 - scaling*int_center[1] + + detector_shape[1]*scaling) + # This generates a phase ramp to use for the final subpixel + # shift. + Is, Js = np.mgrid[0:wavefront.shape[-3],0:wavefront.shape[-2]] + Is = Is - np.mean(Is) + Js = Js - np.mean(Js) + locs = np.stack((Is, Js), axis=-1) + # Move from pixel to frequency units + phase_ramp_freq = 2 * np.pi * subpixel_shift / wavefront.shape[-3:-1] + phase_ramp = np.exp(-1j * np.dot(locs, phase_ramp_freq)) + phase_ramp = complex_to_torch(phase_ramp).to(device=wavefront.device, + dtype=wavefront.dtype) + + ramped_wavefront = cmult(phase_ramp[None,...], wavefront) + + sims = fftshift(fft(ifftshift(ramped_wavefront, dims=(-2,-3)), + 2, normalized=True), dims=(-2,-3)) + + return sims[:,sel[0]:sel[1],sel[2]:sel[3]] + +def inverse_far_field(wavefront, detector_shape, detector_center=None, + scaling=1): + if detector_center is None: + # Default is the exact center. This is subtly different from the + # default for a shifted FFT, where for even-sized arrays, the + # zero frequency pixel is placed at (shape-1)//2, not (shape-1)/2 + detector_center = (np.array(detector_shape)-np.array([1,1]))/2 + + +def inverse_near_field(): + pass + + + +def near_field(wavefront, spacing, wavelength, z): + """Implements an angular-spectrum based near-field propagator in torch + + This function is an angular-spectrum based near field + propagator that will work on torch Tensors. The function is structured + this way - to generate the propagator first - because the + generation of the propagation mask is a bit expensive and if this + propagator is used in a reconstruction program, then it will be best + to calculate this mask once and close over it. + The resulting function accepts an 3d torch tensor, where the + last dimension represents the real and imaginary components of + the wavefield, and returns the near-field propagated version of it. + + Args: + wavefront (torch.Tensor) : The JxNxMx2 stack of complex wavefronts to be propagated + spacing (iterable) : The pixel size in each dimension of the arrays to be propagated + wavelength (float) : The wavelength of light to simulate propagation of + z (float) : The distance to simulate propagation over + Returns: + function : A function to propagate a torch tensor. + """ + + ki = fftpack.fftfreq(shape[0],spacing[0]) + kj = fftpack.fftfreq(shape[1],spacing[1]) + Ki, Kj = np.meshgrid(ki,kj) + propagator = np.exp(1j*np.sqrt((2*np.pi/wavelength)**2 + - Ki**2 - Kj**2) * z) + propagator = complex_to_float(propagator).astype(np.float32) + propagator = t.from_numpy(propagator).cuda() + + return t.ifft(propagator * t.fft(wavefront,2),2) + + + +def get_exit_waves(probe, object, translations): + """Returns a stack of exit waves accounting for subpixel shifts + + This function returns a collection of exit waves, with the first + dimension as the translation index and the final dimensions + corresponding to the detector. The exit waves are calculated by + shifting the object with each translation in turn, using linear + interpolation. + Args: + probe (torch.Tensor) : An MxM probe function for the exit waves + object (torch.Tensor) : The object function to be probed + translations (torch.Tensor) : The Nx2 array of translations to simulate + Returns: + torch.Tensor : An NxMxM tensor of the calculated exit waves + """ + + # Separate the translations into a part that chooses the window + # And a part that defines the windowing function + integer_translations = t.floor(translations) + subpixel_translations = translations - integer_translations + integer_translations = integer_translations.to(dtype=t.int32) + + selections = [] + for tr, sp in zip(integer_translations, + subpixel_translations): + + sel00 = object[tr[0]:tr[0]+probe.shape[0], + tr[1]:tr[1]+probe.shape[1]] + + sel01 = object[tr[0]:tr[0]+probe.shape[0], + tr[1]+1:tr[1]+1+probe.shape[1]] + + sel10 = object[tr[0]+1:tr[0]+1+probe.shape[0], + tr[1]:tr[1]+probe.shape[1]] + + sel11 = object[tr[0]+1:tr[0]+1+probe.shape[0], + tr[1]+1:tr[1]+1+probe.shape[1]] + + selections.append(sel00 * (1-sp[0])*(1-sp[1]) + \ + sel01 * (1-sp[0])*sp[1] + \ + sel10 * sp[0]*(1-sp[1]) + \ + sel11 * sp[0]*sp[1]) + + return t.stack([cmult(probe,selection) for selection in selections]) diff --git a/dist/CDTools-0.0.1-py2.7.egg b/dist/CDTools-0.0.1-py2.7.egg new file mode 100644 index 0000000000000000000000000000000000000000..443ee231d6d546a6cca410ff6184815227f62ec7 GIT binary patch literal 26496 zcma&NQ;=xUvMpG)xy!a~+qP}nwr$(CZQHKeWm~;ZM8~_g-)%*#{8(R^F>=nFiIMVB zz#u39000mGTyoB`WTsv^Qvd$6{|WLx5fWCmv$J-h)6=uCwQ$zcqqX;tqMDpj36GPb zqM@LXjEPelpQah7qf`)|pq`?hot|4#ASa)krUINiBT4grkB}I7=O`hmxJO9Q*0#V$ z&ejopcyx+1*p%YOrUj7yU%KJ*J!D7Sp`FkF>0|>00QiqZ|2MjgRO03&hv^YU$tG0O zK%YUq{dHS+Kq+vAK}L$gE-E07UgD(iEhL!|o>{zw^#;)GMPCfXeW#^Qv*;N3%&Kn3 z>ek?5r5-!Hln~gkpSn@m0qXJ<_F0%Qy8}l36S-h*0eOQ+kGA`rdi4`h7&;8*6;0RS zv26aFXD3;-5*>J9awRT;rn%a{&@T9|${CbuN6;e&?;k4p5*|bbR2MPBhqtgVFRD*D zohLi|-VuJcyq(^`HyAcgboN?%VO899(GCxJc2QJ5A9N00RARq;KCl-GHy4o%8-9F! z4#dFwI1bhP5oPFJ?)xujlGL4LDdR?k)&D^f^iMGU8#K=UyK{V^rJ(%jVL)GbMBwms z1>^R4mp2XbZ6p&tef`>99`_lBh!23mE40$WA(Szr21LOr9B}d>Rl3|^U%dM$&ukI( zmO7l#OFax`IX8Rsh*9i5;peI`S91Q#od0-APSa!rqjdB_%sd4-=Zh z3&y+%0KN~fBMpsIsUAHqUT>;^5Cui?G=e~3QMoBH1j-Jeb^^4Xvw=X;_pmi##`UE` z)2HXFi^X2+KldU7X#2kSpsEj;9jF)poe$F=A|IMhjP4)PKO#R0jyLNIcoOIXA}2A8 zK(tP5T%cr-DHk;Q>O5-tM@H_$=K?Eno+GRhvnQu~X%yDCSa?G0fVs9XE=Z?dJ#PKL z$dNFBB+3*{#;DAN02SIpW1PL&u*1*&W!F%$FvM9&coF>zdAoi&sli6I4UED?<;bD0 zw^#D<`GushaumX#$w zXO?ni>{B<4Z^i(n3z}sA=^obg!w3`i5pYS57 z)NMNe^8D<*tcE^*@7DnA351|Zp(z=|?OB?P*e|}H_u;{n5x@OxK~j8^UGa=4FNn%! zpDlKB1os`1^lE68*`r+WIw=V7Dio+643cw!C`0K3DQL|B4YI~e+|2Yk0TL}yG$=~O zz;fdn0EV`}9Iy{!+~v$$5QmjD1}X602#>MU@@52`!1?M#YE_bb_=}mfHVLv~jFwX# z1UNf8+gQ?M)u2xOInULio{t#Xz3{pjlqTpbn^iCF6nH>ox5Zly`Z`@GkR63yDkVel zMrf?jsy^Nd`s;Lm&kdvt(v{D9Hp0QvK=1=5!xDXbdkD>opO`lrC4Kcfd+D+l`Us>8 zK+n^a@&vV)I<>FnDbaUva0Hg$J+2!LPj!Mi?ots`I8HB<3j|J0goT2c2ponB!&HwD zc4xNrl_@X4?U#Zh?-6P6z)65lRJ2zpu3~;l4NS7T$j zeqJEFAD9NUY8QU0Q~-~YU#^Q`D3N=ain5rg1Ks&Jvh?=W^W2?zS-koN21`qfD|je2+@JI=Y6F2FDux841* zJ?=nmT@H;IBjx+y1T|bs)L}@i_3W8IHIFC`*bg4Pl}v0Om@1*q$$<4(N7KmN=>frn zBDj6#+^#{Kk7qRzQ05pTVbopDi3D(Wg}Ov7p)pqN&yUzK0UaNZ8}O*v@E*ep_(SZg zcf`;k;X)=lg53742?e?XLEc26XSTe0epWABGq>&4!25ZhA!DHA7Y^VjDbNWzq%=97p>`|KdV~D4&$GLzenIto2fg=QK$v z&lL1&9t$KnG<#tL_e00Z9mu_*P))icibBO2P=2Ztls(NF(40zLdA+&-SRZ znJinEd}L=MkEe2{dpOv-a(%(z#TP)0!jM?ngWY=)Y&?!75p6)@2Nfaimul+aW4e>8 zkwprKstUIVa=;jXvLy4m72hJ0k8Vaq$@ zcFh=E(XZi!BNWrDRIq?T*hAh${D!eINSRL2vqq>i@1Fv-j-*|wrXZkcn&*_T2Y}+b zArT35uh`hes5mKQ?}X$nZ|F}zpXuYWhA6YB6M&Y&)6Bz4sx^dK3zl-{xmMyOqVNoWb0^rEUH(BBY0hH23kF8l*Il7ZPj<_M*M z;kz;h%t77g3X8lQ5YoaAA=&ytAn^``ANFa6kVpZTUgga|;3hpa>rT0R8{*Q8T*M(6htlLVdlmO*GB&vsomiK(~YnuQ6|v z*HN(#@t3Niy5|pJO^Gk`YFXoYG0mBsl^)HtZneyHM7-_B2Rux$>qqDBzMljaI0WwJ z7(3e7n8|A;Z)*;(VdQyP9~U5N2do^u1f^4TNC$ZT@zO~ec*(olCDjvpF9SGw36P%Kfr{ev}>dTYVDGJ zAf3P>zin>S$2_$IG>S`Fhr+?yMHUU~&KoOOflU6C7%f3ch<$w%-vWeBGQg}L!aWH**j{ib_s*CS&Uo8Og&!_B?9mo`;caI1$Ry1mnPHl|5ulp7=4bgF z)Q{p`qC8?nZU{IXx-i%NlLU{7zFkbY?9+^y8BQ@GlVD;9<9W@G#Vc)$OOcCZl0Voo z499u^=3J@awE1P0T2;gx05@|^0-5wCZ-S;)Qg$U?(Z0!n_v~otEdqC%Z2Xla_8Egq--)6Q`_KKZdhf^BtT+m3J}9@vvOcihcEeDz&gKF z%TF-lY~AHc>55t@JV&fxtAdUaq?NnLQcewayICwXS&>LwO}V@2VGa2`6JUc3~Ydsv3je;H1A`DUXbrY0m9DsB?KDh(Bu;y#MFvMuVwV<;>o%$ zH7ct;qfPK?HDb#R2}FJ~Pxxm&^XiNcR4@C4(788VbHFrd=oTGx5?g(Hu@e^GU2lx= zP%#MI-o#_^3|>L5bWs|^9@c^#wA3tkz(VSZUF9rYNOF7JKapHgG@KQPv4Bn#ozsYK z_GMa(_~16YO?FMkJ5w$EN*=Xb)hs(6TmH50+!y1(LSrJ+RCf+Id z9Y1RVyr&}87H%XxgV<}|>&{Ejr`Vp~vd_~2a7)|k{(X{}1X2cgYA(S|kt^p7J}IjL zEZm%`WHSr*=1r=(Q#q$8Ie?zsk>8^=*1DI`(#k*0@+%S|kn_RyxM{xrkt=0)sv|H2wv7PjK{T>kaH;DRofO9?=0n6m6lD_B z6W(T(u}j_aTcnFif5#ZIKE98}Xq=m^>xU+WJ(Kl{w z@&a{dj&8Br6*`KU94r8Y?tY$NIIFzuYt-Q7#>|~r4T&BN!PIVv5+y-$HKHoR&g zR1|au4X%Zy)^!*WA-Ba@q9Dtvy|7#`oIp)0X!p>5We+KJVTF~;&VB#pxC@?FKP@O} zw@RMg?)uegr;dIw=Qdf#8c^s)itRPIjV|9K2T1qEdCToJ?hQnUl9(NAZHigkNax97 zqzhlr6{v%IhK=)tI`+)DRM+jB*71jVw>rIrv2qLK0=x$4!>VjsPBR{piag&hr+pxb zDHK&aQFAqiJC!@j|G<~44h`{b)4#6y*c*MWhR@anH#g=10=_TgI$%&r9kT~ixCn6^ zFq{Eybh>0z2q*e3I6hPn#0PR5Ai9ei9|1%bsc7o~I>sRF;*ymY!4_&gKc%;4nv&z;cn5zZw+Y6)>f+<5^S~_A;VX(K|-8j&qc~kT*G+ z3|(qzanxanHwJGw*g15+5;+|gXCzqt2;O{C=4xhfE>UWxHhWo=tG0vCdb$VX* z4A|Qn6;JQ6!THP;v@r4+t!?cjdI)8fG$oa#=dcc=(*4!@rqp4sr&MIfJj(ua#YcN?jdca zV-rlTc?a4Qi%>0*64_&K@WhhX0$yj`GT-0;Os5I&8@@VJ>P6K+Ovs2;RNFrMuZ~`K zLwk?&=iYx6`Q!c`vWaHc!ufCj0PFMs0I2^XA2K#@HuzWL&uH4jVYMOt`1S>|`^AHd zKj`R24*M&x_Y%`X5sMxbg2ILNZ$TnPT~9BuGO?kOINX}(ii#q8%5$f2?W`1D z5^!FjG-#->$6Bgkn-!D_1)G*-C!!>h67Yg3U(`=*m9NrhbL`UBeEYo?jZj5hA8f&? z+Ju6(RH0yv$BaS=zmV+WNzks+Y+0lrwMuEiC-!ADwBuB`tJJAtpcWGbZ7#F&)=Lv^ zMB9wOnd?)^$G@c%Ij-4tRUw8q0zg!|9^idMG_RJP`ZVx0L&a$Wi8$~TP${V(lPw1q z0de+`GtO_)tXN$9;5tlRF1}>p&5Dc|ql*-yQSDSjYrat?uwawp*DBniRR41Y8#G&? zy8px$O1y~`{_bR^Q`kem5RuLYt~Xi=>Qo7y%y((>nR9lTUgZ9LCcy$pFYXm}5;ZD6soV=VNe^p% zMYa)|wXr?T>ZGN-XY_pe==~mSQB%oW@elcR^H%@ii~6Qu%w0X}RrM zBrk3KK#*$J{>4vn;iF*mlr9~As%cZ-6ZJe2gV^+>oRPHo3SExv1X?j($cSmShk5UWdQ3Wh z7|LD0mq5uaJi;a=-Rf%PMLbSHy;~^egM#f9L>G}SIJbfp8*@93^jq|~C?BCl4oK1x|r$T8k@GT`O zout@&3rKwu7U`yN@v*kra69_;-K4P_eQy5`*VDMv!!gShNx#w<;_zBX3ol<|lc4jV z@CLgz?L-sKwF@p|Oee-KsVf=@skASK8v#IbLE7{0GJL01TU=OnO%90{d@UPJpl&gB z`9{5UVOp(tU}-I~OTz5&X|eOPcrrdLC~)y)PS#IWuq%YP)zPu^-8`Fxi%IRuI7R%;X8ZO`fUK(!XGB8&diYHs zGwRTc4HT4AW~-O%roySIslKKdbza$1k=D;+!`wJ`QC#N2z;L+Xr{y~2Cg}{n9_zIi zk38Rzo{?n1$3redA?~qJj4c8S|I#UF-9PJ6%HdMgo*w--;VP{5aFf8HbcNerZJnYY znD3oAAPVs?i~7VZ#1qtK9kN{G+~eer$BYU7Z#Xkg8tbQWW_3JsxdZS`dgShmXx4-MK_&^m!&34;@caeD z%{&*t5!(`Gcc$zZ={LnyEm~+$DF>Jh-5gf>FKx{9STNCkHSw^W(JI2xfgnj`IM#uf zi@ARsEw|(~LWM7uVTcW0TR0^MDoy}ZTyNbw$}#57hrxh&@#14E?>U3hU7GYU^ZZ!n zudWK&QgWjv()H!2vyjavTT{FW)8S{(78P}wWU@ki;+cQkABkXksYhJlwW=w0QyXa# zfSs=;(aKY$yHZw690h1wmhMF`>x_<5Ofo)2uRz|k&jA%WrAbefK!Szk&SQb>iVh^e zJVZ8+m=Q{BhcU1TZ5M~rOEs2SMkMxN(;3W$ern{)$ejGS&&}(9uF3AWf-{fk3pnb| z{-6~EQqLfdy02~$Ja9I3MCSzDTc&R|Jl7j*mR1gc@z!zS-}Lbcu~$P2!@qUtdVjCz z^}ZcW)vYYe^`guH_4NDCFsN7(d;2MU+HI{iFHOKpDuI&yWKeXW?s>Y}CKe)2g54`t0Mf{m=CN|^x!{#l&L zfC4xjYG%h&$OtqFDq~)j)|+q~Z#@j?YpDpHB%>2&LG+Q1aReOE*C>H(U8C!d`mif+ zYM<3<>pbS}6ChC95McyUlsL}xxR<6QS_<~Z0V(B$Vcl5n+!qW1h-UC^Xt;}!i@CUT6WP)I;bsN{;I)eJ zk53Zq1lQH;4d_Lt>DlJ@mv>C5y0j*i_#{vI}$Uz8Ty#cee& zhDSNyG6T+l2k>{N8ym=yXV>88C9Na)ADW$&m0jk=-lMuzskBk-i;K+ zu!RG6doeB-D2A5tak`!0{!Nhg#Zxu|)Gh$eI2b$@Ffpz!EbU;}J?Z|baN2Ky1ExJI z6rx$~tPMR-{Twk|pbKjao7g`}bklv94ayh0m=JpmFT;giX>UsY&^-9t`q{K~!)Lhl zonZ4w?cR?zn?OomeH~$c?*oW>n0vZIErZElfNi z9hW{^$gT$}ZF=igvb9z9zLXvEzS5<{6a01da(<+bm|7sX&%f{|s<5`q!9ib?M;A9; z4*UGu!`;JCEC_Bb$$fF%HpjvNQ#TP1p z>IMaj@j%Jflb1e?W5kys33EX%8O(qk@fboY3`J$D8D0;w(FPO~X&{H%8nJ-0Y<=re zokI+#e#~9(_1%eF)9Bh+KeBC(Z+>^!eZwV`_Gl;f-ZiDxq(yh$ZD}ET15#@`ZQsaZ z+nM;ryL-dh`206%Q&t7dtPuiA&?rf&EZayx*Q6NRIDhRueWn_$0Q&V+d5p)2`86mK z%1;7I^Kkh@?B~MzBDzL=j8<1Dw_Ht$W}tva_iy?~hDJ?#n!6LB&#rXCE9ae<+KzDW zEig~MtPBtIE6(>*6*CT{6q8z&#yIVP3CikiGtu`_j?{iiihx_G3d@<)M|_2;2;h=2 zrLxP~3h$~OkM6-4re8PRJM`c*OnzETVoP0acVA> zcwU0rXrD~`-LC3(-MeY*6jRmVuC_fs9^4v&Gwd+1pSWF6KsIMp%iV7?eXDwGLni!Y ze-O$Gawn!n<0PHzCr^57m|2 zV!c6lDpk(qPVPDttBu09Hr5wZf9~jmX>JcV5a^b*LUhxioPe6m-HO#ND^xnt4i$-J zh^;`7pI$6vBNj7spQO-|x!Qp5BhXbg<@LSR%^2&!8k0r))qwEiH_IA2BI^C35)pSp>At~{{aW@BP+2XBp4Wgaj;BQ##qtR`pP z{ItnioKQ&Wq2ROry(N0+Zj!iywXdwh*ndxp-WgC*YqHgz(yJMNx5`|wYx;4ST&2wSJBT5mKoy#?7&yxeXA4m< z5pmsK+(6TNRA~i_)jtlsdB2noD%J==+Smdr6J?Lx&bLpjb|c8~p-sSd$}pTy{Pm=|mrKi8zjz)^=mD09cXU)4 zwjkt*-YJKAOYU(i7R!F4Om@(TK&FVs)TALby_qS5-vP9zB7Q5Rrh8FIDdf zmm;7p%GljWFg06U>PciJ(V<6RSKxQZZy&^qhgQL6x5pRz0j{+dKjBxo~>J`PWYMLe|&(aF7S;AjpLIF7KEW`u#G67 zS+!xP4u6iUH=jl`P^nfE*{qFt?!2QK6qA#AGoOR_4Zg!F3TG+BS#q zOU05om~%KPzLRq>Qhj^ePDlQ+u+1(cv2DW*QU`?GtlTeEW0yF21&Si z8BHjwu)@{335^_TNZrwAk3{x~I{#2~Gbv7zU<>zTSW-xH3@eh)?2Wb_+es~m{_xD7 zE79p)UR-imyqI)puDb1?80p2k!21D*1X=*06d?r#K#3# z+h1^}PUbNS*j~oWW<4>dg^qxb4yN7p;7V*9F04g~;!8FULrUF&kmBKLHMj^ber3Iv zK=BV*1Q|%{c1y-RWqTbWU`S{3_x0l&anb>Vpq>6ProW(!3(!kMy~?^s(d*8p8<~f` zhj3j9QO{Y$n%AVARRS|%)I@oGcLJMl^OMn(|mY z=m+fRY%~sf;h)kER9f3rM<@|GQ8n4=7f^NC=8VXCo(*&|d_0>}osy!lN{1HP6lxx= zov~x1>=554uzd2JJkCAp*Su|=i@h(~H_ui(R&=yezJBj3EnE)Rl_oQCicVh43e9_2 z&X~H~*ga1;hhkJL!XuV4UizkhSBx7HE)v8eM;d7vqKdAwJeBS^n%^dej-l*JP^+pHct){|`nf+%kz}_`_6dmd|nvfSmmaemQ&^8e=S#h0_MY+i6=kYclu8+H0Kz??B5K0$(OfDh4g-PUWRTRKX&^>HK z2HDxm)v_96XU9_%H$FGrJZT#c?#-`NOChgkT!y#$zdRGWO6=?gKkQJ>Ye@#{*JStZ zPv;*FVS#FA(m&7CJh*e(fjDY7eMggxM;9E9YwQo=GT6cDeBqjr&qXH)`hO3bi%aDT z!2{Y*C8Be47gt#cFXrQqbK0h5FoysK zAP3JFgsff11WgbHo0k&NTX;MO!sZLi*pvxkJ)a}O58@_Wg66Sg?!cuxsi)DD?u0`5x&1lRqLsR8_ zj1fcvu2HlkG&dhZD(g_8C~FlgrL>;`bt~@8g2=?R|*O~I1epy5WfOW>3DZ_&{Adpd?SOD|f{EuQZRwpkW?~{Pme7s~^ zR-Wbd^y?la;a@{c!MNs8?&>9YElsv2X2TH;wNdj;g!BE$GtH|D{~YYYq_|?zJxvWw z4Q8-(4$Il6zwWaKjC3{CmC#ChX=N`%5m(OeI_qv--S%n&*BwTR7xF0fYm7A*kQI5? zg~961x;C_oCC?;9y1n$p7k@?~8%|aE8B~3Wl9jL>_>+x`O*0tqo4DJZt1{01s5J@~ zJ&JfzCrW3hB>@M79b0co85;jh)(I9k4)t+5i%23+gZol<(iksZn=?I+=-aeV8e`;9 z{o5cJcV^+;TkTE0gB_W2eV|W5fG6Zc4iW zK6;5>9rdPDPaEXe+kU0iBz<-UKQXh%m7GPisiY>7&YiQV0(Diobmkj5m?9*!VqHGT z!xCGlav|Zd8nMs1Nb-}8^5gd643&oksRGBRu-(n%KVyt4yRHNh4mNqJc7mIe=DIBI zseJ5duc#9&)uXM#>!TRag2zAWfg-i~U)F*K%l{AkZn? zs~8V#12w5@#{eqN&jK11!=HrLpUI@5fw2x^h#+4BaHP2XuU8TBPK)XS-*Gv?y_zKs z5juf(>hQ*|KJ=BcgZfo{*fgx|YQKV$yo2b8A<4!F33Si$08g(5>hqRyWQ)(eD8IB{ z(8u0x6o zlN!)Za;$ne%Hd4_xrHB4vTR;)GikxicZjxx~OedPK64RI|Bu!t(7N% zwbL9pWRoqimObGxfI{=gb|`pW&n|^tRaQ-Q=T!la)s1i$FisWapYF5__IDpssNIM& zuY}cQ7K!nXvcQWCrMFjUq;npXf2Ne5MNdO!44_W<44!zJ2)fRKJt(~knKp>@1c3Ys zdMT*~L5;+rCqUs0DNmnWlvm%Wl}^r$xw*VTqu^T17R_echZtnh_Oz|z>qX2#6omOM z+PfO!acuc4UVz@!hBLhVGCM=WDt{mgcbvB}(s&su z;4)$ne`}N;Mn>uxl*w`e=}a+39m3Ccf?&l8YBTU-2L!W&36q8f<~obBEDZZqp5!6e zeWD|f!5h`t)#4O~-XM@uTm*wuiV9a>z{4G~>I*^YBD->X@Wd&bNeZwV@rPC^Oj?sy z3^F!l1|y&yatu0NGb;4}9@0*{+P{Ev)k<by9tdiyWpce1TK*K`utr|NhNMM4^*10Xe^?p@kn|5OY%yge7`SJYC@(M& z8NqgpBizqWMi<@;3hyeMJW^r1WTzX2p7YLE0I#pj+4fq5Cr=vNUucGk8_O|Aghm$& zv)GI*NUg&gO0z90c?Cz1L6%DN3Ct>OTF%E2lj-g?X-@psyZT6oWtC}zf*deNxqQ(pwX;L`GmmjJ-FuzV z+x~@Te{TP0HCokGEe}q>|7)dpY07N+^v?$Qf9!6;|39`E@x)5W4$&ioj!8h5A36ZS z#mgYa%P)fh5qFv@sB4Im!Xs0^BaRp;$1XZK0-6>=P5nI0x#Rl%cr|cxj9zgZ!H4qvNO9FO%WRX>_A;uwdb`QE?w>oNfstbL%S7 z%Z?U(NcRfK>L^rN7}Z0zFxtW(_o8eEQw0_ZN++pWjuE*94)8ii+)c#rjThjZfrM6k z3epZ5w%*0HlQd7P_gZ=(H#p0nANu-uDXm_9m#RMNl!O&{?)HZ08-*~%Sb^1yQaE#p z9)%^YHne%6e;7Xl+eJbTwy>l7{TF6Fc^R_5v?PAlgSP#inwpyD$r=V$A|$jy zCwb_s1;CRuag3g5Xb~WxKU$Bt+fLMa>0r&|h{)u6g-T-&c@uq?5jBBUUBt#j)9ijg ztqfgBGz^1cKPd1hcL9x@OO9y*qgi7-brQji5Aa%Anfgi|R3vT8-`;_lF!UHUIbCVW z{z*3*9_<= zdop4R+81hBeTbPleAcjX2#@^O-a30dCf07uf>>fAO*5OrB-gdhy6FJ6q@(pqtz z^i$B#L%RPuL;vZpMEhx;ef~Q`U!VX0xc^^g$f!nnz!sSS;hB^oS)O0;(0_PfOy-=_ zCWKNM5NIHUFeG@WV~P! z`f?wiIx~FHi$r>JaJDAow5|A^;egQ!J$bhmVFH+FfApu=eAVXovcpF&sfH8N% ziA0>Y%}BU903wXEvrjP(T8B`?GAF6Go>a+naDWg%d!1ORgZVgdo&{fT6D9t<>zGch zlJ%rIJjp#(S&MU%_uSg{=A}FH7>H}#7jrIcBByEsN75DOA^vRT?tJfr~SR!d`D&Ip%CYELlUOyL))Rap=&8$Oe}}t3)ezw2qQxtMjOY|)S2hYKB%9ZBe(f91QWYpv zDa0ISnlGmVWz@L@FlvUHG#H_0uT?Wq{)(y|)#6Y+@=(ZgyZ77P^9zh;LyY%+JfC|F|3-ESH1+XrErM#V{$ z{Ezv*L^Gld^9}QV=N|uuYqO$rVYBov?9Tt&SNosbV{1DnCzJmK-D|bIe;YLjzH538 zkGR(04QtO0{7}w5l=Y}X^KG>V=I6%9r|>DoSH*tcv0@X6*d7v{g(Lodz_R1JbOyWR01(*l>X_fm&dTGRRE>yjA$8j5N-}(M+N^R#g3R;W6@!ExTvU_ zf;^4?kAL2T!znm=*?b?CbJc_#Db^*hw8#*MuLwO+@1tz>1`92qAq#B*{DFzJA1(>F zKr55=4c3-S25yYM^cR+s)H4nNSh7RhYwFU0c|izZIy#MkWlSUTJ-as_}{3`RAm#CA019uDOX{OmGXvi;wCoROC`h6 zx%g^O^qCLfvDRR9ANn|T+M>Rn{F8nVdZdX&tv9~6_~yBe8# zL#2?-dkX3Etg(Tf5Rg>Yi$Dx&IuZ6Uo)1&kMk4E#7xF#xh@G~ARqV$^PBONlGM5&T zi|__ADAZ_Aj~gK)8c@^*E}DEl&6iZ33k-1B z5jO@LligW};$8hkh1^=hzo~MOhEV$h<2PO|PBt`=+?}=u)^tyN{FMjk- z+qR(-GbBB^z9{uwbaC0QB_5IGVveUpLafk;1-Gu@a118Z;mX|&1U!&d1_E7@W}16W zRI^0pEs#!W`be?^61;MVB~)I|#4D}S-L3iElXuPKbLOUc7H2of=*0yXxz92ZiKwL^ zFJv7b6>owh5pR~pU-)eMDPW=$w4MyIvZx$-tz#cKpQ7LzWRvbek2_)x2ZwIyk^GTf zFOQ`l9HaV>m0C1Eb3Dk_y501WaJ(Qy{<{Y~n zk4}z?q7@>}Ic6TA^p#{WwIu?wP_WghFx%_j|80FXi1AM#$iKxxo_}%h&;Ljt{#P6{ zx=!D-!xBLld-Zl*bbe$ zkVxSjfWE!!&l~I;rt$-rp5A$%l?9>^78h?heDr)h?P6-`!kZiU{g(V)^Xl`)YPdKV z`b@gzHYH&;9DO;sXlNfoiYvis-7~0|d+|wi-SAEGU_Wgt9(!3+-!^fG8?my6a;rA- z>z?(!TmLb!)GKt&+f1wI3AiR3;+;0G^N7B_n1)Qz;6)cbAAAtjnj#ZH<0vWEvfxea zu4g<&wrQWukz9ce2>Sj(kDV&g)`EPUiB2}AzQ(C#nuu`@nNVg8VVF>&V6@FcAM{TY ze-!Q@Ox4CLg5BA3V>+ND^4b(Oku_EAINB68W{_#JT$LM1kc>(uyWl6aEOmjAoGlLr zLfFRDB2)4$-~LDINom0}bwEgPp=>Pj)jP<9aZuPM(MDiVLmPCl8%W#sCo;&vZG+S? zw02Qs$Eajp8@=wJTGtT!Q>v;sM+k7S!C&L}KrfXDVRor{tFsNoCv;RCbm?UukujI? zxZVV(69S%CpVc7cj)lZr5NQ(VnZ{YCwZ@kk5oB%YWQv)AZH4cXqdlDk034L&cJMBR>iA zRB8#AmmlK;pwQ1H5HAP>(rIhQ9}f> zkv^i{+b`tQBqKS2LIFZ2sl}&WCHc)Lm?Vqx1 zk^sJpQ>a<$U_b5JP%tj>^J4K-M3qb)R;AB)%{NOjtG~?=rbT5aTWEgccxn!*Am*-d z&ihEtb0Ui3UGss$A$6tT7FO3Unzbff^Tl7;1j66PXz??Pq!TL%(hx(GQAon1rk^3q zo-W$B3+@?l7v`)K1=%?aB(}N2sRBv(q+4hfkyGYWEGq*3-+*+6N4xeo2jSJ5z+wTI zcQXWo{qddGPDkYoY7g^>4EAFP5DelNVXCC@`2@Sd%+@Cn^__`sib6_mEeqx`@Q|DC z4e$qyw8}`g(MJ9s);@ZD9rs%A@FCp7riYw^ClAm?$E&AA9+TjDQP#BN&RH#LOs!&aODoi+e&a)9e! zSO9IW#~L8DYroI@-nf2kR=<-$su|fX{BN!iSxfz-QBBWcQ-8=p)&$=Z*MHNp(Zl;o zg1Z3^^!)8%YXNl-A^=yT09RK532FEP5e|Sr^!!1M{mmhFYYkioh6gcl^WkLsZDG1E zcVY2YRqTq4RqH&Fn|!O&Pmr5%X&9<>baneSmXmaKd9{;Kj?#l?)_S<{o^7fl^leAlK%YVq7l>=@2vU zdPiX20l?$x3o-uj;ei&F&TRK+d(3cEAHN`UQb#G}Om&-dQeG58r=`QLP}e7Z12`+8 z=2=Ku?48`sbyb0tRI!dxo*ZNjnrxaTr=F=t9RSzj2TUr|dZQQpL5sI{mj!|LjwSs4 zH4Ulaj#sQMnXyfi6=O<}A(`po_ns~j9hNkmb%b{A&m!;u3qk{bek9{=m4snkq%h|1 z9r(hUW=pTj9e!8SQ!VQ&{FWT)0`jbMH`uEM>T`p;tC+5-(o50RpNM4g4dP^n^>W+S z8O;OmTiSX8+eOH6{Oj^(loLv5t?Z@PKY4;4W}dVkleb@wt#XyZ4#mO*Vmg?M!Z$k8 z*U@~FG(0Li_~+n(LBC4;xahePLho_+5$vawZyocF@4q#C8bUl|&lrhJ|Bu4X0w}I# zTidupa0%}2?ykWd1{ol@yA#~qg1ZNTyGzhOa0%`f+$G5+=ltIx;pW_bS53{Xn(1fW z-m|M$@6~IqZhd<2B!mS4dCE%vPvM1sR~LA4J6RYTIobZpb#bI7XS>LR(EgwS(@sX+ z2E?s;hUz?@tb zax07MI<~P`R!gL0#7dfsxEI+z#A9t9Z3$tKw9hYzmzDPQ znYsorj3tC0{7m8cpd)7~oD+ZPTpw`GM(aE%Es>NIg@0T6#o_e=TqrihIaZ+uf8II* z3*32-`iZfPXkb$J(V5ewFpRh4SK4c#Em1wno-9^(3W|o3gn7#WFVFSuMjOY~Eiygm z2M=MGk~FW(o%$1^4WDefXw2d5?ei%S)Y7Qc08JKHtFqAG+p#q_b$`U`rCF2&Y;O^KW8c=m0LW2YL7ucg3dFMf{Q*JfVLeT+?LC;qvf!tQ3v#vw^- z)H2(8{bhN1gYr*jOY3(zn^kZSkVQBU5WN5Uo;NyDAGKMOK-N?KA+1`>&n4uS7eeo0 z00)o{?}G7Xvt9!ehwFnL5(2FcY9Z0s8Xh71hGhLT=DUZr>1OBfUqbS*7LZBKRy3nm z(0zR0%Dd`y(W~(7`&G`*UCtAqF!3AMv-{&lq1$Y`_%%*s^-w=2TZFna}v zIxM(W@-eenG#5amt5^+rymj10x30FVOCNx(oiSebOg+>-8&h(EzmH>!{$L%%DLn;# zQ%>M;-fK*wBYbVf0vHME%Sjh94s%1r7DP&6Tt~XlPuz z0Jls`=C#z{yWqos4IW$vBK^c!zNGs;#8LSRGE2`@`ljnnhAL6%NLQ z%8)#9D_iy33xzG?tRZ1xkNM#JZJ+C=tbj?-sn}bh3p$_;TN&!=Wh8MF8l}{Uts+?| zSP$I6IU!F+{K_D}*ewThZK{e6hZZ%+5nR7m&34g+rEpxja5{xk zL0hGxoDcfxzVRBW${w?$S;&yut+Tvt_41b;;1A39$N>}Fk%Zu~$VJ9>&LywQKuioe z46=n^3pf%nBEGd!X{Q9AECgRBE5EkvVMTxFG?NFn81-Xk#LStVOz2Vxm zHjC+CUDD8MCx4Wbx^J+Gmk$t^l9_|G3u&}WmbenSKsE&?Q?{)kbRKTS8!tI2AVvyN z!v;G(l0;qV_2V17qWbEg2{vW|sOIoDt9SzLd|U!nS@;Uk*#!G+{5*F*%IPWu!3bJ+ z(ciF_d<+g8-jKTN>uq0~VnX8WeuV`l?WJ_HIzs@!$I29rG@A^wWQk%Y8)|ZC(ns^!?g_no%Tka@Je4b2r2sY|f`sjWybH&I$&Zw#|H*f? zg96Q%ULY%Du4Tg{g$;aBO`vmQtd<6Z|sr&e=Ymvu?Z??|^1DC(MJUgp%Bj=ML zJa`%v{yH$&8JHUUV-2<5#4K{5c06bx4-vwITW?Q>9yL zO~$G&E_eErg(r&LB@IE;7F$F3gk62SFAJwuZ8tL4t}9X_kxN}-uV|3^kfLpj+t1FLAf0wz17YrJ zb(gzFLtj93@e1WZ5%y<#xnAHT3H{i2E`b(I3JhrTN&a^z_q|hZS%$>BJ%&NvLF>|@ zRqJ5K`SFyloJzQ@vMS-^9(T+pg;IQ*%pP;2{&1y>%CI-y>UwhQDo5=iRDmOz${G*b z@d3AErz>q}m=AAIIaBCCRF5;;d8DTjV%Ok1eUlUw(=Pu6Mir$ zBMTO{y>^5LWn}ysGE-Fv1_VLaZJ7-s{)rB5ttWW(+>v~aqP#+wBqiE@7d2vm!5z;Y z#Br=iz_*&H!3K52H!ET9$OVMbXVxeYZP>)#zKS1*S{V3dY27EZA_cD3DB#cm$?Me1 z+TEo5h~ciC=?H@Y2xCQsLi2PkfA5@b?WF~6;GPc&JP2j^fB}J#mu!hHCTu=3c^+p= zg|z{yU9L9(TJO=%Q99mR_)%f$+v1D{=+`1)4yK^#t)y)+uVfYv9|tX}FZ8;_fIhWj^?Z%Y`~4;^~XWP**jO$phzu+77M1jy#-U((I{l*R;Vr1BKXrmxU!_A&PrgJ zU+Zl?hT8SPMr5dlhyfVN7Lz+LIWSvblvI7@2jfM=kkf#FS(vm5)3q_zAt7A?l{;J< z;@Q5m0`*b8ong|C*~CeIsoh<;sV?&v7!V|bw9uE_lk+J1Aaz$q6(hY5QQ1w zG9yIfP-fBjneLH2VJds9Q+ZDNs~{4n87(F%LZ!aSAgEQD4b_y(GW7wXD3px=midIl)H^H0FpU;iYS?vA%m&zIihO#UE zmJb%>?b^D@%FLDG`FSw0gNyq<-W3(em!d%QP9bwZ`{rP!hP4V~J@h6?KThgL!)&tloFIuPW(-c{ z0)?4(cuDlkMg|TqzdpMVyk&O)+ zt)od}o}zVvG$&0f@w!%!Q?G~!7}k7S>&&po?&tBi+Cv=R&*;;re>L zss;ka_8itF88LMb9!G)_B7#_(oxEyK1y>JRF!*8>^S#u6p96LUBHj;8Q_gUjOY>?1%-8>%5P!57g?EQd0 z-5=#@2?-ydom2n%5&C$s|1{bc$H~ZWs|}9OdUaJmpDz=k7|r$G6{mVNd-*f@CDg!2 zAZ`AUqnnP&cbIrS1$)hV(Nv-o1@o#qHaG4J`CScP+T`xgh_(qoAelxFRSRdjdv_p zY1umQDorg31XzwNN=xdB1?TuG2hey7H#2>$;|2AFimF)iy4mnm2=q;zt4U)CLfMU1 zhL5T61pfu?qT%##FQp*W55>k{Z8auFe3h7U#YmE{dH69GJD=&&S(VPo!vzZyTU$BU&Y+$5{%iE5kNRYaAn`WRqE zQ)H3hU-m6nnDpxN;H*@5%K$-fe`ie_nA+nSh8v__z(CF#XW&KVthE(ox+-LU^7NeJ z;gomBRQ)^p8ECoXqJ89^$nV#>yuawFjIa&u^3u~?v(UYb@OR?}Apbzc!VK#C z97GeR{`$rnX|;L1yb6r$Y_AdFhvVqid` zC%7v{0ydc zgSEJWIGw9`f<;_qGRCU+9EmLZzz+q#Cf&2mwKg=}*VD`IelpjM+;yMjgJXiPk3Mu) zS!PS+i7BK=v*9~QGR|b^;-6?Gn(b%3+=?vy={cu>)ohC*AKnv>%@}%?=VlziK5rvE zejfdPev;NsL$(uwP-4Btjod1-*`N6bxw0DYV83SkgGeQ8J464-snt)GKw{~^<||3I zevwU#WIORljWzU3g_fAz$eG=FC#?YM*XohZgJus=YUeaFh1U0Gv4IXmq!oB`#HPz8 zl}_<&@=IACHm6}pLIm_W^<4LE>IlcBz9hhYStF1fk@~cs?6~y4)X7PoC-+0$CbSqf zZ=`7vi*%u(Ja(K6r_<7)tH0kOJP9YieOn5Lb|2{m!>b-2pkG4@L~yb4tEs_pF2>4B z{BKMfu&916!8D;Bd?~?rkR6w)$_oA9rBhv`#Qr6?#o$;EEB3r+iWheJ7W;7?a=URZ zI)%aDutf|vEKwPhJ8??!1~k11$`mY;IS)j;SB0)-eN0o1*Oam$J08pn(s4NK3M}BM z=3wV?ri`;S8Cei%hU5sW%@M>EgkmgR0>7iUGecnniSEXdQ8lN}c6_F3J$hzoI0IMW6%_{{!T z?$wPnOx|XZ0A?QP1d14Bl0p37bZ)A5CiNZBQR#zDPo-*cw`a*+;>|eyu7EZ3z8Cjb zTz(vHID!>oY@ltRgX);dLQis*{`iA~-)F~m2d7%3(KqzUB(s1J@ zrc;TGB7-CbOu+(3qn}MyGnBArrv~M=b_!=+dFKgQB?oSGo%U!z@cG=z85{`ZbiLO% zvloeK(hdK)&_*g^5Ov%@-pe@(C2~iepUia&mG(=aFe9QDg$z`*Eal_-p{dlMJBP&6 zxtMyWwxP$CicOBil!1dL$Wd*0s(?sIn^K%Q=T_~jfR835ss?13=oP>pZ1}@aW{&#R z2H56u3^=`QtPCGPn^*=6wrh7swL}2x2 z8TL!Bl2fI$1QaaGZNu*-KdLXx@;!d~Ubw*ZRn+X1F22*N9Go!`0az0(%+vHP= z&Cisqqp!M_Nk2X=5&{ST>9cw57L|Z>j^611loE&*HpV*+JGP0$ zGY{J%2anlj5U?;A(peb>S~x_OuqcVmPQmcX%7}&<`CBZ8%JAb1P0z4MmT|wgiA9ZF zDH1wq7A1#wi;JoZoV4h|YdVO%a^QjR+aH8P5R8YiaZSC_ayGkEG>#1iVy%&kjKncV zc{-um>qLQ+X<7FW;7bz^21PZQNMlaYOwL$?oLC-u49S_P4 zpC)lQzPFWOOzCY?(v{i#;dZmGiBTv=ZXKQPMr8e}k8R*8X?$0wo00S?{Q*5~C@r3qn* zo#_YU55&s=8cMw2#+rw-N&XyVu7cYmq@(s($?oIPRPDS5A}iyr%^+k`4n*@_1ltHi z%p+rtBcFi&j5s=YP&XpsZ%H}WlK+5MqdXLm348u|N< z$(5SBwHq#&l5{0z6my!YrIo@3pm2d=QT_qq%lv*!?DroE-MdmmdyX$UHN{Kv+=79D zgUQ|y0}M3QYxhRMwqfTngB!1J*0Vyv##@zgO!Kw*6wEtBs0pnVB7HPc4^E`_7mWcz zfyGkMmN~EW_mmEIwADZ6@?~Z34S${E=)Y>SH@ts)Qcs5#>5im@2`f`h7YxJ zkXMq0dC51tI8t6qrag4ZJ`@}xC3qLc3xr-Niz3MYynv{&RiRj z#DshiwjF~2e`p3}KBbne!Doie%(6YmOxX^2aNSuA{tCIJeN{twaYAS2RQ~X!9d&qv za-PR;QJ)$Mj|%WbP7>(EwhXFkAWYe|5b7d{cN8Upr+V<;9rBWlC=4w*wQY8J;*jBS z$gX4{pPvKLu^)6l>km^_p}0jd$LonHh;z1^3=U_>E~4S*=)KKC|bFE(3TSkGuCI``Xz%-f?5CALS?P%r2C6ozBDqFq>bPQjn8SmU0M zrSJ{sq1w6EhUm(u{p-;&@>=t*#{FGEV$gT$7h-TMMVW@%Ab*zAS!r=4wu{s9Mwz2mS*B%un^v9E|oPzv6zyEwO zl~%^_??8Nt+D?CZp*$^Wh)YV+OUp|rJO!Hm`$R2PSfQ5*A?V}|Be`kPrK^In=$NT* zR_?2B?Ocns*rsY2{N?IqwS^Fp($>v_4L(?|6z=B`+T@d}iW7Anik) z6c!DE#4eL;*+cYGFzzsiT$c>G|C^DoyfRQ$&N7rNeS9O@w99TXnMsKU%n=yjid9P{ z7g?&XwRnXW#39>ZKv6>VSX+Y#_A4||zhqMq;*f?fpbzkuTwmC60VtBo56S}XTN@|K zsnolX&J)&lUc>DeUGvLq20R>5c|=ZdBu@_r+;ZeN9-wgy9+jjL&D|F;+&+5T4xB*2 zQJLN|i4(h0xInvYq(MzaY**$MoYHB(iJrh!MsPMQ9j-}l#Jh}?G?-NL79 zmhyD%epR3U&$UxgP*WBaS7C5+b8;SggVaNW?7QxLqU-@RrEU}%Q`bTPA%R!uD-?Bg zS+kc1>-6vroK?5gJenRyY@ImrseEEd`UmD=jl{fuq5V9R*cvf!&YA5h+MD7d?xY`H zhX)Om>$uA3KwW@o&@PwLf1zJ^Ili}bALKUn} zpVPCV*(0(93EBy$8MlSV#s~`dK__^kjj{}xJie&RngBJ&^0dqa^3X`kfwS%#}v$?}lVC7S*5!zAu zk}VqfF^Kr)Fmg!lJ96$+X&sy*A%Pb(nHveNr{%V%tp@&`LjE6`ldYY;m9dMl)i0c* z3+&LBFVs6Pwfx?bKW=Uhb31xR0~6zasRso^2mkNe89e>)zm5XGUvGceU_kNBGvL=< z@p4ase>VpKDe&|9JHY67!rw~wKNH?R{U`rO_+N$mf5-e1^0MUqGlt@6<^I27{*d#x z^7GFe!>8NiPtHHZ=wE7iS?%|k!TI#}{+;oMkiV74dxl6p1+e`I`G3^~Ug~&RlHeKf zn&iJA{*dvv#G7Y^8{Rq`XO_uJ zX)n_Yo>5+`{~y#Jvi=rx{miQ6c**)tU+n(`{UPaZONh@TFzLULe)Dg=)b!Ff_DphA{=Z3oi25~vyqp_6qdIhd z8`%GuDZEtj(mXylzt{V%`ESPZCE@q+`xi=Gj>za##zOfMO~_vn{J`sXO{O9%cY;~)E^_y+PRJoi^{E;dN!)7j|Lmq0-N E2dWrN9{>OV literal 0 HcmV?d00001 diff --git a/tests/tools/test_initializers.py b/tests/tools/test_initializers.py new file mode 100644 index 0000000..d4dab18 --- /dev/null +++ b/tests/tools/test_initializers.py @@ -0,0 +1,24 @@ +from __future__ import division, print_function, absolute_import + +from CDTools.tools import initializers +import numpy as np +import torch as t + + + +def test_gaussian(): + # Generate gaussian as a numpy array (square array) + shape = [10, 10] + sigma = [2.5, 2.5] + center = ((shape[0]-1)/2, (shape[1]-1)/2) + y, x = np.mgrid[:shape[0], :shape[1]] + np_result = 10*np.exp(-((x-center[1])/sigma[1])**2-((y-center[0])/sigma[0])**2) + assert(np.allclose(initializers.gaussian([10, 10], 10, [2.5, 2.5]), np_result)) + + # Generate gaussian as a numpy array (rectangular array) + shape = [10, 5] + sigma = [2.5, 2.5] + center = ((shape[0]-1)/2, (shape[1]-1)/2) + y, x = np.mgrid[:shape[0], :shape[1]] + np_result = 10*np.exp(-((x-center[1])/sigma[1])**2-((y-center[0])/sigma[0])**2) + assert(np.allclose(initializers.gaussian([10, 5], 10, [2.5, 2.5]), np_result)) diff --git a/tests/tools/test_projectors.py b/tests/tools/test_projectors.py new file mode 100644 index 0000000..36ce13b --- /dev/null +++ b/tests/tools/test_projectors.py @@ -0,0 +1,21 @@ +from __future__ import division, print_function, absolute_import + +from CDTools.tools import cmath +from CDTools.tools import projectors +import numpy as np +import torch as t +from scipy.fftpack import fftshift, ifftshift + +def test_modulus(): + # Create a complex array with modulus 12 and phase pi/4 + np_result = 6**.5*(np.ones((10,10))+1j*np.ones((10,10))) + + assert(np.allclose(cmath.torch_to_complex(projectors.modulus(t.ones((10,10,2)), 12*t.ones((10,10)))), np_result)) + + +def test_support(): + # Test masking + support = t.zeros((10,10)) + np_result = np.zeros((10,10)) + + assert(np.allclose(cmath.torch_to_complex(projectors.support(t.ones((10,10,2)), support)), np_result)) diff --git a/tests/tools/test_propagators.py b/tests/tools/test_propagators.py new file mode 100644 index 0000000..b038d37 --- /dev/null +++ b/tests/tools/test_propagators.py @@ -0,0 +1,48 @@ +from __future__ import division, print_function, absolute_import + +from CDTools.tools import cmath +from CDTools.tools import initializers +from CDTools.tools import propagators + +import numpy as np +import torch as t +import pytest +import scipy.misc +from scipy.fftpack import fftshift, ifftshift + +@pytest.fixture(scope='module') +def exit_waves_1(): + # Import scipy test image and add a random phase + object = scipy.misc.ascent()[0:64,0:64].astype(np.complex128) + arr = np.random.random_sample((64,64)) + object *= (arr+(1-arr**2)**.5*1j) + + # Construct wavefront from image + probe = initializers.gaussian([64, 64], 1e3, [5, 5])*(1+1j) + return cmath.complex_to_torch(probe*object) + + +def test_far_field(exit_waves_1): + # Far field diffraction patterns calculated by numpy with zero frequency in center + np_result = np.fft.fftshift(np.fft.fft2(cmath.torch_to_complex(exit_waves_1))) + + assert(np.allclose(np_result, cmath.torch_to_complex(propagators.far_field(exit_waves_1)))) + + +def test_inverse_far_field(exit_waves_1): + # We want the inverse far field to map back to the exit waves with no intensity corrections + np_result = exit_waves_1 + # Far field result for exit waves calculated with numpy + far_field_np_result = cmath.complex_to_torch(np.fft.fftshift(np.fft.fft2(cmath.torch_to_complex(exit_waves_1)))) + + assert(np.allclose(np_result, propagators.inverse_far_field(far_field_np_result))) + + + +def test_near_field(exit_waves_1): + pass + + + +def test_get_exit_waves(exit_waves_1): + pass