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 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 analyze 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

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_obj()

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:shift – if True, the data array will be fft-shifted so that the center of the data is in the center of the array, rather than in the corner (the default).
Returns:the 3D numpy data array (nb object modes, nyo, nxo)
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
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:

reset_history()

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

save_obj_probe_cxi(filename, sample_name=None, experiment_id=None, instrument=None, note=None, process=None, append=False, shift_phase_zero=False, params=None)

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
  • 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
Returns:

Nothing. a CXI file is created

set_background(background)

Set the incoherent background data array.

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_backgroung=False, update_pos=False, verbose=False, **kwargs)

Update the history record. :param 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
Parameters:
  • verbose – if True, print some info about current process (only if mode==’llk’)
  • kwargs – other parameters to be recorded, e.g. probe_inertia=xx, dt=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)

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. Negative values will be used to mark masked pixels.
  • 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.
  • detector_distance – detector distance in meters
  • mask – 2D mask (>0 means masked pixel) for the observed data. Can be None. Will be fft-shifted like iobs. Masked pixels are stored as negative values in the iobs array.
  • 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.
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.save_ptycho_data_cxi(file_name, iobs, pixel_size, wavelength, detector_distance, x, y, z=None, monitor=None, mask=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
  • 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, update_pos_mult=1, update_pos_max_shift=2, update_pos_history=False)

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
  • update_pos_max_shift – maximum allowed shift (in pixels) per scan position (default=2)
  • update_pos_mult – multiply the calculated position shifts by this value. Useful since the calculated shifts usually are a fraction of the actual shift.
  • update_pos_history – if True, save the position history (for debugging)
op(p: pynx.ptycho.ptycho.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, sum_icalc=False)

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 according to Marchesini et el, Inverse Problems 29(11), 115009 (2013). This actually updates arrays from which the final background update can be evaluated. Note that the temporary arrays (_cl_vd, _cl_vd2, _cl_vz2, _cl_vdz2) must be created first in the Ptycho object beforehand (e.g. in the calling AP() operator), where the new background will be calculated once all stack of frames are processed.
  • 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
op(p: pynx.ptycho.ptycho.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)

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

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: pynx.ptycho.ptycho.Ptycho)
Parameters:p – the Ptycho object this operator applies to
Returns:the updated Ptycho object
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: pynx.ptycho.ptycho.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, 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, update_pos_mult=1, update_pos_max_shift=2, update_pos_history=False)

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 ?
  • 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).
  • update_pos_max_shift – maximum allowed shift (in pixels) per scan position (default=2)
  • update_pos_mult – multiply the calculated position shifts by this value. Useful since the calculated shifts usually are a fraction of the actual shift.
  • update_pos_history – if True, save the position history (for debugging)
op(p: pynx.ptycho.ptycho.Ptycho)
Parameters:p – the Ptycho object this operator applies to
Returns:the updated Ptycho object
class pynx.ptycho.cl_operator.DM1(processing_unit=None)

Equivalent to operator: 2 * ObjProbe2Psi() - 1

op(p: pynx.ptycho.ptycho.Ptycho)
Parameters:p – the Ptycho object this operator applies to
Returns:the updated Ptycho object
class pynx.ptycho.cl_operator.DM2(processing_unit=None)

# Psi(n+1) = 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_v

On output Psi(n+1) is the current Psi, and Psi_fourier has been swapped to p._cl_psi_v

op(p: pynx.ptycho.ptycho.Ptycho)
Parameters:p – the Ptycho object this operator applies to
Returns:the updated Ptycho object
class pynx.ptycho.cl_operator.FT(scale=True)

Forward 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: pynx.ptycho.ptycho.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 gpyfft data reference in self.processing_unit is removed, as well as any OpenCL pyopencl.array.Array attribute in the supplied wavefront.

op(p: pynx.ptycho.ptycho.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: pynx.ptycho.ptycho.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: pynx.ptycho.ptycho.Ptycho)
Parameters:p – the Ptycho object this operator applies to
Returns:the updated Ptycho object
class pynx.ptycho.cl_operator.LLK(processing_unit=None)

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

op(p: pynx.ptycho.ptycho.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, 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.
  • copy – make a copy of the original p._cl_psi swapped in as p._cl_psi_copy, and delete it after applying the operations. This is useful for operations requiring the previous value.
op(p: pynx.ptycho.ptycho.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, update_pos_mult=1, update_pos_max_shift=2, update_pos_history=False)

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).
  • update_pos_max_shift – maximum allowed shift (in pixels) per scan position (default=2)
  • update_pos_mult – multiply the calculated position shifts by this value. Useful since the calculated shifts usually are a fraction of the actual shift.
  • update_pos_history – if True, save the position history (for debugging)
op(p: pynx.ptycho.ptycho.Ptycho)
Parameters:p – the Ptycho object this operator applies to
Returns:the updated Ptycho object
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: pynx.ptycho.ptycho.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, sum_icalc=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 True, update the background according to Marchesini et el, Inverse Problems 29(11), 115009 (2013).
  • 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
op(p: pynx.ptycho.ptycho.Ptycho)
Parameters:p – the Ptycho object this operator applies to
Returns:the updated Ptycho object
class pynx.ptycho.cl_operator.PropagateNearField(forward=True)

Near field propagator

Parameters:forward – if True, propagate forward, otherwise backward. The distance is taken from the ptycho data this operator applies to.
op(p: pynx.ptycho.ptycho.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: pynx.ptycho.ptycho.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: pynx.ptycho.ptycho.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, 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 are capped to this value
  • save_position_history – if True, save the position history in the ptycho object (slow, for debugging)
op(p: pynx.ptycho.ptycho.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: pynx.ptycho.ptycho.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: pynx.ptycho.ptycho.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: pynx.ptycho.ptycho.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: pynx.ptycho.ptycho.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: pynx.ptycho.ptycho.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: pynx.ptycho.ptycho.Ptycho)
Parameters:p – the Ptycho object this operator applies to
Returns:the updated Ptycho object
class pynx.ptycho.cl_operator.ScaleObjProbe(verbose=False)

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
op(p: pynx.ptycho.ptycho.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)

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.
op(p: pynx.ptycho.ptycho.Ptycho)
Parameters:p – the Ptycho object this operator applies to
Returns:the updated Ptycho object

Ptychography Runner class

class pynx.ptycho.runner.runner.PtychoRunner(argv, params, ptycho_runner_scan_class)

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.
  • ptycho_runner_scan_class – the class to use to run the analysis.
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

parse_arg()

Parses the arguments given on a command line

Returns: nothing

parse_arg_beamline(k, v)

Parse argument in a beamline-specific way. This function only parses single arguments. If an argument is recognized and interpreted, the corresponding value is added to self.params

This method should be superseded in a beamline/instrument-specific child class.

Returns:True if the argument is interpreted, false otherwise
process_scans()

Run all the analysis on the supplied scan list

Returns:Nothing
exception pynx.ptycho.runner.runner.PtychoRunnerException
class pynx.ptycho.runner.runner.PtychoRunnerScan(params, scan)

Abstract class to handle ptychographic data. Must be derived to be used.

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
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:

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_data so that the size of the dataset is known

Returns: nothing. Creates self.processing_unit

print_probe_fwhm()

Analyze probe shape and print estimated FWHM

Returns:Nothing
run()

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

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)

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’
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:

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 analyze, 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)

Analyze 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 analyze. If number of dimensions is 3, only the first mdoe is analyzed
  • 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=10, show_obj_probe=10) ** 100 * p
p = ML(update_object=True, update_probe=True, calc_llk=10, show_obj_probe=10) ** 40 * p

if True:
    print("###########################################################################################################")
    # Add probe modes
    pr = p.get_probe()
    nb_probe, ny, nx = pr.shape
    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=10, show_obj_probe=10) ** 100 * p
    p = AP(update_object=True, update_probe=True, calc_llk=10, show_obj_probe=10) ** 40 * p
    p = ML(update_object=True, update_probe=True, calc_llk=10, show_obj_probe=10) ** 100 * p

if False:
    p = ShowObjProbe(fig_num=1) * DM(update_object=True, update_probe=False, calc_llk=10, show_obj_probe=10) ** 40 * p
    p = ShowObjProbe(fig_num=1) * DM(update_object=True, update_probe=True, calc_llk=10, show_obj_probe=10) ** 20 * p
    p = ShowObjProbe(fig_num=1) * AP(update_object=True, update_probe=True, calc_llk=10, show_obj_probe=10) ** 200 * p
    # p = ScaleObjProbe(verbose=True) * p  # Important to avoid overflow during ML
    # p = ShowObjProbe(fig_num=1) * DM(update_object=True, update_probe=True, calc_llk=1) ** 100 * p
    p = ShowObjProbe(fig_num=1) * ML(update_object=True, update_probe=True, calc_llk=10, show_obj_probe=10) ** 40 * p

if False:
    n = 50
    p = DM(update_object=True) ** 10 * p
    t0 = timeit.default_timer()
    p = DM(update_object=True) ** n * p
    print("DM dt/cycle=%5.3fs" % ((timeit.default_timer() - t0) / n))
    p = AP(update_object=True) ** 10 * p
    t0 = timeit.default_timer()
    p = AP(update_object=True) ** n * p
    print("AP dt/cycle=%5.3fs" % ((timeit.default_timer() - t0) / n))
    p = ML(update_object=True) ** 10 * p
    t0 = timeit.default_timer()
    p = ML(update_object=True) ** n * p
    print("ML dt/cycle=%5.3fs" % ((timeit.default_timer() - t0) / n))

if False:
    # Timing vs stack size
    n = 50
    vx = []
    vy = []
    for stack_size in range(1, 32 + 1, 1):
        default_processing_unit.cl_stack_size = np.int32(stack_size)
        p = DM(update_object=True) ** 10 * p
        t0 = timeit.default_timer()
        p = DM(update_object=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('Timd 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 = ShowObjProbe(fig_num=1) * DM(update_object=True, update_probe=True, calc_llk=10, show_obj_probe=10) ** 100 * p
p = ShowObjProbe(fig_num=1) * AP(update_object=True, update_probe=True, calc_llk=10, show_obj_probe=10) ** 40 * p
p = ShowObjProbe(fig_num=1) * ML(update_object=True, update_probe=True, calc_llk=10, show_obj_probe=10) ** 40 * p

if False:
    p = ShowObjProbe(fig_num=1) * DM(update_object=True, update_probe=False, calc_llk=10, show_obj_probe=10) ** 40 * p
    p = ShowObjProbe(fig_num=1) * DM(update_object=True, update_probe=True, calc_llk=10, show_obj_probe=10) ** 20 * p
    p = ShowObjProbe(fig_num=1) * AP(update_object=True, update_probe=True, calc_llk=10, show_obj_probe=10) ** 200 * p
    # p = ScaleObjProbe(verbose=True) * p  # Important to avoid overflow during ML
    # p = ShowObjProbe(fig_num=1) * DM(update_object=True, update_probe=True, calc_llk=1) ** 100 * p
    p = ShowObjProbe(fig_num=1) * ML(update_object=True, update_probe=True, calc_llk=10, show_obj_probe=10) ** 40 * p

if False:
    n = 50
    p = DM(update_object=True) ** 10 * p
    t0 = timeit.default_timer()
    p = DM(update_object=True) ** n * p
    print("DM dt/cycle=%5.3fs" % ((timeit.default_timer() - t0) / n))
    p = AP(update_object=True) ** 10 * p
    t0 = timeit.default_timer()
    p = AP(update_object=True) ** n * p
    print("AP dt/cycle=%5.3fs" % ((timeit.default_timer() - t0) / n))
    p = ML(update_object=True) ** 10 * p
    t0 = timeit.default_timer()
    p = ML(update_object=True) ** n * p
    print("ML dt/cycle=%5.3fs" % ((timeit.default_timer() - t0) / n))

if False:
    # Timing vs stack size
    n = 50
    vx = []
    vy = []
    for stack_size in range(1, 32 + 1, 1):
        default_processing_unit.cl_stack_size = np.int32(stack_size)
        p = DM(update_object=True) ** 10 * p
        t0 = timeit.default_timer()
        p = DM(update_object=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('Timd for a DM cycle (s)')