From 59d103f6d61c3338ccea689d30e52cfc70cb5627 Mon Sep 17 00:00:00 2001 From: Abe Levitan Date: Mon, 1 Jul 2024 17:16:07 +0200 Subject: [PATCH] update RPI to work with the latest data saving method, change dependence to latest pytorch for LBFGS to work on complex leaf tensors, and fix an issue with the test for the FRCs --- setup.py | 4 +- src/cdtools/models/__init__.py | 1 + src/cdtools/models/fancy_ptycho.py | 2 +- src/cdtools/models/rpi.py | 492 ++++++++++++++++++----------- tests/tools/test_analysis.py | 4 +- 5 files changed, 315 insertions(+), 188 deletions(-) diff --git a/setup.py b/setup.py index 8868f61..59d4600 100644 --- a/setup.py +++ b/setup.py @@ -6,7 +6,7 @@ with open("README.md", "r") as fh: setuptools.setup( name="cdtools", version="0.2.0", - python_requires='>3.7', # recommended minimum version for pytorch + python_requires='>3.8', # recommended minimum version for pytorch 2.3.0 author="Abe Levitan", author_email="abraham.levitan@psi.ch", description="Tools for coherent diffractive imaging and ptychography", @@ -18,7 +18,7 @@ setuptools.setup( "scipy>=1.0", "matplotlib>=2.0", # 2.0 has better colormaps which are used by default "python-dateutil", - "torch>=1.9.0", #1.9.0 supports autograd on indexed complex tensors + "torch>=2.3.0", #2.3.0 is the earliest release for which L-BFGS works directly on complex-valued leaf tensors "h5py>=2.1"], extras_require={ 'tests': [ diff --git a/src/cdtools/models/__init__.py b/src/cdtools/models/__init__.py index 9a10e91..85a190b 100644 --- a/src/cdtools/models/__init__.py +++ b/src/cdtools/models/__init__.py @@ -33,6 +33,7 @@ from cdtools.models.bragg_2d_ptycho import Bragg2DPtycho from cdtools.models.multislice_2d_ptycho import Multislice2DPtycho from cdtools.models.multislice_ptycho import MultislicePtycho from cdtools.models.rpi import RPI +from cdtools.models.bragg_rpi import BraggRPI from cdtools.models.multimode_rpi import MultimodeRPI from cdtools.models.time_resolved_ptycho_calibration import TimeResolvedPtychoCalibration from cdtools.models.time_resolved_rpi import TimeResolvedRPI diff --git a/src/cdtools/models/fancy_ptycho.py b/src/cdtools/models/fancy_ptycho.py index 252eda8..4249c05 100644 --- a/src/cdtools/models/fancy_ptycho.py +++ b/src/cdtools/models/fancy_ptycho.py @@ -79,7 +79,7 @@ class FancyPtycho(CDIModel): self.register_buffer('mask', t.tensor(mask, dtype=t.bool)) - probe_guess = t.tensor(probe_guess, dtype=t.complex64) + probe_guess = t.tensor(probe_guess, dtype=t.complex4) obj_guess = t.tensor(obj_guess, dtype=t.complex64) # We rescale the probe here so it learns at the same rate as the diff --git a/src/cdtools/models/rpi.py b/src/cdtools/models/rpi.py index 3ebe77b..d8ae029 100644 --- a/src/cdtools/models/rpi.py +++ b/src/cdtools/models/rpi.py @@ -12,32 +12,24 @@ import time __all__ = ['RPI'] # -# This model has to work a bit differently from a ptychography model +# This model has works a bit differently from a ptychography model # because a typical RPI dataset will have lots of images, each of which -# can be reconstructed on it's own. I can see a few ways to approach this: +# can be reconstructed on it's own. # -# 1) Enfore 1 image per dataset as a restriction for internal consistency -# of the "model" idea +# I made the decision that the model itself will only store one object guess, +# so it can't (for example) simultaneously reconstruct images from all the +# diffraction patterns within a large dataset. This could be possible to +# implement in the future, but is not implemented at present. # -# 2) Override some of the base functions of the CDIModel to accept optional -# parameters that make it work on larger datasets +# To support the approach I chose, I have added a "subset" option to all of the +# optimization functions. To do an RPI reconstruction, the "subset" option +# should be set to a container (list, set, etc) which includes only the index +# of the diffraction pattern to be reconstructed. A few additional decisions: # -# I think #1 is basically untenable, because it would require creating -# a new dataset and model for each reconstruction - which means instantiating -# tons of stuff and storing all sorts of excess information for each frame, -# when that could easily be reused for other frames. As a result, I think I -# will try the following pattern: -# -# 1) A model will contain a single object "guess" at all times -# 2) Constructing a model from a data will automatically instantiate the object guess from the first diffraction pattern in the dataset -# 3) All the optimization functions will get an additional argument for the index of the diffraction pattern to reconstruct. This could either be handled by editing the CDIModel base class to pass through some kwargs, or by overriding all the optimization functions explicitly. -# 4) A few convenience functions can be written to re-initialize the object array from any image / index in the dataset +# 1) Constructing a model from a dataser will automatically instantiate +# the object guess from the first diffraction pattern in the dataset # 5) A new function can be written to reconstruct the entire dataset by running through each pattern one at a time. # -# The advantage of this approach is that I can start by writing a class that -# will only reconstruct the first diffraction pattern from a dataset and then -# extend it. -# # Final note: It is worth seeing whether it is possible to include the # probe propagation explicitly as a parameter which can be reconstructed # via gradient descent. @@ -47,93 +39,137 @@ __all__ = ['RPI'] __all__ = ['RPI'] class RPI(CDIModel): - - @property - def obj(self): - return t.complex(self.obj_real, self.obj_imag) - - @property - def weights(self): - ws = t.complex(self.weights_real, self.weights_imag) - return ws / 10# / self.obj_real.size().numel() - - - def __init__(self, wavelength, detector_geometry, probe_basis, - probe, obj_guess, - background=None, mask=None, saturation=None, - obj_support=None, oversampling=1, weight_matrix=False, - simulate_finite_pixels=False): - + def __init__( + self, + wavelength, + detector_geometry, + probe_basis, + probe, + obj_guess, + background=None, + mask=None, + saturation=None, + obj_support=None, + oversampling=1, + weight_matrix=False, + propagation_distance=0, + units='um', + dtype=t.float32, + ): + super(RPI, self).__init__() - self.wavelength = t.tensor(wavelength) - self.detector_geometry = copy(detector_geometry) + complex_dtype = (t.ones([1], dtype=dtype) + + 1j * t.ones([1], dtype=dtype)).dtype - det_geo = self.detector_geometry - if hasattr(det_geo, 'distance'): - det_geo['distance'] = t.tensor(det_geo['distance']) - if hasattr(det_geo, 'basis'): - det_geo['basis'] = t.tensor(det_geo['basis']) - if hasattr(det_geo, 'corner'): - det_geo['corner'] = t.tensor(det_geo['corner']) + self.register_buffer('wavelength', + t.tensor(wavelength, dtype=dtype)) + self.store_detector_geometry(detector_geometry, + dtype=dtype) - self.probe_basis = t.tensor(probe_basis) + # NOTE: It is required that the probe basis match the exit wave + # basis. If, for example, the probe reconstruction from ptychography + # used a bandlimiting constraint and had a larger basis, the user is + # expected to upsample it explicitly before doing RPI. + self.register_buffer('probe_basis', + t.tensor(probe_basis, dtype=dtype)) scale_factor = t.tensor([probe.shape[-1]/obj_guess.shape[-1], probe.shape[-2]/obj_guess.shape[-2]]) - self.obj_basis = self.probe_basis * scale_factor + self.register_buffer('obj_basis', + (self.probe_basis * scale_factor).to(dtype=dtype)) - # Maybe something to include in a bit - # self.surface_normal = t.tensor(surface_normal) - - self.saturation = saturation - - if mask is None: - self.mask = mask + if saturation is None: + self.saturation = None else: - self.mask = t.tensor(mask, dtype=t.bool) + self.register_buffer('saturation', + t.tensor(saturation, dtype=dtype)) - - self.probe = t.tensor(probe, dtype=t.complex64) - + # not sure how to make this a buffer, or if I have to + self.units = units + + if mask is None: + self.mask = None + else: + self.register_buffer('mask', + t.tensor(mask, dtype=t.bool)) + + self.register_buffer('probe', t.tensor(probe, dtype=complex_dtype)) + + # We always use multi-modes to store the object, so we convert it + # if we just get a single 2D array as an input if obj_guess.dim() == 2: obj_guess = obj_guess[None, :, :] - - obj_guess = t.tensor(obj_guess, dtype=t.complex64) - self.obj_real = t.nn.Parameter(obj_guess.real) - self.obj_imag = t.nn.Parameter(obj_guess.imag) + self.obj = t.nn.Parameter(t.tensor(obj_guess, dtype=complex_dtype)) - self.weights_real = t.nn.Parameter(t.eye(probe.shape[0])* 10)# * self.obj_real.size().numel()) - self.weights_imag = t.nn.Parameter(t.zeros(probe.shape[0])) + self.weights = t.nn.Parameter( + t.eye(probe.shape[0], dtype=complex_dtype)* 10) if not weight_matrix: - self.weights_real.requires_grad=False - self.weights_imag.requires_grad=False + self.weights.requires_grad=False - # Wait for LBFGS to be updated for complex-valued parameters - # self.obj = t.nn.Parameter(obj_guess.to(t.float32)) - if background is None: background = 1e-6 * t.ones(self.probe[0].shape, - dtype=t.float32) - - self.background = t.tensor(background, dtype=t.float32) + dtype=dtype) - if obj_support is not None: - self.obj_support = obj_support - self.obj.data = self.obj * obj_support[None, ...] - else: - self.obj_support = t.ones_like(self.obj[0, ...]) + self.register_buffer('background', + t.tensor(background, dtype=t.float32)) - self.oversampling = oversampling + if obj_support is None: + obj_support = t.ones_like(self.obj[0, ...], dtype=int) - self.simulate_finite_pixels=simulate_finite_pixels + self.register_buffer('obj_support', + t.tensor(obj_support, dtype=int)) + + self.obj.data = self.obj * self.obj_support[None, ...] + + self.register_buffer('oversampling', + t.tensor(oversampling, dtype=int)) + + self.register_buffer('propagation_distance', + t.tensor(propagation_distance, dtype=dtype)) + + # The propagation direction of the probe. For now it's fixed, + # but perhaps it would need to be updated in the future + self.register_buffer('prop_dir', + t.tensor([0, 0, 1], dtype=dtype)) + + # This propagator should be able to be multiplied by the propagation + # distance each time to get a propagator + #universal_propagator = t.angle(ggasp( + # self.probe.shape[-2:], + # self.probe_basis, self.wavelength, + # t.tensor([0, 0, self.wavelength/(2*np.pi)], dtype=dtype), + # propagation_vector=self.prop_dir, + # dtype=complex_dtype, + # propagate_along_offset=True)) + + # TODO: probably doesn't support non-float-32 dtypes + #self.register_buffer('universal_propagator', + # universal_propagator) @classmethod - def from_dataset(cls, dataset, probe, obj_size=None, background=None, mask=None, n_modes=1, saturation=None, scattering_mode=None, oversampling=1, initialization='random', weight_matrix=False, probe_threshold=0, simulate_finite_pixels=False): + def from_dataset( + cls, + dataset, + probe, + obj_size=None, + background=None, + mask=None, + n_modes=1, + saturation=None, + scattering_mode=None, + oversampling=1, + initialization='random', + weight_matrix=False, + probe_threshold=0, + dtype=t.float32, + ): + complex_dtype = (t.ones([1], dtype=dtype) + + 1j * t.ones([1], dtype=dtype)).dtype wavelength = dataset.wavelength det_basis = dataset.detector_geometry['basis'] @@ -143,10 +179,10 @@ class RPI(CDIModel): # always do this on the cpu get_as_args = dataset.get_as_args dataset.get_as(device='cpu') + # We only need the patterns here, not the inputs associated with them. _, patterns = dataset[:] dataset.get_as(*get_as_args[0],**get_as_args[1]) - # Then, generate the probe geometry from the dataset ewg = tools.initializers.exit_wave_geometry @@ -156,8 +192,7 @@ class RPI(CDIModel): distance, oversampling=oversampling) - if not isinstance(probe,t.Tensor): - probe = t.as_tensor(probe) + probe = t.as_tensor(probe) # Potentially need all of this orientation stuff later @@ -192,37 +227,16 @@ class RPI(CDIModel): and dataset.mask is not None: mask = dataset.mask.to(t.bool) - # Now we initialize the object + # This will be superceded later by a call to init_obj, but it sets + # the shape if obj_size is None: - # This is a standard size for a well-matched probe and detector - obj_size = (np.array(probe_shape) // 2).astype(int) + obj_size = (np.array(self.probe.shape[-2:]) // 2).astype(int) - if initialization.lower().strip() == 'random': - # I think something to do with the fact that the object is defined - # on a coarser grid needs to be accounted for here that is not - # accounted for yet - scale = t.sum(patterns[0]) / t.sum(t.abs(probe)**2) - obj_guess = scale * t.exp(2j * np.pi * t.rand([n_modes,]+list(obj_size))) - elif initialization.lower().strip() == 'uniform': - # I think something to do with the fact that the object is defined - # on a coarser grid needs to be accounted for here that is not - # accounted for yet - scale = t.sum(patterns[0]) / t.sum(t.abs(probe)**2) - obj_guess = scale * t.ones([n_modes,]+list(obj_size), dtype=probe.dtype) - elif initialization.lower().strip() == 'spectral': - if background is not None: - obj_guess = initializers.RPI_spectral_init( - patterns[0], probe, obj_size, mask=mask, - background=background**2, n_modes=n_modes) - else: - obj_guess = initializers.RPI_spectral_init( - patterns[0], probe, obj_size, mask=mask, - n_modes=n_modes) - - else: - raise KeyError('Initialization "' + str(initialization) + \ - '" invalid - use "spectral", "uniform", or "random"') + dummy_init_obj = t.ones([n_modes, obj_size[0], obj_size[1]], + dtype=complex_dtype) + # This defines an object support in real space using the probe + # intensities, if requested probe_intensity = t.sqrt(t.sum(t.abs(probe)**2,axis=0)) probe_fft = tools.propagators.far_field(probe_intensity) pad0l = (probe.shape[-2] - obj_size[-2])//2 @@ -235,60 +249,187 @@ class RPI(CDIModel): obj_support = probe_lr > t.max(probe_lr) * probe_threshold obj_support = t.as_tensor(binary_dilation(obj_support)) - return cls(wavelength, det_geo, ew_basis, - probe, obj_guess, - background=background, mask=mask, saturation=saturation, - obj_support=obj_support, oversampling=oversampling, - weight_matrix=weight_matrix, - simulate_finite_pixels=simulate_finite_pixels) + rpi_object = cls(wavelength, det_geo, ew_basis, + probe, dummy_init_obj, + background=background, mask=mask, + saturation=saturation, + obj_support=obj_support, oversampling=oversampling, + weight_matrix=weight_matrix) + + # I don't love this pattern, where I do the "real" obj initialization + # after creating the rpi object. But, I chose this so that I could + # have a function to reinitialize the obj, for repeat reconstructions + # using the same rpi object, without repeating code. There probably + # is a better pattern for doing this. + rpi_object.init_obj(initialization, + pattern=dataset.patterns[0]) + + return rpi_object @classmethod - def from_calibration(cls, calibration, obj_size=None, n_modes=1, saturation=None): + def from_calibration( + cls, + calibration, + obj_size=None, + n_modes=1, + saturation=None, # TODO can we get this from the calibration? + initialization='random', + dtype=t.float32 + ): - wavelength = calibration['wavelength'] - probe_basis = t.as_tensor(calibration['basis']) - probe = t.as_tensor(calibration['probe']) + complex_dtype = (t.ones([1], dtype=dtype) + + 1j * t.ones([1], dtype=dtype)).dtype + + + wavelength = t.as_tensor(calibration['wavelength'], dtype=dtype) + probe_basis = t.as_tensor(calibration['obj_basis'], dtype=dtype) + # TODO this will fail if the probe from the calibration was restricted + # in Fourier space + probe = t.as_tensor(calibration['probe'], dtype=complex_dtype) if 'background' in calibration: - background = t.as_tensor(calibration['background']) + background = t.sqrt(t.as_tensor(calibration['background'], + dtype=dtype)) else: background = t.zeros_like(probe.real) if 'mask' in calibration: - mask = t.as_tensor(calibration['mask']) + mask = t.as_tensor(calibration['mask'], dtype=t.bool) else: mask = t.ones_like(probe.real, dtype=t.bool) - # Now we initialize the object + # This will be superceded later by a call to init_obj, but it sets + # the shape if obj_size is None: - # This is a standard size for a well-matched probe and detector - obj_size = (np.array(probe.shape[-2:]) // 2).astype(int) + obj_size = (np.array(self.probe.shape[-2:]) // 2).astype(int) - obj_guess = t.exp(2j * np.pi * t.rand([n_modes,]+obj_size)) + dummy_init_obj = t.ones([n_modes, obj_size[0], obj_size[1]], + dtype=complex_dtype) + # Pretty sure that this will not work for reflection-mode det_geo = {'distance': 1, - 'basis': wavelength / probe_basis} # Is this even right? + 'basis': wavelength / probe_basis} - return cls(wavelength, det_geo, probe_basis, - probe, obj_guess, - background=background, mask=mask) + rpi_object = cls( + wavelength, + det_geo, + probe_basis, + probe, + dummy_init_obj, + background=background, + mask=mask, + ) + + rpi_object.init_obj(initialization) + + return rpi_object - def random_init(self, pattern): - scale = t.sum(pattern) / t.sum(t.abs(self.probe)**2) - self.obj.data = scale * t.exp( - 2j * np.pi * t.rand(self.obj.shape)).to( + def init_obj( + self, + initialization_type, + obj_shape=None, + n_modes=None, + pattern=None, + ): + + # I think something to do with the fact that the object is defined + # on a coarser grid needs to be accounted for here that is not + # accounted for yet + + if initialization_type.lower().strip() == 'random': + self.random_init(obj_shape=obj_shape, n_modes=n_modes) + elif initialization_type.lower().strip() == 'uniform': + self.uniform_init(obj_shape=obj_shape, n_modes=n_modes) + elif initialization_type.lower().strip() == 'spectral': + if pattern is None: + raise KeyError( + 'A pattern must be supplied for spectral initialization') + + self.spectral_init(pattern=pattern, + obj_shape=obj_shape, + n_modes=n_modes) + else: + raise KeyError('Initialization "' + str(initialization) + \ + '" invalid - use "spectral", "uniform", or "random"') + + + def get_obj_shape_and_n_modes(self, obj_shape=None, n_modes=None): + """Sets defaults for obj shape and n modes""" + + if obj_shape == None: + if hasattr(self, 'obj'): + obj_shape = self.obj.shape[-2:] + else: + obj_size = (np.array(self.probe.shape[-2:]) // 2).astype(int) + obj_shape = [obj_size, obj_size] + + if n_modes == None: + if hasattr(self, 'obj'): + n_modes = self.obj.shape[0] + else: + n_modes = 1 + + return n_modes, obj_shape + + + def uniform_init(self, obj_shape=None, n_modes=None): + """Sets a uniform object initialization""" + + n_modes, obj_shape = self.get_obj_shape_and_n_modes( + obj_shape=obj_shape, n_modes=n_modes) + + obj_guess = t.ones( + [n_modes,]+list(obj_shape), + dtype=self.probe.dtype, + device=self.probe.device, + ) + + if hasattr(self, 'obj'): + self.obj.data = obj_guess + else: + self.obj = t.nn.Parameter(obj_guess) + + + def random_init(self, obj_shape=None, n_modes=None): + """Sets a uniform amplitude object initialization with random phase""" + + n_modes, obj_shape = self.get_obj_shape_and_n_modes( + obj_shape=obj_shape, n_modes=n_modes) + + obj_guess = t.exp( + 2j * np.pi * t.rand([n_modes,]+list(obj_shape))).to( + dtype=self.probe.dtype, device=self.probe.device) + + if hasattr(self, 'obj'): + self.obj.data = obj_guess + else: + self.obj = t.nn.Parameter(obj_guess) + + + def spectral_init(self, pattern, obj_shape=None, n_modes=None): + """Initializes the object with a spectral method""" + + n_modes, obj_shape = self.get_obj_shape_and_n_modes( + obj_shape=obj_shape, n_modes=n_modes) + + if self.background is not None: + background = self.background**2 + else: + background = None + + obj_guess = initializers.RPI_spectral_init( + pattern, + self.probe, + obj_shape, + mask=self.mask, + background=background, + n_modes=n_modes).to( dtype=self.obj.dtype, device=self.obj.device) - def spectral_init(self, pattern): - if self.background is not None: - self.obj.data = initializers.RPI_spectral_init( - pattern, self.probe, self.obj.shape[-3:-1], mask=self.mask, - background=self.background**2, n_modes=self.obj.shape[0]).to( - dtype=self.obj.dtype, device=self.obj.device) + if hasattr(self, 'obj'): + self.obj.data = obj_guess else: - self.obj.data = initializers.RPI_spectral_init( - pattern, self.probe, self.obj.shape[-3:-1], mask=self.mask, - n_modes=self.obj.shape[0]).to( - dtype=self.obj.dtype, device=self.obj.device) + self.obj = t.nn.Parameter(obj_guess) + # Needs work def interaction(self, index, *args): @@ -297,7 +438,6 @@ class RPI(CDIModel): # "input" parameters (such as translations for a ptychography dataset). # This makes it seamless to use such a dataset even though those # extra arguments will not be used. - all_exit_waves = [] @@ -348,13 +488,11 @@ class RPI(CDIModel): self.background, measurement=tools.measurements.incoherent_sum, saturation=self.saturation, - oversampling=self.oversampling, - simulate_finite_pixels=self.simulate_finite_pixels) + oversampling=self.oversampling) return m def loss(self, sim_data, real_data, mask=None): return tools.losses.amplitude_mse(real_data, sim_data, mask=mask) - #return tools.losses.poisson_nll(real_data, sim_data, mask=mask) def regularizer(self, factors): if self.obj.shape[0] == 1: @@ -363,29 +501,6 @@ class RPI(CDIModel): return factors[0] * t.sum(t.abs(self.obj[0,:,:])**2) \ + factors[1] * t.sum(t.abs(self.obj[1:,:,:])**2) - def to(self, *args, **kwargs): - super(RPI, self).to(*args, **kwargs) - self.wavelength = self.wavelength.to(*args,**kwargs) - # move the detector geometry too - det_geo = self.detector_geometry - if hasattr(det_geo, 'distance'): - det_geo['distance'] = det_geo['distance'].to(*args,**kwargs) - if hasattr(det_geo, 'basis'): - det_geo['basis'] = det_geo['basis'].to(*args,**kwargs) - if hasattr(det_geo, 'corner'): - det_geo['corner'] = det_geo['corner'].to(*args,**kwargs) - - if self.mask is not None: - self.mask = self.mask.to(*args, **kwargs) - - self.probe = self.probe.to(*args,**kwargs) - self.probe_basis = self.probe_basis.to(*args,**kwargs) - self.obj_basis = self.obj_basis.to(*args,**kwargs) - self.obj_support = self.obj_support.to(*args,**kwargs) - self.background = self.background.to(*args, **kwargs) - - # Maybe include in a bit - #self.surface_normal = self.surface_normal.to(*args, **kwargs) def sim_to_dataset(self, args_list): raise NotImplementedError('No sim to dataset yet, sorry!') @@ -412,22 +527,33 @@ class RPI(CDIModel): ] - def save_results(self, dataset=None, full_obj=False): + def save_results(self, dataset=None): # dataset is set as a kwarg here because it isn't needed, but the # common pattern is to pass a dataset. This makes it okay if one # continues to use that standard pattern + + # This will save out everything needed to recreate the object + # in the same state, but it's not the best formatted. For example, + # "background" stores the square root of the background, etc. + base_results = super().save_results() + probe_basis = self.probe_basis.detach().cpu().numpy() obj_basis = self.obj_basis.detach().cpu().numpy() probe = self.probe.detach().cpu().numpy() - # Provide the option to save out the subdominant objects or - # just the dominant one - if full_obj: - obj = self.obj.detach().cpu().numpy() - else: - obj = self.obj[0].detach().cpu().numpy() + + # We only save out the top object mode, which is the final result + # The full object is still saved in the state_dict that is + # included in the base_results + obj = self.obj[0].detach().cpu().numpy() background = self.background.detach().cpu().numpy()**2 - return {'probe_basis': probe_basis, 'obj_basis': obj_basis, - 'probe': probe,'obj': obj, - 'background': background} + results = { + 'probe_basis': probe_basis, + 'obj_basis': obj_basis, + 'probe': probe, + 'obj': obj, + 'background': background + } + + return {**base_results, **results} diff --git a/tests/tools/test_analysis.py b/tests/tools/test_analysis.py index e367483..2c7d1f8 100644 --- a/tests/tools/test_analysis.py +++ b/tests/tools/test_analysis.py @@ -272,13 +272,13 @@ def test_calc_frc(): n_pix, bins = np.histogram(Rs,bins=nbins) bins = bins[:-1] - frc = np.abs(numerator / np.sqrt(denominator_F1*denominator_F2)) + frc = numerator / np.sqrt(denominator_F1*denominator_F2) # This moves from combined-image SNR to single-image SNR snr /= 2 threshold = (snr + (2 * snr + 1) / np.sqrt(n_pix)) / \ (1 + snr + (2 * np.sqrt(snr)) / np.sqrt(n_pix)) - + test_bins, test_frc, test_threshold = analysis.calc_frc( obj1, obj2, basis, im_slice=np.s_[10:-10,20:-20], nbins=100, snr=2, limit='corner')