Difference map implementation and switchable interpolation

[hruth]

Updates

• Added difference map reconstructor
• Added new objects option called patch_interpolation_method that selects the type of interpolation used to extract object patches in PlanarObject.extract_patches. The default value is fourier. Other options are bilinear and nearest. This setting also determines what kind of interpolation is used in the update_object method of DMReconstructor.

Notes/Issues

• Difference map reconstruction does not work with fourier interpolation. I did some tests trying different interpolators for the patch extraction and placement, and reconstructions where fourier interpolation was used always blew up after a while. Here is a plot showing that:

• The difference map error (defined as $||\Psi_{n+1} - \Psi_n||$) does not converge, instead it oscillates. Here is an example difference map error:

I have more detailed notes on the performance in a powerpoint we can go over sometime.

[mdw771]

Not sure if this should be in ForwardModel. Maybe we can put it in AnalyticalIterativePtychographyReconstructor.

[hruth]

Ok, it is now a static method of AnalyticalIterativePtychographyReconstructor

[mdw771]

You can use the existing routine to calculate the object preconditioner, which uses less GPU memory. Add this to DMReconstructor:

    def run_pre_epoch_hooks(self) -> None:
        self.update_preconditioners()

then you can access the preconditioner here as self.parameter_group.object.preconditioner.

[hruth]

I just added this, but self.update_preconditioners() had to be put directly before the object update in order for it to work properly

[mdw771]

1e-10 is too small. The regions not illuminated by the probe can have zero values, and they can make the algorithm unstable. In LSQML this is what we do:

        preconditioner = self.parameter_group.object.preconditioner
        delta_o_hat = delta_o_hat / torch.sqrt(
            preconditioner**2 + (preconditioner.max() * alpha_mix) ** 2
        )

alpha_mix is 0.05 by default and it is used in pty-chi and Tike, but PtychoShelves doesn't use it (i.e., their alpha_mix is 1).

[hruth]

I tried using this:

        updated_object = object_numerator / torch.sqrt(
            object_denominator**2 + (0.05 * object_denominator.max()) ** 2
        )

and it did prevent the instability when using Fourier interpolation, but the difference map error still oscillates.

The 1e-10 value does work fine for bilinear and nearest interpolation, I think it is because they are more precise than place_patches_fourier_shift.

Anyway, I'll use the new code snippet since it is working with all three of the interpolators.

[mdw771]

The reason why Fourier shift is unstable with 1e-10 is likely that it produces small non-zero values where it should be exactly 0. Bilinear and especially nearest are not more precise in well illuminated areas but they preserve the zeros.

[hruth]

Yes, sorry, I wasn't being precise in my language. That's what I meant.

[mdw771]

The object update denominator might have contributed to the instability you saw. Could you test it again after changing it to sqrt(preconditioner ** 2 + preconditioner.max() ** 2)?

[stevehenke]

Nice work, Hanna!

[stevehenke]

Can this conditional be removed in favor of always doing the imports? The quoting business on type hints can be avoided by adding: "from __ future __ import annotations" to the top of the file.

[mdw771]

I ran into circular import issues in quite a few files without using TYPE_CHECKING. With this conditional the intention is to keep the import as unidirectional as possible (API imports from backend but not the other way round).

[stevehenke]

How about addressing this by creating an issue to (1) identify where the circular import issues are, and (2) refactor to fix the circular imports?

[mdw771]

I looked into this sometime in the past. Most of the circular imports are caused by importing ptychi.api in other modules such as reconstructors, forward model, etc. The API imports are mainly used for (1) type annotation and (2) getting the enums in ptychi.api.enums to use them in conditionals (e.g., if options.noise_model == api.NoiseModels.GAUSSIAN).

(1) can't be helped with quoting or __future__.annotations as long as the imports are run in the global namespace. We could alternatively move the imports inside classes, but not sure that's a good way.

(2) could probably be fixed by moving enums out of api/, but users of the Python API will have to import the enums from somewhere else which I feel isn't intuitive.

I used this TYPE_CHECKING thing because it just prevents the import of api in the backend at runtime, so at runtime we have a clean one-way dependency. I think this also reduces the chance of circular imports for future developments.

[stevehenke]

Would it be beneficial to use an Enum instead of a Literal for function_type?

[mdw771]

function_type is internal and not intended to be used by users, so I think this is fine.

[stevehenke]

This is probably a bigger discussion, but when base options are added what should be done to ensure that they are supported by all of the reconstruction algorithms?

[mdw771]

Some of the base options are used in the abstract classes of the reconstructors (like IterativePtychographyReconstructor). For this patch_interpolation_method, looks like it sets the function used in PlanarObject.extrract/place_patches. The PlanarObject class is used in all reconstructors, so it will apply to all.

[hruth]

So this actually is an issue. I added the option to switch interpolators because of some issues I was having with the Fourier interpolation method, but that is resolved now. You can use bilinear, nearest, or Fourier interpolation with difference map... but it doesn't work with the ePIE reconstructor (and probably won't with the LSQML reconstructor) without replacing all calls to place_patches_fourier_shift with a call to the output of the get_patch_interp_function_by_enum function. Mixing and matching interpolators doesn't give good results. So I will have to update this.

This new feature won't actually change the default behavior of any of the algorithms since everything is defaulting to using the fourier interpolation anyway. But, it won't be hard to finish this update in case we want to change the interpolators in the future.

[mdw771]

It looks good to me now and we can merge it after the conflicts are resolved. Well done!