Skip to main content
`Ctrl`+`K`

# `pynx.ptycho`: 2D Ptychography

## Description

This modules allows the simulation and analysis of ptychography experiments, with the following features:

- 2D ptychography using a variety of algorithms: alternating projections, difference map (Thibault et al.), maximum likelihood conjugate gradient

- Works with any type of illumination (sharp or not)

- Object and/or Probe reconstruction

- Probe analysis (modes and focus)

- Incoherent background optimization

- GPU implementation using CUDA/OpenCL is available (and recommended), and is the main focus of current development:

  - example speed on single V1OO GPU card as of 2019/06: 13 ms/cycle for 1000 frames of 256x256 pixels and a simple alternating projection (17 ms/cycle for DM and 34 ms/cycle for ML)

  - GPU implementation allows using modes for probe and object

  - Maximum likelihood minimization (Poisson noise, regularisation)

- simple usage scripts are provided to analyse data from CXI files, ESRF beamlines (id01, id13, id16a), and ptypy files.

## API Reference

Note that the Python API is quickly evolving.

For regular data analysis, it is recommended to use the scripts which are stable, independently of the underlying Python API.

### 2D Ptychography (operator-based)

This is the Ptychography class, which can be used with operators

class pynx.ptycho.ptycho.OperatorPtycho(lazy=False)

  Base class for an operator on Ptycho2D objects.

  timestamp_increment(p)

    Increment the timestamp counter corresponding to the processing language used (OpenCL or CUDA) Virtual, must be derived.

    Parameters:

      w – the object (e.g. wavefront) the operator will be applied to.

    Returns:

class pynx.ptycho.ptycho.Ptycho(probe=None, obj=None, background=None, data=None, nb_frame_total=None)

  Class for two-dimensional ptychographic data: object, probe, and observed diffraction. This may include only part of the data from a larger dataset.

  Parameters:

    - probe – the starting estimate of the probe, as a complex 2D numpy array - can be 3D if modes are used. the probe should be centered in the center of the array.

    - obj – the starting estimate of the object, as a complex 2D numpy array - can be 3D if modes are used.

    - background – 2D array with the incoherent background.

    - data – the PtychoData object with all observed frames, ptycho positions

    - nb_frame_total – total number of frames (used for normalization)

  calc_regularisation_scale()

    Calculate the scale factor for object and probe regularisation. Calculated according to Thibault & Guizar-Sicairos 2012 :return: nothing

  calc_scan_area()

    Compute the scan area for the object and probe, using scipy ConvexHull. The scan area for the object is augmented by twice the average distance between scan positions for a more realistic estimation. scan_area_points is also computed, corresponding to the outline of the scanned area.

    Returns:

      Nothing. self.scan_area_probe and self.scan_area_obj are updated, as 2D arrays with the same shape as the object and probe, with False outside the scan area and True inside.

  from_pu()

    Get all relevant arrays from processing unit, if necessary :return: Nothing

  get_background()

    Get the background data array. This will automatically get the latest data, either from GPU or from the host memory, depending where the last changes were made.

    Returns:

      the 2D numpy data array

  get_illumination_obj()

    Get the sum of the probe intensity for all illuminations, which is used for object normalisation.

    Returns:

      the array of the illumination norm, with the same 2D shape as the object

  get_iobs(i=None, shift=False)

    Get the observed intensity array.

    Parameters:

      - i – the desired frame. if None, all data is returned.

      - shift – if True, the array will be shifted so that the array is centered on the frame center (and not the corner).

    Returns:

      a masked array of the observed intensity for frame i

  get_llk(noise=None, norm=True)

    Get the log-likelihood.

    Parameters:

      - noise – noise model, either ‘poisson’, ‘gaussian’ or ‘euclidian’. If None, a dictionary is returned.

      - norm – if True (the default), the LLK values are normalised

    Returns:

      either a single LLK value, or a dictionary

  get_obj(remove_obj_phase_ramp=False)

    Get the object data array. This will automatically get the latest data, either from GPU or from the host memory, depending where the last changes were made.

    Parameters:

      remove_obj_phase_ramp – if True, the object will be returned after removing the phase ramp coming from the imperfect centring of the diffraction data (sub-pixel shift). Calculated diffraction patterns using such a corrected object will present a sub-pixel shift relative to the diffraction data. The ramp information comes from the PtychoData phase_ramp_d{x,y} attributes, and should have been calculated beforehand using a ZroPhaseRamp(obj=True) operator.

    Returns:

      the 3D numpy data array (nb object modes, nyo, nxo)

  get_obj_coord()

    Get the object coordinates :return: a tuple of two arrays corresponding to the x (columns) and y (rows coordinates)

  get_probe()

    Get the probe data array. This will automatically get the latest data, either from GPU or from the host memory, depending where the last changes were made.

    Returns:

      the 3D probe numpy data array

  get_probe_coord()

    Get the probe coordinates :return: a tuple of two arrays corresponding to the x (columns) and y (rows coordinates)

  get_psi(shift=False)

    Get the current Psi array. Note that this is not not intended for an efficient data access, but only for occasional ones.

    Parameters:

      shift – if True, the array will be shifted so that the array is centered on the frame center (and not the corner).

    Returns:

      the raw psi numpy array for the currently processed stack of frames. Can also be None if no calculation took place.

  get_scan_area_obj()

    Return the scan_area_obj (2D array with the object shape, True inside the area scanned, and False outside). It is computed if necessary. :return: self.scan_area_obj

  get_scan_area_points()

    Return the scan_area_points (outside polygon points delimiting the scan area). It is computed if necessary. :return: self.scan_area_points

  get_scan_area_probe()

    Return the scan_area_probe (2D array with the probe shape, True inside the area scanned, and False outside). It is computed if necessary. :return: self.scan_area_probe

  load_obj_probe_cxi(filename, entry=None, verbose=True)

    Load object and probe from a CXI file, result of a previous optimisation. If no data is already present in the current object, then the pixel size and energy/wavelength are also loaded, and a dummy (one frame) data object is created.

    Parameters:

      - filename – the CXI filename from which to load the data

      - entry – the entry to be read. By default, the last in the file is loaded. Can be ‘entry_1’, ‘entry_10’… An integer n can also be supplied, in which case ‘entry_%d’ % n will be read

    Returns:

  print(*args, **kwargs)

    Print function which can be muted e.g. for non-master MPI processes :param args: arguments passed to the print function :param kwargs: keyword arguments passed to the print function :return: nothing

  reset_history()

    Reset history, and set current cycle to zero :return: nothing

  save_obj_probe_cxi(filename, sample_name=None, sample_nx_dict=None, experiment_id=None, instrument=None, note=None, process=None, append=False, shift_phase_zero=False, params=None, remove_obj_phase_ramp=False, cxi_data='object_probe')

    Save the result of the optimisation (object, probe, scan areas) to an HDF5 CXI-like file.

    Parameters:

      - filename – the file name to save the data to

      - sample_name – optional, the sample name. Ignored if sample_nx_dict is given.

      - sample_nx_dict – NXsample object as a dictionary, from silx.io.dictdump.nxtodict.

      - experiment_id – the string identifying the experiment, e.g.: ‘HC1234: Siemens star calibration tests’

      - instrument – the string identifying the instrument, e.g.: ‘ESRF id10’

      - note – a string with a text note giving some additional information about the data, a publication…

      - process – a dictionary of strings which will be saved in ‘/entry_N/data_1/process_1’. A dictionary entry can also be a ‘note’ as keyword and a dictionary as value - all key/values will then be saved as separate notes. Example: process={‘program’: ‘PyNX’, ‘note’:{‘llk’:1.056, ‘nb_photons’: 1e8}}

      - append – by default (append=False), any existing file will be overwritten, and the result will be saved as ‘entry_1’. If append==True and the file exists, a new entry_N will be saved instead. This can be used to export the different steps of an optimisation.

      - shift_phase_zero – if True, remove the linear phase ramp from the object

      - params – a dictionary of parameters to be saved into process_1/configuration NXcollection

      - remove_obj_phase_ramp – if True, the object will be saved after removing the phase ramp coming from the imperfect centring of the diffraction data (sub-pixel shift). Calculated diffraction patterns using such a corrected object will present a sub-pixel shift relative to the diffraction data. The ramp information comes from the PtychoData phase_ramp_d{x,y} attributes, and are not re-calculated.

      - cxi_data –

        a string describing what large data arrays are actually saved. This only affects the large object and probe arrays.

        - ’object_probe’: all object and probe modes are saved as complex objects.

        - ’object_phase’: only the first mode of the object’s phase is saved (as float16). This is the format which saves the most space

        - ’object’: save only the complex object

        - ’probe’: save only the complex probe

    Returns:

      Nothing. a CXI file is created

  set_background(background)

    Set the incoherent background data array. It will be shifted so that the center of the diffraction image is at (0,0), like the stored intensity.

    Parameters:

      background – the background (float32 numpy array)

    Returns:

      nothing

  set_obj(obj)

    Set the object data array.

    Parameters:

      obj – the object (complex64 numpy array)

    Returns:

      nothing

  set_obj_zero_phase_mask(mask)

    Set an object mask, which has the same 2D shape as the object, where values of 1 indicate that the area corresponds to vacuum (or air), and 0 corresponds to some material. Values between 0 and 1 can be given to smooth the transition. This mask will be used to restrain the corresponding area to a null phase, dampening the imaginary part at every object update. :param mask: a floating-point array with the same 2D shape as the object, where values of 1 indicate that the area corresponds to vacuum (or air), and 0 corresponds to the sample. :return: nothing

  set_probe(probe)

    Set the probe data array.

    Parameters:

      probe – the probe (complex64 numpy array)

    Returns:

      nothing

  update_history(mode='llk', update_obj=False, update_probe=False, update_background=False, update_pos=False, verbose=False, **kwargs)

    Update the history record.

    Parameters:

      - mode – either ‘llk’ (will record new log-likelihood and number of photons) or ‘algorithm’ (will only update the algorithm) - for the latter case, algorithm should be given as a keyword argument.

      - verbose – if True, print some info about current process (only if mode==’llk’)

      - kwargs – other parameters to be recorded, e.g. probe_inertia=xx, algorithm=’DM’

    Returns:

      nothing

class pynx.ptycho.ptycho.PtychoData(iobs=None, positions=None, detector_distance=None, mask=None, pixel_size_detector=None, wavelength=None, near_field=False, padding=0, vidx=None)

  Class for two-dimensional ptychographic data: observed diffraction and probe positions. This may include only part of the data from a larger dataset.

  Parameters:

    - iobs – 3d array with (nb_frame, ny,nx) shape observed intensity (assumed to follow Poisson statistics). The frames will be stored fft-shifted so that the center of diffraction lies in the (0,0) corner of each image. The supplied frames should have the diffraction center in the middle of the frames. Intensities must be >=0, except for masked pixels which should be <0. For the purpose of Paganin and CTF operators, masked values should be set to inter/extrapolated values and stored as -1-I_interp

    - positions – (x, y, z) tuple or 2-column array with ptycho probe positions in meters. For 2D data, z is ignored and can be None or missing, e.g. with (x, y) The orthonormal coordinate system follows the CXI/NeXus/McStas convention, with z along the incident beam propagation direction and y vertical towards ceiling. Values for x or y larger than 1e10 indicate that the direct beam is measured, allowing an absolute reference for the probe. Alternatively, those values may also be masked in a numpy masked array.

    - detector_distance – detector distance in meters

    - mask – obsolete, should not be used anymore. Masked intensity should be supplied as negative intensities.

    - pixel_size_detector – in meters, assuming square pixels

    - wavelength – wavelength of the experiment, in meters.

    - near_field – True if using near field ptycho

    - padding – an integer value indicating the number of zero-padded pixels to be used on each side of the observed frames. This can be used for near field ptychography. The input iobs should already padded, the corresponding pixels will be added to the mask.

    - vidx – array of indices of the positions. This is useful when only a subset of positions is studied, e.g. when the dataset is split between parallel processes.

    - kwargs – used to detect obsoleted parameters such as mask

  get_required_obj_shape(margin=16)

    Estimate the required object shape

    Parameters:

      margin – number of pixels on the border to avoid stepping out of the object, e.g. when refining positions.

    Returns:

      (nyo, nxo), the 2D object shape

  pixel_size_object()

    Get the x and y pixel size in object space after a FFT. :return: a tuple (pixel_size_x, pixel_size_y) in meters

pynx.ptycho.ptycho.algo_string(algo_base, p, update_object, update_probe, update_background=False, update_pos=False)

  Get a short string for the algorithm being run, e.g. ‘DM/o/3p’ for difference map with 1 object and 3 probe modes.

  Parameters:

    - algo_base – ‘AP’ or ‘ML’ or ‘DM’

    - p – the ptycho object

    - update_object – True if updating the object

    - update_probe – True if updating the probe

    - update_background – True if updating the background

    - update_pos – True if updating the positions

  Returns:

    a short string for the algorithm

pynx.ptycho.ptycho.calc_throughput(p: Ptycho = None, cxi=None, verbose=False)

  Analyse the throughput after a series of algorithms, either from a Ptycho object or from a CXI file. A few things like object & probe smoothing are not taken into account, or the use of bilinear interpolation, background… :param p: the Ptycho object the timings are extracted from. :param cxi: the CXI file the history of cycles will be obtained from :param verbose: if True, print the average throughput per algorithm step :return: the average throughput in Gbyte/s

pynx.ptycho.ptycho.save_obj_probe_cxi(filename, obj, probe, wavelength, detector_distance, pixel_size_detector, llk_poisson, llk_gaussian, llk_euclidian, nb_photons_calc, history, pixel_size_object, posxy, posxy_c, posxy0, scale=None, obj_zero_phase_mask=None, scan_area_obj=None, scan_area_probe=None, obj_illumination=None, background=None, sample_name=None, sample_nx_dict=None, experiment_id=None, instrument=None, note=None, process=None, append=False, shift_phase_zero=False, params=None, extra_data=None, obj_phase_ramp=None)

  Save the result of the optimisation (object, probe, scan areas) to an HDF5 CXI-like file.

  Note that object and probed are flipped (up/down) to have a (top, left) array origin.

  Parameters:

    - filename – the file name to save the data to

    - obj – the object to save

    - probe – the probe to save

    - wavelength – the wavelength (SI unit)

    - detector_distance – detector distance

    - pixel_size_detector – the detector’s pixel size

    - llk_euclidian (llk_poisson, llk_gaussian,) – normalised log-likelihood values

    - pixel_size_object – the object pixel size

    - posxy – the final scanning positions (SI units)

    - posxy0 – the initial scanning positions (SI units)

    - posxy_c – xy coordinates of the object center

    - scale – array of per-frame scaling factir

    - obj_zero_phase_mask – the area used for the object phase optimisation (maybe obsolete..)

    - scan_area_obj – scan area mask on the object array

    - scan_area_probe – scan area mask on the probe array

    - obj_illumination – integrated incident intensity on the object area

    - background – incoherent background

    - sample_name – optional, the sample name. Ignored if sample_nx_dict is given.

    - sample_nx_dict – NXsample object as a dictionary, from silx.io.dictdump.nxtodict.

    - experiment_id – the string identifying the experiment, e.g.: ‘HC1234: Siemens star calibration tests’

    - instrument – the string identifying the instrument, e.g.: ‘ESRF id10’

    - note – a string with a text note giving some additional information about the data, a publication…

    - process – a dictionary of strings which will be saved in ‘/entry_N/data_1/process_1’. A dictionary entry can also be a ‘note’ as keyword and a dictionary as value - all key/values will then be saved as separate notes. Example: process={‘program’: ‘PyNX’, ‘note’:{‘llk’:1.056, ‘nb_photons’: 1e8}}

    - append – by default (append=False), any existing file will be overwritten, and the result will be saved as ‘entry_1’. If append==True and the file exists, a new entry_N will be saved instead. This can be used to export the different steps of an optimisation.

    - shift_phase_zero – if True, centre the object phase around zero

    - params – a dictionary of parameters to be saved into process_1/configuration NXcollection

    - extra_data – a dictionary of data which will be saved as entry/_extra_data, and may be useful for debugging. each value may itself be a dictionary of values to save

    - obj_phase_ramp – the shifts (dx, dy) of the average calculated intensity from the array centre

  Returns:

    Nothing. a CXI file is created

pynx.ptycho.ptycho.save_obj_probe_nxtomo(filename, obj, probe, wavelength, detector_distance, pixel_size_detector, llk_poisson, llk_gaussian, llk_euclidian, nb_photons_calc, history, tomo_angle, crop_shape, pixel_size_object, posxy, posxy_c, posxy0, scale=None, obj_zero_phase_mask=None, scan_area_obj=None, scan_area_probe=None, obj_illumination=None, background=None, sample_name=None, experiment_id=None, instrument=None, note=None, process=None, append=False, shift_phase_zero=False, params=None, extra_data=None, obj_phase_ramp=None, deriv_phase=True, retry=True)

  Save the result of the optimisation (object, probe, scan areas) to an NXTomo (https://manual.nexusformat.org/classes/applications/NXtomo.html) file. This is used for ptycho-tomo experiments, with each projection (tomo angle) being written independently and possibly concurrently.

  Note that object and probed are flipped (up/down) to have a (top, left) array origin.

  Parameters:

    - filename – the file name to save the data to

    - obj – the object to save

    - probe – the probe to save (will not be saved if None)

    - wavelength – the wavelength (SI unit)

    - detector_distance – detector distance

    - pixel_size_detector – the detector’s pixel size

    - llk_euclidian (llk_poisson, llk_gaussian,) – normalised log-likelihood values

    - nb_photons_calc – total number of photons, calculated

    - history – the pynx.utils.history.History record of the optimisation

    - tomo_angle – the rotation angle for this projection, in radians

    - crop_shape – the 2D shape (ny, nx) the object should be cropped to before saving, so all projections have the same shape.

    - pixel_size_object – the object pixel size

    - posxy – the final scanning positions (SI units)

    - posxy0 – the initial scanning positions (SI units)

    - posxy_c – xy coordinates of the object center

    - background – incoherent background

    - sample_name – the sample name. If not given, the filename without extension will be used.

    - experiment_id – the string identifying the experiment, e.g.: ‘HC1234: Siemens star calibration tests’

    - instrument – the string identifying the instrument, e.g.: ‘ESRF id10’

    - note – a string with a text note giving some additional information about the data, a publication…

    - process – a dictionary of strings which will be saved in ‘/entry_N/data_1/process_1’. A dictionary entry can also be a ‘note’ as keyword and a dictionary as value - all key/values will then be saved as separate notes. Example: process={‘program’: ‘PyNX’, ‘note’:{‘llk’:1.056, ‘nb_photons’: 1e8}}

    - params – a dictionary of parameters to be saved into process_1/configuration NXcollection

    - extra_data – a dictionary of data which will be saved as entry/_extra_data, and may be useful for debugging. each value may itself be a dictionary of values to save

    - obj_phase_ramp – the shifts (dx, dy) of the average calculated intensity from the array centre

    - deriv_phase – if True, the phase derivative (along the horizontal direction) will be computed from the complex object.

    - retry – if True, the opening of the .nx file in append mode will be re-tried until successful - this is needed for independent, asynchronous parallel optimisations.

  Returns:

    Nothing. A nexus file is created or appended

pynx.ptycho.ptycho.save_ptycho_data_cxi(file_name, iobs, pixel_size, wavelength, detector_distance, x, y, z=None, monitor=None, mask=None, dark=None, instrument='', overwrite=False, scan=None, params=None, verbose=False, **kwargs)

  Save the Ptychography scan data using the CXI format (see http://cxidb.org)

  Parameters:

    - file_name – the file name (including relative or full path) to save the data to

    - iobs – the observed intensity, with shape (nb_frame, ny, nx)

    - pixel_size – the detector pixel size in meters

    - wavelength – the experiment wavelength

    - x – the x scan positions

    - y – the y scan positions

    - z – the z scan positions (default=None)

    - monitor – the monitor

    - mask – the mask for the observed frames

    - dark – the incoherent background

    - instrument – a string with the name of the instrument (e.g. ‘ESRF id16A’)

    - overwrite – if True, will overwrite an existing file

    - params – a dictionary of parameters which will be saved as a NXcollection

    - verbose – if True, print something.

  Returns:

### 2D Ptychography Operators

This section lists the operators, which can be imported automatically using from pynx.ptycho import * or from pynx.ptycho.operator import *. When this import is done, either CUDA (preferably) or OpenCL operators will be imported. The documentation below corresponds to OpenCL operators, but this should be identical to CUDA operators.

class pynx.ptycho.cl_operator.AP(update_object=True, update_probe=False, update_background=False, floating_intensity=False, nb_cycle=1, calc_llk=False, show_obj_probe=False, fig_num=-1, obj_smooth_sigma=0, obj_inertia=0.01, probe_smooth_sigma=0, probe_inertia=0.001, update_pos=False, pos_mult=1, pos_max_shift=2, pos_min_shift=0, pos_threshold=0.2, pos_history=False, zero_phase_ramp=True, background_smooth_sigma=0)

  Perform a complete Alternating Projection cycle: - forward all object*probe views to Fourier space and apply the observed amplitude - back-project to object space and project onto (probe, object) - update background optionally

  Parameters:

    - update_object – update object ?

    - update_probe – update probe ?

    - update_background – update background ?

    - floating_intensity – optimise floating intensity scale factor

    - nb_cycle – number of cycles to perform. Equivalent to AP(…)**nb_cycle

    - calc_llk – if True, calculate llk while in Fourier space. If a positive integer is given, llk will be calculated every calc_llk cycle

    - show_obj_probe – if a positive integer number N, the object & probe will be displayed every N cycle. By default 0 (no plot)

    - fig_num – the number of the figure to plot the object and probe, as for ShowObjProbe()

    - obj_smooth_sigma – if > 0, the previous object array (used for inertia) will convoluted by a gaussian array of this standard deviation.

    - obj_inertia – the updated object retains this relative amount of the previous object.

    - probe_smooth_sigma – if > 0, the previous probe array (used for inertia) will convoluted by a gaussian array of this standard deviation.

    - probe_inertia – the updated probe retains this relative amount of the previous probe.

    - update_pos – if True, update positions. This automatically inhibits object and probe update

    - pos_max_shift – maximum allowed shift (in pixels) per scan position (default=2)

    - pos_min_shift – minimum required shift (in pixels) per scan position (default=0)

    - pos_threshold – if the integrated grad_obj*probe along dx or dy is smaller than (grad_obj*probe).mean()*threshold, then the shift is ignored. This allows to prevent position shifts where there is little contrast in the object. It applies independently to x and y.

    - pos_mult – multiply the calculated position shifts by this value. Useful since the calculated shifts usually are a fraction of the actual shift.

    - pos_history – if True, save the position history (for debugging)

    - zero_phase_ramp – if True, the conjugate phase ramp in the object and probe will be removed by centring the FT of the probe, at the end and before every display. Ignored for near field.

    - background_smooth_sigma – gaussian convolution parameter for the background update

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.ApplyAmplitude(calc_llk=False, update_background=False, scale_in=1, scale_out=1, sum_icalc=False, background_smooth_sigma=0)

  Apply the magnitude from observed intensities, keep the phase. Applies to a stack of N views.

  Parameters:

    - calc_llk – if True, the log-likelihood will be calculated for this stack.

    - update_background – if True, update the background. The new background is automatically updated once the last stack is processed.

    - scale_in – a scale factor by which the input values should be multiplied, typically because of FFT

    - scale_out – a scale factor by which the output values should be multiplied, typically because of FFT

    - sum_icalc – if True, will sum the calculated intensity on each frame, and store this in p._cl_icalc_sum, to be used for floating intensities

    - background_smooth_sigma – sigma for gaussian smoothing of the background

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.ApplyPhaseRamp(dx, dy, probe=False, obj=False)

  Apply a given phase ramp to the object and/or probe. The actual phase factor is: For the probe: np.exp(-2j*np.pi*(x * dx + y*dy)) For the object: np.exp(2j*np.pi*(x * dx + y*dy)) Where (x,y) are reduced pixel coordinates going from -.5 to .5 for the probe, and +/-0.5 * nxo/nx (or nyo/ny) for the object.

  If the phase ramp is applied to both object and probe, the calculated intensity patterns remain unchanged.

  Parameters:

    - dy (dx,) – relative shifts from the centre, calculated in reciprocal space (probe array pixel coordinates)

    - obj – if True, apply the correction to the object

    - probe – if True, apply the correction to the probe

  op(p: Ptycho)

    Applies the operator to one object. Virtual function, must be derived. By default this is the identity operator.

    Returns:

      the result of the operation, usually the same type as the input, which is modified in-place. But this can also be a scalar (reduction).

class pynx.ptycho.cl_operator.BackgroundFilter(sigma)

  Apply a Gaussian filter to the background array. This operator will perform a FT on the backround, multiply it by a Gaussian, and back-FT the array. The resulting array is normalised so that its sum is equal to the original array sum.

  Parameters:

    sigma – standard deviation for the Gaussian kernel in the detector space, corresponding to a convolution of the background by a Gaussian of FWHM=2.3548*sigma

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.CLObsDataStack(cl_obs, cl_x, cl_y, i, npsi)

  Class to store a stack (e.g. 16 frames) of observed data in OpenCL space

  Parameters:

    - cl_obs – pyopencl array of observed data, with N frames

    - cl_y (cl_x,) – pyopencl arrays of the positions (in pixels) of the different frames

    - i – index of the first frame

    - npsi – number of valid frames (others are filled with zeros)

class pynx.ptycho.cl_operator.CLOperatorPtycho(processing_unit=None)

  Base class for a operators on Ptycho objects using OpenCL

  apply_ops_mul(pty)

    Apply the series of operators stored in self.ops to a Ptycho2D object. The operators are applied one after the other.

    Parameters:

      pty – the Ptycho2D to which the operators will be applied.

    Returns:

      the Ptycho2D object, after application of all the operators in sequence

  init_cl_vobs(p)

    Initialize observed intensity and scan positions in OpenCL space

    Parameters:

      p – the Ptycho object the operator will be applied to.

    Returns:

  prepare_data(p: Ptycho)

    Make sure the data to be used is in the correct memory (host or GPU) for the operator. Virtual, must be derived.

    Parameters:

      p – the Ptycho object the operator will be applied to.

    Returns:

  timestamp_increment(p)

    Increment the timestamp counter corresponding to the processing language used (OpenCL or CUDA) Virtual, must be derived.

    Parameters:

      w – the object (e.g. wavefront) the operator will be applied to.

    Returns:

  view_copy(pty, i_source, i_dest)

    Create a new view of the object by copying the original data. This will make a copy of all relevant data, which can be a wavefront, CDI object, Ptychography object, probe and psi arrais, etc…

    This (virtual) function is used to make temporary copies of the data Operators apply to. This is used to make linear combination (additions) of operators, which requires several copies of the data. As the copying part depends on the processing unit used (GPU, CPU) and the exact data to duplicate, this is a pure virtual class, which must be derived to be used. .. note:

        - this should only copy the 'active' data (which is affected by calculations)
        - index 0 corresponds to the original array, to which subsequent operators will be applied to

    Parameters:

      - obj – the object where the data will be duplicated

      - i_source – the index (integer) of the source object data

      - i_dest – the index (integer) of the destination object data

    Returns:

      nothing. The object is modified in-place

  view_purge(pty, i)

    Purge the different views of an object (except the main one). As the purging depends on the processing unit used (GPU, CPU) and the exact data to purge, this is a pure virtual function, which must be derived to be used. :param obj: the object where the view will be purged :param i_view: the index of the view to purge. If None, all views are purged. This de-registers the view. :return: nothing

  view_register(obj)

    Creates a new unique view key in an object. When finished with this view, it should be de-registered using view_purge. Note that it only reserves the key, but does not create the view. :return: an integer value, which corresponds to yet-unused key in the object’s view.

  view_sum(pty, i_source, i_dest)

    Add the view data from one index into another. As the summing depends on the processing unit used (GPU, CPU) and the exact data to sum, this is a pure virtual function, which must be derived to be used. :param obj: the object where the data will be duplicated :param i_source: the index (integer) of the source object data :param i_dest: the index (integer) of the destination object data :return: nothing. The object is modified in-place

  view_swap(pty, i1, i2)

    Swap the object view between index i1 and i2. As the swapping part depends on the processing unit used (GPU, CPU) and the exact data to swap, this is a pure virtual function, which must be derived to be used. :param obj: the object where the data will be duplicated :param i1: the index (integer) of the first object data :param i2: the index (integer) of the second object data :return: nothing. The object is modified in-place

class pynx.ptycho.cl_operator.CLOperatorPtychoPower(op, n)

class pynx.ptycho.cl_operator.CLOperatorPtychoSum(op1, op2)

class pynx.ptycho.cl_operator.CLProcessingUnitPtycho

  Processing unit in OpenCL space, for 2D Ptycho operations.

  Handles initializing the context and kernels.

  cl_init_kernel_n(n)

    Init kernels specifically written for reduction with N-sized array (N being normally the stack_size) :param n: the size for the float_n arrays used :return: Nothing.

  cl_init_kernels()

    Initialize opencl kernels :return: nothing

  set_stack_size(s)

    Change the number of frames which are stacked to perform all operations in //. If it is larger than the total number of frames, operators like AP, DM, ML will loop over all the stacks. :param s: an integer number (default=16) :return: nothing

class pynx.ptycho.cl_operator.Calc2Obs(nb_photons_per_frame=None, poisson_noise=False, scale_probe_direct_beam=True)

  Copy the calculated intensities to the observed ones, optionally with Poisson noise. This operator will loop other all stacks of frames, multiply object and probe and propagate the wavefront to the detector, and compute the new intensities.

  Parameters:

    - nb_photons_per_frame – average number of photons per frame, to scale the images. If None, no scaling is performed.

    - poisson_noise – if True, Poisson noise will be applied on the calculated frames. This uses numpy.random.poisson and is not GPU-optimised.

    - scale_probe_direct_beam – if True, poisson_noise is True, and the dataset includes at least one direct beam frame (without sample), then the probe will also be scaled so that it corresponds to the direct beam intensity.

  op(p: Ptycho)

    Applies the operator to one object. Virtual function, must be derived. By default this is the identity operator.

    Returns:

      the result of the operation, usually the same type as the input, which is modified in-place. But this can also be a scalar (reduction).

  timestamp_increment(p)

    Increment the timestamp counter corresponding to the processing language used (OpenCL or CUDA) Virtual, must be derived.

    Parameters:

      w – the object (e.g. wavefront) the operator will be applied to.

    Returns:

class pynx.ptycho.cl_operator.Calc2Obs1

  Copy the calculated intensities to the observed ones. Can be used for simulation. Applies to a stack of N views, assumes the current Psi is already in Fourier space.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.CalcIllumination(processing_unit=None)

  Compute the integrated illumination of the object by all probe positions

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

  timestamp_increment(p)

    Increment the timestamp counter corresponding to the processing language used (OpenCL or CUDA) Virtual, must be derived.

    Parameters:

      w – the object (e.g. wavefront) the operator will be applied to.

    Returns:

class pynx.ptycho.cl_operator.CenterObjProbe(max_shift=5, power=2, verbose=False)

  Operator to check the center of mass of the probe and shift both object and probe if necessary.

  Parameters:

    - max_shift – the maximum shift of the probe with respect to the center of the array, in pixels. The probe and object are only translated if the shift is larger than this value.

    - power – the center of mass is calculated on the amplitude of the array elevated at this power.

    - verbose – print deviation if verbose=True

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.DM(update_object=True, update_probe=True, update_background=False, nb_cycle=1, calc_llk=False, show_obj_probe=False, fig_num=-1, obj_smooth_sigma=0, obj_inertia=0.01, probe_smooth_sigma=0, probe_inertia=0.001, center_probe_n=0, center_probe_max_shift=5, loop_obj_probe=1, update_pos=False, pos_mult=1, pos_max_shift=2, pos_min_shift=0, pos_threshold=0.2, pos_history=False, zero_phase_ramp=True, background_smooth_sigma=0, alpha=0.02)

  Operator to perform a complete Difference Map cycle, updating the Psi views for all stack of frames, as well as updating the object and/or probe.

  Parameters:

    - update_object – update object ?

    - update_probe – update probe ?

    - update_background – update background ?

    - nb_cycle – number of cycles to perform. Equivalent to DM(…)**nb_cycle

    - calc_llk – if True, calculate llk while in Fourier space. If a positive integer is given, llk will be calculated every calc_llk cycle

    - show_obj_probe – if a positive integer number N, the object & probe will be displayed every N cycle. By default 0 (no plot)

    - fig_num – the number of the figure to plot the object and probe, as for ShowObjProbe()

    - obj_smooth_sigma – if > 0, the previous object array (used for inertia) will convoluted by a gaussian array of this standard deviation.

    - obj_inertia – the updated object retains this relative amount of the previous object.

    - probe_smooth_sigma – if > 0, the previous probe array (used for inertia) will convoluted by a gaussian array of this standard deviation.

    - probe_inertia – the updated probe retains this relative amount of the previous probe.

    - center_probe_n – test the probe every N cycle for deviation from the center. If deviation is larger than center_probe_max_shift, probe and object are shifted to correct. If 0 (the default), the probe centering is never calculated.

    - center_probe_max_shift – maximum deviation from the center (in pixels) to trigger a position correction

    - loop_obj_probe – when both object and probe are updated, it can be more stable to loop the object and probe update for a more stable optimisation, but slower.

    - update_pos – positive integer, if >0, update positions every ‘update_pos’ cycle. Note that object and probe are not updated during the same cycle as positions. Still experimental, recommended are 5, 10 (default=False, positions are not updated).

    - pos_max_shift – maximum allowed shift (in pixels) per scan position (default=2)

    - pos_min_shift – minimum required shift (in pixels) per scan position (default=0)

    - pos_threshold – if the integrated grad_obj*probe along dx or dy is smaller than (grad_obj*probe).mean()*threshold, then the shift is ignored. This allows to prevent position shifts where there is little contrast in the object. It applies independently to x and y.

    - pos_mult – multiply the calculated position shifts by this value. Useful since the calculated shifts usually are a fraction of the actual shift.

    - pos_history – if True, save the position history (for debugging)

    - zero_phase_ramp – if True, the conjugate phase ramp in the object and probe will be removed by centring the FT of the probe, at the end and before every display. Ignored for near field.

    - background_smooth_sigma – gaussian convolution parameter for the background update

    - alpha – mixing parameter between AP and DM, to bring more stability to DM. Default is 0.02, equivalent to 0.02*AP+0.98*DM.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.DM1(alpha=0.0)

  Equivalent to operator: (2-alpha) * ObjProbe2Psi() - (1-alpha)

  Parameters:

    alpha – alpha parameter to mix some AP within DM

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.DM2(alpha=0.0)

  # Psi(n+1) = (1-alpha) * (Psi(n) - P*O) + Psi_fourier

  This operator assumes that Psi_fourier is the current Psi, and that Psi(n) is in p._cl_psi_copy

  On output Psi(n+1) is the current Psi

  Parameters:

    alpha – alpha parameter to mix AP with DM

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.FT(scale=True, copy_psi=False)

  Forward Fourier-transform a Psi array, i.e. a stack of N Obj*Probe views

  Parameters:

    - scale – if True, the FFT will be normalized.

    - copy_psi – if True, perform the FT out-of-place, keeping the original Psi array in the Ptycho object’s _cl_psi_copy

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.FreePU(processing_unit=None)

  Operator freeing OpenCL memory. The fft plan/app in self.processing_unit is removed, as well as any OpenCL pyopencl.array.Array attribute in the supplied wavefront.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

  timestamp_increment(p)

    Increment the timestamp counter corresponding to the processing language used (OpenCL or CUDA) Virtual, must be derived.

    Parameters:

      w – the object (e.g. wavefront) the operator will be applied to.

    Returns:

class pynx.ptycho.cl_operator.Grad(update_object=True, update_probe=False, update_background=False, floating_intensity=False, reg_fac_obj=0, reg_fac_probe=0, calc_llk=False)

  Operator to compute the object and/or probe and/or background gradient. The gradient is stored in the ptycho object. It is assumed that the GPU gradient arrays have been already created, normally by the calling ML operator.

  Parameters:

    - update_object – compute gradient for the object ?

    - update_probe – compute gradient for the probe ?

    - update_background – compute gradient for the background ?

    - floating_intensity – optimise floating intensity scale factor

    - calc_llk – calculate llk while in Fourier space

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.IFT(scale=True)

  Backward Fourier-transform a Psi array, i.e. a stack of N Obj*Probe views

  Parameters:

    scale – if True, the FFT will be normalized.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.LLK(scale=False)

  Log-likelihood reduction kernel. Can only be used when Psi is in diffraction space. This is a reduction operator - it will write llk as an argument in the Ptycho object, and return the object. If _cl_stack_i==0, the llk is re-initialized. Otherwise it is added to the current value.

  The LLK can be calculated directly from object and probe using: p = LoopStack(LLK() * FT() * ObjProbe2Psi()) * p

  Parameters:

    scale – if True, will scale the calculated amplitude to calculate the log-likelihood. The amplitudes are left unchanged.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.LoopScan(op)

  Operator to apply a given operator sequentially to different scans,

  Parameters:

    op – the operator to apply, which can be a multiplication of operators

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.LoopStack(op, keep_psi=False, keep_psi_copy=False)

  Operator to apply a given operator sequentially to the complete stack of frames of a ptycho object.

  Make sure that the current selected stack is in a correct state (i.e. in sample or detector space,…) before starting such a loop with keep_psi=True.

  Parameters:

    - op – the operator to apply, which can be a multiplication of operators

    - keep_psi – if True, when switching between stacks, store psi in p._cl_psi_v.

    - keep_psi_copy – if True, when switching between stacks, store and restore psi_copy in p._cl_psi_copy_v. Used for algorithms (RAAR) requiring two psi copies across cycles.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.ML(nb_cycle=1, update_object=True, update_probe=False, update_background=False, floating_intensity=False, reg_fac_obj=0, reg_fac_probe=0, calc_llk=False, show_obj_probe=False, fig_num=-1, update_pos=False, pos_mult=1, pos_max_shift=2, pos_min_shift=0, pos_threshold=0.2, pos_history=False, zero_phase_ramp=True)

  Operator to perform a maximum-likelihood conjugate-gradient minimization.

  Parameters:

    - update_object – update object ?

    - update_probe – update probe ?

    - update_background – update background ?

    - floating_intensity – optimise floating intensity scale factor

    - reg_fac_obj – use this regularization factor for the object (if 0, no regularization)

    - reg_fac_probe – use this regularization factor for the probe (if 0, no regularization)

    - calc_llk – if True, calculate llk while in Fourier space. If a positive integer is given, llk will be calculated every calc_llk cycle

    - show_obj_probe – if a positive integer number N, the object & probe will be displayed every N cycle. By default 0 (no plot)

    - fig_num – the number of the figure to plot the object and probe, as for ShowObjProbe()

    - update_pos – positive integer, if >0, update positions every ‘update_pos’ cycle. Note that object and probe are not updated during the same cycle as positions. Still experimental, recommended are 5, 10 (default=False, positions are not updated).

    - pos_max_shift – maximum allowed shift (in pixels) per scan position (default=2)

    - pos_min_shift – minimum required shift (in pixels) per scan position (default=0)

    - pos_threshold – if the integrated grad_obj*probe along dx or dy is smaller than (grad_obj*probe).mean()*threshold, then the shift is ignored. This allows to prevent position shifts where there is little contrast in the object. It applies independently to x and y.

    - pos_mult – multiply the calculated position shifts by this value. Useful since the calculated shifts usually are a fraction of the actual shift.

    - pos_history – if True, save the position history (for debugging)

    - zero_phase_ramp – if True, the conjugate phase ramp in the object and probe will be removed by centring the FT of the probe, at the beginning and end. Ignored for near field.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.MemUsage(verbose=True)

  Print memory usage of current process (RSS on host) and used GPU memory

  op(p: Ptycho)
    Parameters:

      p – the ptycho object this operator applies to

    Returns:

      the updated ptycho object

  prepare_data(p: Ptycho)

    Make sure the data to be used is in the correct memory (host or GPU) for the operator. Virtual, must be derived.

    Parameters:

      p – the Ptycho object the operator will be applied to.

    Returns:

  timestamp_increment(p)

    Increment the timestamp counter corresponding to the processing language used (OpenCL or CUDA) Virtual, must be derived.

    Parameters:

      w – the object (e.g. wavefront) the operator will be applied to.

    Returns:

class pynx.ptycho.cl_operator.ObjProbe2LLK(use_copy=False)

  Compute the LLK from the current object and probe. This includes the object and probe multiplication and propagation

  Parameters:

    use_copy – if True, will work on the ptycho object’s _cl_psi_copy and leave _cl_psi untouched

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho with the computed LLK’s

class pynx.ptycho.cl_operator.ObjProbe2Psi(processing_unit=None)

  Computes Psi = Obj(r) * Probe(r-r_j) for a stack of N probe positions.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.PropagateApplyAmplitude(calc_llk=False, update_background=False, background_smooth_sigma=0, sum_icalc=False, copy_psi=False)

  Propagate to the detector plane (either in far or near field, perform the magnitude projection, and propagate back to the object plane. This applies to a stack of frames.

  Parameters:

    - calc_llk – if True, calculate llk while in the detector plane.

    - update_background – if >0, update the background with the difference between the observed and calculated intensity (with a damping factor), averaged over all frames.

    - background_smooth_sigma – sigma for the gaussian smoothing of the updated background

    - sum_icalc – if True, update the background.

    - copy_psi – if True, perform the propagation out-of-place, keeping the original Psi array in the Ptycho object’s _cl_psi_copy

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.PropagateNearField(forward=True, copy_psi=False)

  Near field propagator

  Parameters:

    - forward – if True, propagate forward, otherwise backward. The distance is taken from the ptycho data this operator applies to.

    - copy_psi – if True, perform the propagation out-of-place, keeping the original Psi array in the Ptycho object’s _cl_psi_copy

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.Psi2Obj(floating_intensity=False)

  Computes updated Obj(r) contributions from Psi and Probe(r-r_j), for a stack of N probe positions.

  Parameters:

    floating_intensity – if True, the intensity will be considered ‘floating’ from frame to frame, and a scale factor will be adjusted for each frame.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.Psi2ObjMerge(inertia=0.01, smooth_sigma=0, floating_intensity=False)

  Call this when all stack of probe positions have been processed, and the final update of the object can be calculated. Temporary arrays are cleaned up

  Parameters:

    - inertia – a regularisation factor to set the object inertia.

    - floating_intensity – if True, the floating scale factors will be corrected so that the average is 1.

    - smooth_sigma – if > 0, the previous object array (used for inertia) will be convolved by a gaussian with this sigma.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.Psi2PosMerge(multiplier=1, max_displ=2, min_displ=0, threshold=0.0, save_position_history=False)

  Merge scan position shifts, once the entire stqck of frames has been processed.

  Parameters:

    - multiplier – the computed displacements are multiplied by this value

    - max_displ – the displacements (at each iteration) are capped to this value (in pixels), after applying the multiplier.

    - min_displ – the displacements (at each iteration) are ignored if smaller than this value (in pixels), after applying the multiplier.

    - threshold – if the integrated grad_obj*probe along dx or dy is smaller than (grad_obj*probe).mean()*threshold, then the shift is ignored. This allows to prevent shifts where there is little contrast in the object. It applies independently to x and y.

    - save_position_history – if True, save the position history in the ptycho object (slow, for debugging)

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.Psi2PosShift(processing_unit=None)

  Computes scan position shifts, by comparing the updated Psi array to object*probe, for a stack of frames.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.Psi2Probe(processing_unit=None)

  Computes updated Probe contributions from Psi and Obj, for a stack of N probe positions.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.Psi2ProbeMerge(inertia=0.001, smooth_sigma=0)

  Call this when all stack of probe positions have been processed, and the final update of the probe can be calculated. Temporary arrays are cleaned up.

  Parameters:

    - inertia – a regularisation factor to set the object inertia.

    - smooth_sigma – if > 0, the previous object array (used for inertia) will be convolved by a gaussian with this sigma.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.PurgeStacks(processing_unit=None)

  Operator to delete stored psi stacks in a Ptycho object’s _cl_psi_v.

  This should be called for each main operator using LoopStack(), once it is finished processing, in order to avoid having another operator using the stored stacks, and to free memory.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.QuadraticPhase(factor, scale=1)

  Operator applying a quadratic phase factor

  Application of a quadratic phase factor, and optionally a scale factor.

  The actual factor is: \(scale * e^{i * factor * ((ix/nx)^2 + (iy/ny)^2)}\) where ix and iy are the integer indices of the pixels

  Parameters:

    - factor – the factor for the phase calculation.

    - scale – the data will be scaled by this factor. Useful to normalize before/after a Fourier transform, without accessing twice the array data.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.RAAR(update_object=True, update_probe=True, update_background=False, beta=0.9, nb_cycle=1, calc_llk=False, show_obj_probe=False, fig_num=-1, obj_smooth_sigma=0, obj_inertia=0.01, probe_smooth_sigma=0, probe_inertia=0.001, center_probe_n=0, center_probe_max_shift=5, loop_obj_probe=1, update_pos=False, pos_mult=1, pos_max_shift=2, pos_min_shift=0, pos_threshold=0.2, pos_history=False, zero_phase_ramp=True, background_smooth_sigma=0)

  Operator to perform a complete Relaxed Averaged Alternating Reflectors (RAAR) cycle, updating the Psi views for all stack of frames, as well as the object and/or probe.

  Parameters:

    - update_object – update object ?

    - update_probe – update probe ?

    - update_background – update background ?

    - beta – RAAR beta parameter.

    - nb_cycle – number of cycles to perform. Equivalent to DM(…)**nb_cycle

    - calc_llk – if True, calculate llk while in Fourier space. If a positive integer is given, llk will be calculated every calc_llk cycle

    - show_obj_probe – if a positive integer number N, the object & probe will be displayed every N cycle. By default 0 (no plot)

    - fig_num – the number of the figure to plot the object and probe, as for ShowObjProbe()

    - obj_smooth_sigma – if > 0, the previous object array (used for inertia) will convoluted by a gaussian array of this standard deviation.

    - obj_inertia – the updated object retains this relative amount of the previous object.

    - probe_smooth_sigma – if > 0, the previous probe array (used for inertia) will convoluted by a gaussian array of this standard deviation.

    - probe_inertia – the updated probe retains this relative amount of the previous probe.

    - center_probe_n – test the probe every N cycle for deviation from the center. If deviation is larger than center_probe_max_shift, probe and object are shifted to correct. If 0 (the default), the probe centering is never calculated.

    - center_probe_max_shift – maximum deviation from the center (in pixels) to trigger a position correction

    - loop_obj_probe – when both object and probe are updated, it can be more stable to loop the object and probe update for a more stable optimisation, but slower.

    - update_pos – positive integer, if >0, update positions every ‘update_pos’ cycle. Note that object and probe are not updated during the same cycle as positions. Still experimental, recommended are 5, 10 (default=False, positions are not updated).

    - pos_max_shift – maximum allowed shift (in pixels) per scan position (default=2)

    - pos_min_shift – minimum required shift (in pixels) per scan position (default=0)

    - pos_threshold – if the integrated grad_obj*probe along dx or dy is smaller than (grad_obj*probe).mean()*threshold, then the shift is ignored. This allows to prevent position shifts where there is little contrast in the object. It applies independently to x and y.

    - pos_mult – multiply the calculated position shifts by this value. Useful since the calculated shifts usually are a fraction of the actual shift.

    - pos_history – if True, save the position history (for debugging)

    - zero_phase_ramp – if True, the conjugate phase ramp in the object and probe will be removed by centring the FT of the probe, at the end and before every display. Ignored for near field.

    - background_smooth_sigma – gaussian convolution parameter for the background update

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.RAARProj(beta=0.9)

  # Psi(n+1) = beta/2 [2PO-Psi(n)+ Psi_Fourier] +(1-beta)Psi_Fourier

  This operator assumes that Psi_fourier is the current Psi, and that Psi(n) is in p._cl_psi_copy

  On output Psi(n+1) is the current Psi

  Parameters:

    alpha – beta RAAR parameter

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.Scale(x, obj=True, probe=True, psi=True)

  Multiply the ptycho object by a scalar (real or complex).

  Parameters:

    - x – the scaling factor

    - obj – if True, scale the object

    - probe – if True, scale the probe

    - psi – if True, scale the all the psi arrays, _cl_psi as well as _cl_psi_v

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.ScaleObjProbe(verbose=False, absolute=True)

  Operator to scale the object and probe so that they have the same magnitude, and that the product of object*probe matches the observed intensity (i.e. sum(abs(obj*probe)**2) = sum(iobs))

  Parameters:

    - verbose – print deviation if verbose=True

    - absolute – if True, the absolute scale is computed by comparing the calculated intensities with the observed ones. If False, only the relative scale of object and probe is modified (to avoid numerical divergence), not affecting the calculated intensities.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.SelectStack(stack_i, keep_psi=False, keep_psi_copy=False)

  Operator to select a stack of observed frames to work on. Note that once this operation has been applied, the new Psi value may be undefined (empty array), if no previous array existed.

  Select a new stack of frames, swapping data to store the last calculated psi array in the corresponding, ptycho object’s _cl_psi_v[i] dictionary.

  What happens is: * keep_psi=False: only the stack index in p is changed (p._cl_stack_i=stack_i)

  - keep_psi=True: the previous psi is stored in p._cl_psi_v[p._cl_stack_i], the new psi is swapped

      with p._cl_psi_v[stack_i] if it exists, otherwise initialized as an empty array.

  Parameters:

    - stack_i – the stack index.

    - keep_psi – if True, when switching between stacks, store and restore psi in p._cl_psi_v.

    - keep_psi_copy – if True, when switching between stacks, store and restore psi_copy in p._cl_psi_copy_v. Used for algorithms (RAAR) requiring two psi copies across cycles.

  op(p: Ptycho)
    Parameters:

      p – the Ptycho object this operator applies to

    Returns:

      the updated Ptycho object

class pynx.ptycho.cl_operator.SumIntensity1(icalc=None, iobs=None)

  Operator to compute the sum of calculated and/or observed frames. When calculating the sum of observed intensities, masked values will be replaced by calculated ones.

  This operator applies to a single stack, and will perform the object*probe multiplication, propagate to the detector, and compute the sums.

  Parameters:

    - icalc – the GPU array holding the sum of calculated intensities. If None, it is not calculated.

    - iobs – the GPU array holding the sum of observed intensities. If None, it is not calculated.

  op(p: Ptycho)

    Applies the operator to one object. Virtual function, must be derived. By default this is the identity operator.

    Returns:

      the result of the operation, usually the same type as the input, which is modified in-place. But this can also be a scalar (reduction).

class pynx.ptycho.cl_operator.ZeroPhaseRamp(obj=False)

  Operator to remove the linear phase ramp on both the object and the probe, for far field ptycho. This first computes the center of mass of the square norm of Fourier transform of the probe, and then corrects both probe and object for the phase ramp corresponding to the shift relative to the center of the array.

  Then, a remaining phase ramp in the object can optionally be computed by: - computing the sum of the calculated intensity along all frames - computing the shift of the center of mass of that intensity w/r to the array center - the phase ramp parameters are stored in the ptycho object, but not applied as the calculated

  Parameters:

    obj – if True, after correcting both object and probe from the probe phase ramp, the object phase ramp is evaluated from the center of the calculated diffraction.

  op(p: Ptycho)

    Applies the operator to one object. Virtual function, must be derived. By default this is the identity operator.

    Returns:

      the result of the operation, usually the same type as the input, which is modified in-place. But this can also be a scalar (reduction).

pynx.ptycho.cl_operator.get_icalc(p: Ptycho, i, shift=False)

  Get the calculated intensity for frame i and a given Ptycho object. Note that this will reset the current psi, so should not be used during algorithms like DM which rely on the previous Psi value.

  Parameters:

    - p – the Ptycho object

    - i – the index of the desired frame

    - shift – if True, the array will be shifted so that the array is centered on the frame center (and not the corner).

  Returns:

    the numpy array of the calculated intensity

### Ptychography Runner class

class pynx.ptycho.runner.runner.PtychoRunner(argv, params, *args, **kwargs)

  Class to process a series of scans with a series of algorithms, given from the command-line.

  Parameters:

    - argv – the command-line parameters.

    - params – parameters for the optimization, with some default values.

    - args – ignored. For backward compatibility only

    - kwargs – ignored. For backward compatibility only

  check_params()

    Check if self.params includes a minimal set of valid parameters

    Returns: Nothing. Will raise an exception if necessary

  check_params_beamline()

    Check if self.params includes a minimal set of valid parameters, specific to a beamline. Derived implementations can also set default values when appropriate.

    Returns: Nothing. Will raise an exception if necessary

  classmethod make_parser(default_par, script_name='pynx-ptycho', description=None, epilog=None)

    Create the parser for the command-line arguments. This should be updated in derived objects.

    Parameters:

      - default_par – dictionary of default parameters

      - script_name – name of the script, e.g. ‘pynx-ptycho-cxi’. By default, will be read from the command-line arguments.

      - description – initial description string (top of help text)

      - epilog – epilog for the parser (usually with examples)

    Returns:

      the parser

  parse_arg()

    Parses the arguments given on a command line

    Returns: nothing

  parse_parameters_file(filename)

    Read parameters from a text file, written with one parameter per line. The parameters are converted to command-line arguments (e.g. `data = myfile.h5` is converted to `--data myfile.h5` and appended to an argv-style list.

    Parameters:

      filename – the file to read the parameters from.

    Returns:

      an argv-style list with the parameters parsed from the file

  print(*args, **kwargs)

    MPI-aware print function. Non-master processes will be muted :param args: args passed to print :param kwargs: kwrags passed to print :return: nothing

  process_scans()

    Run all the analysis on the supplied scan list, unless ‘help’ is given as a command-line argument.

    Returns:

      Nothing

exception pynx.ptycho.runner.runner.PtychoRunnerException

class pynx.ptycho.runner.runner.PtychoRunnerScan(params, scan, mpi_comm=None, timings=None)

  Abstract class to handle ptychographic data. Must be derived to be used. Only the load_scan() and load_data() functions need be derived.

  center_crop_data()

    Once the data has been loaded in self.load_data() [overloaded in child classes), this function can be called at the end of load_data to take care of centering the data (finding the center of diffraction) and cropping it with a size suitable for clFFT. Rebin is performed if necessary.

    Returns:

      Nothing. self.iobs and self.dsize are updated. self.raw_data holds the raw data if needed for CXI export

  get_scan_prefix(run, scan=None)

    Get the save prefix given the parameters, scan and run numbers

    Parameters:

      - run – the run number

      - scan – the scan number - if None, self.scan is used

  init_run_number(_scan=None)

    Determine the current run number, depending on existing files :return: run0

  load_data()

    Loads data, using beamline-specific parameters. Abstract function, must be derived

    Returns:

  load_data_post_process()

    Applies some post-processing to the input data, according to parameters. Also loads the mask. User-supplied mask is loaded if necessary.

    This must be called at the end of load_data() :return:

  load_scan()

    Loads scan positions, using beamline-specific parameters. Abstract function, must be derived. This also filters the set of scan positions according to parameters (xyrange, monitor values, maxframe,…)

    If MPI is used, only the master reads the scan positions. This automatically calls mpi_scan_split() at the end to split the scan if necessary (when MPI is used and mpi=’splitscan’ is used).

    Returns: Nothing. The scan positions, and the scan position indices ()

      to be loaded for that runner are stored in self.x, self.y, self.imgn. If MPI is used, only the master should call this, and then call mpi_scan_split(). This is handled in Runner:process_scan().

  load_scan_post()

    Take care of frames with the direct beam illumination, mask positions. Should be called right after load_scan

  mpi_scan_split()

    This function is called after load_scan(), and will split the scan among all the MPI process. If MPI is not used, it does nothing and just passes the (x, y, imgn) values.

    Returns:

      nothing. The x, y, imgn attributes are updated if necessary after splitting among MPI processes

  prepare()

    Prepare object and probe.

    Returns: nothing. Adds self.obj0 and self.probe0

  prepare_processing_unit()

    Prepare processing unit (CUDA, OpenCL, or CPU). This must be called after load_scan so that the size of the dataset is known

    Returns: nothing. Creates self.processing_unit, and adapts the stack size for CUDA

  print(*args, **kwargs)

    MPI-aware print function. Non-master processes will be muted :param args: args passed to print :param kwargs: kwrags passed to print :return: nothing

  print_probe_fwhm()

    Analyse probe shape and print estimated FWHM. Ignored for near field ptycho.

    Returns:

      Nothing

  run(reuse_ptycho=False, share_probe_scan_data=None)

    Main work function, will run according to the set of algorithms specified in self.params.

    Parameters:

      - reuse_ptycho – if True, will reuse the previous Ptycho and PtychoData objects and skip some initialisation steps

      - share_probe_scan_data – used when processing jointly multiple scans with a shared probe

    Returns:

  run_algorithm(algo_string)

    Run a single or suite of algorithms in a given run.

    Parameters:

      algo_string – a single or suite of algorithm steps to use, e.g. ‘ML**100’ or ‘analysis,ML**100,DM**100,nbprobe=3,DM**100’ or ‘analysis,ML**100*DM**100,nbprobe=3,DM**100’

    Returns:

      Nothing

  save(run, stepnum=None, algostring=None, scan=None)

    Save the result of the optimization, and (if self.params[‘saveplot’] is True) the corresponding plot. This is an internal function.

    Parameters:

      - run – the run number (integer)

      - stepnum – the step number in the set of algorithm steps

      - algostring – the string corresponding to all the algorithms ran so far, e.g. ‘100DM,100AP,100ML’

      - scan – the scan number - normally self.scan, but can be used when processing simultaneously multiple scans

    Returns:

  save_data_cxi(crop=True, verbose=False, **kwargs)

    Save the scan data using the CXI format (see http://cxidb.org) :param crop: if True, only the already-cropped data will be save. Otherwise, the original raw data is saved.

    Returns:

  save_movie()

    Create a movie of the scan with all scan positions and diffraction frames. Requires matplotlib and ffmpeg. :return:

  save_plot(run, stepnum=None, algostring=None, display_plot=False)

    Save the plot to a png file.

    Parameters:

      - run – the run number (integer)

      - stepnum – the step number in the set of algorithm steps

      - algostring – the string corresponding to all the algorithms ran so far, e.g. ‘100DM,100AP,100ML’

      - display_plot – if True, the saved plot will also be displayed

    Returns:

  test_no_rerun(verbose=True)

    Test if this scan # should be skipped; :param scan: the scan number :return: True if the scan should be skipped, False otherwise

### Ptychography Analysis

pynx.ptycho.analysis.modes(d, pixel_size, do_plot=True, show_plot=True, verbose=False)

  Determine complex modes for a given object or probe.

  Parameters:

    - d – the probe or object to analyse, with the modes along the first axis/

    - pixel_size – the pixel size in meters

    - do_plot – plot if True (default)

    - show_plot – if True, the plot will be done in a visible figure, otherwise in an offline one. Ignored unless do_plot is True.

  Returns:

    a tuple of (orthogonalized modes with the same shape as the input, a vector of the relative intensities, figure)

pynx.ptycho.analysis.probe_fwhm(probe, pixel_size, verbose=True)

  Analyse probe shape and estimated FWHM using different methods:

  1. Full width at half maximum

  2. Full width at 20%

  3. Full width at half maximum using a statistical gaussian analysis - the width is the same as a gaussian with a maximum intensity of 1 and the same integrated intensity.

  Parameters:

    - probe – the probe to analyse. If number of dimensions is 3, only the first mdoe is analysed

    - pixel_size – the pixel size in meters

    - verbose – if True (the default), will print the corresponding sizes

  Returns:

    ((fwhm_x, fwhm_y), (fwhm_x@20%, fwhm_y@20%), (fwhm_gauss_stat, fwhm_gauss_stat_x, fwhm_gauss_stat_y))

pynx.ptycho.analysis.probe_propagate(probe, propagation_range, pixel_size, wavelength, do_plot=True, show_plot=True, fig_size=(15, 10))

  Propagate a probe along a given range, plot it and return the probe at the focus (position where the size is minimal)

  Parameters:

    - probe – the 2D complex probe to propagate

    - propagation_range – either given as a range/array of dz values (in meters), or a tuple with (zmin,zmax,nb)

    - pixel_size – the probe pixel size, in meters

    - do_plot – plot if True (default)

    - show_plot – if True, the plot will be done in a visible figure, otherwise in an offline one. Ignored unless do_plot is True

    - fig_size – the figure size, default = (15,10)

  Returns:

    - the 3d array of the propagated probe along z, with size (nz, ny, nx)

    - the z coordinates along the 3d array

    - the index of the found focus point along the z direction

    - the figure if it was plotted

  Return type:

    A tuple with

## Examples

### Operator-based API, far field

    ########################################################################
    #
    # Example of the ptychograpic reconstruction using OpenCL on simulated data
    # (c) ESRF 2017-present
    # Authors: Vincent Favre-Nicolin <favre@esrf.fr>
    #
    ########################################################################

    import timeit
    from pylab import *

    if False:
        from pynx.processing_unit import default_processing_unit

        # This can be used to select the GPU and/or the language, and must be called before other pynx imports
        # Otherwise a default GPU will be used according to its speed, which is usually sufficient
        default_processing_unit.select_gpu(language='OpenCL', gpu_name='R9')

    # The following statement will import either CUDA or OpenCL operators
    from pynx.ptycho import *
    from pynx.ptycho import simulation

    ##################
    # Simulation of the ptychographic data:
    n = 256
    pixel_size_detector = 55e-6
    wavelength = 1.5e-10
    detector_distance = 1
    obj_info = {'type': 'phase_ampl', 'phase_stretch': pi / 2, 'alpha_win': .2}
    # probe_info = {'type': 'focus', 'aperture': (30e-6, 30e-6), 'focal_length': .08, 'defocus': 100e-6, 'shape': (n, n)}
    probe_info = {'type': 'gauss', 'sigma_pix': (40, 40), 'defocus': 100e-6, 'shape': (n, n)}

    # 50 scan positions correspond to 4 turns, 78 to 5 turns, 113 to 6 turns
    scan_info = {'type': 'spiral', 'scan_step_pix': 30, 'n_scans': 120}
    data_info = {'num_phot_max': 1e9, 'bg': 0, 'wavelength': wavelength, 'detector_distance': detector_distance,
                 'detector_pixel_size': pixel_size_detector, 'noise': 'poisson'}

    # Initialisation of the simulation with specified parameters
    s = simulation.Simulation(obj_info=obj_info, probe_info=probe_info, scan_info=scan_info, data_info=data_info)
    s.make_data()

    # Positions from simulation are given in pixels
    posx, posy = s.scan.values

    ampl = s.amplitude.values  # square root of the measured diffraction pattern intensity

    pixel_size_object = wavelength * detector_distance / pixel_size_detector / n

    ##################
    # Size of the reconstructed object (obj)
    nyo, nxo = shape.calc_obj_shape(posx, posy, ampl.shape[1:])

    # Initial object
    # obj_init_info = {'type':'flat','shape':(nx,ny)}
    obj_init_info = {'type': 'random', 'range': (0, 1, 0, 0.5), 'shape': (nyo, nxo)}
    # Initial probe
    probe_init_info = {'type': 'focus', 'aperture': (20e-6, 20e-6), 'focal_length': .08, 'defocus': 50e-6, 'shape': (n, n)}
    data_info = {'wavelength': wavelength, 'detector_distance': detector_distance,
                 'detector_pixel_size': pixel_size_detector}
    init = simulation.Simulation(obj_info=obj_init_info, probe_info=probe_init_info, data_info=data_info)

    init.make_obj()
    init.make_probe()

    data = PtychoData(iobs=ampl ** 2, positions=(posx * pixel_size_object, posy * pixel_size_object), detector_distance=1,
                      mask=None, pixel_size_detector=55e-6, wavelength=1.5e-10)

    p = Ptycho(probe=s.probe.values, obj=init.obj.values, data=data, background=None)

    # Initial scaling is important to avoid overflows during ML
    p = ScaleObjProbe(verbose=True) * p
    p = DM(update_object=True, update_probe=True, calc_llk=20, show_obj_probe=20) ** 100 * p
    p = ML(update_object=True, update_probe=True, calc_llk=20, show_obj_probe=20) ** 40 * p

    if True:
        print("###########################################################################################################")
        # Add probe modes
        pr = p.get_probe()
        ny, nx = pr.shape[-2:]
        nb_probe = 3
        pr1 = np.empty((nb_probe, ny, nx), dtype=np.complex64)
        pr1[0] = pr[0]
        for i in range(1, nb_probe):
            n = abs(pr).mean()
            pr1[i] = np.random.uniform(0, n, (ny, nx)) * exp(1j * np.random.uniform(0, 2 * np.pi, (ny, nx)))

        p.set_probe(pr1)

        p = DM(update_object=True, update_probe=True, calc_llk=20, show_obj_probe=20) ** 100 * p
        p = AP(update_object=True, update_probe=True, calc_llk=20, show_obj_probe=20) ** 40 * p
        p = ML(update_object=True, update_probe=True, calc_llk=20, show_obj_probe=20) ** 100 * p

    if False:
        # Timing vs stack size
        n = 50
        vx = []
        vy = []
        for stack_size in range(1, 32 + 1, 1):
            default_processing_unit.set_stack_size(stack_size)
            p = DM(update_object=True, update_probe=True) ** 10 * p
            t0 = timeit.default_timer()
            p = DM(update_object=True, update_probe=True) ** n * p
            dt = (timeit.default_timer() - t0) / n
            print("DM dt/cycle=%5.3fs [stack_size=%2d]" % (dt, stack_size))
            vx.append(stack_size)
            vy.append(dt)
        figure()
        plot(vx, vy, '-')
        xlabel('stack size')
        ylabel('Time for a DM cycle (s)')

    if False:
        # Export data and object, probe to CXI files
        p.save_obj_probe_cxi('obj_probe.cxi')
        save_ptycho_data_cxi('data.cxi', ampl ** 2, pixel_size_detector, wavelength, detector_distance,
                             posx * pixel_size_object, posy * pixel_size_object, z=None, monitor=None,
                             mask=None, instrument='simulation', overwrite=True)

### Operator-based API, near field

    ########################################################################
    #
    # Example of the ptychograpic reconstruction using OpenCL on simulated data
    # (c) ESRF 2017-present
    # Authors: Vincent Favre-Nicolin <favre@esrf.fr>
    #
    ########################################################################

    import timeit
    from pylab import *
    from pynx.ptycho import simulation, shape

    # from pynx.ptycho import *
    from pynx.ptycho.ptycho import *
    from pynx.ptycho.cpu_operator import *  # Use only CPU operators ?
    from pynx.ptycho.operator import *  # Use CUDA > OpenCL > CPU operators, as available
    import pynx.ptycho.cpu_operator as cpuop
    import pynx.ptycho.cl_operator as clop

    ##################
    detector_distance = 1.5
    wavelength = 1.5e-10
    pixel_size_detector = 1e-6
    # Simulation of the ptychographic data:
    n = 256
    obj_info = {'type': 'phase_ampl', 'phase_stretch': pi / 2, 'alpha_win': .2}
    probe_info = {'type': 'near_field', 'aperture': (80e-6, 80e-6), 'defocus': 0.5, 'shape': (n, n)}

    # 50 scan positions correspond to 4 turns, 78 to 5 turns, 113 to 6 turns
    scan_info = {'type': 'spiral', 'scan_step_pix': 20, 'n_scans': 120}
    data_info = {'num_phot_max': 1e9, 'bg': 0, 'wavelength': wavelength, 'detector_distance': detector_distance,
                 'detector_pixel_size': pixel_size_detector, 'noise': 'poisson', 'near_field': True}

    # Initialisation of the simulation
    s = simulation.Simulation(obj_info=obj_info, probe_info=probe_info, scan_info=scan_info, data_info=data_info)
    s.make_data()

    # Positions from simulation are given in pixels
    posx, posy = s.scan.values
    ampl = s.amplitude.values  # square root of the measured diffraction pattern intensity

    ##################
    # Size of the reconstructed object (obj)
    nyo, nxo = shape.calc_obj_shape(posx, posy, ampl.shape[1:])

    # Initial object
    # obj_init_info = {'type':'flat','shape':(nx,ny)}
    obj_init_info = {'type': 'random', 'range': (0, 1, 0, 0.5), 'shape': (nyo, nxo)}
    # Initial probe
    probe_init_info = {'type': 'near_field', 'aperture': (90e-6, 90e-6), 'defocus': 0.3, 'shape': (n, n)}
    init = simulation.Simulation(obj_info=obj_init_info, probe_info=probe_init_info, data_info=data_info)

    init.make_obj()
    init.make_probe()

    data = PtychoData(iobs=ampl ** 2, positions=(posx * pixel_size_detector, posy * pixel_size_detector),
                      detector_distance=detector_distance, mask=None,
                      pixel_size_detector=pixel_size_detector, wavelength=wavelength, near_field=True)

    p = Ptycho(probe=init.probe.values, obj=init.obj.values, data=data, background=None)

    # Initial scaling is important to avoid overflows during ML
    p = ScaleObjProbe(verbose=True) * p
    p = DM(update_object=True, update_probe=True, calc_llk=20, show_obj_probe=20) ** 100 * p
    p = AP(update_object=True, update_probe=True, calc_llk=20, show_obj_probe=20) ** 40 * p
    p = ML(update_object=True, update_probe=True, calc_llk=20, show_obj_probe=20) ** 40 * p

    if False:
        # Timing vs stack size
        n = 50
        vx = []
        vy = []
        for stack_size in range(1, 32 + 1, 1):
            default_processing_unit.set_stack_size(stack_size)
            p = DM(update_object=True, update_probe=True) ** 10 * p
            t0 = timeit.default_timer()
            p = DM(update_object=True, update_probe=True) ** n * p
            dt = (timeit.default_timer() - t0) / n
            print("DM dt/cycle=%5.3fs [stack_size=%2d]" % (dt, stack_size))
            vx.append(stack_size)
            vy.append(dt)
        figure()
        plot(vx, vy, '-')
        xlabel('stack size')
        ylabel('Time for a DM cycle (s)')