.. _concepts: ******************** Concepts and Classes ******************** .. _ptypyclasses: Tutorial: Data access - Storage, View, Container ================================================ .. note:: This tutorial was generated from the python source :file:`[ptypy_root]/tutorial/ptypyclasses.py` using :file:`ptypy/doc/script2rst.py`. You are encouraged to modify the parameters and rerun the tutorial with:: $ python [ptypy_root]/tutorial/ptypyclasses.py This tutorial explains how |ptypy| accesses data / memory. The data access is governed by three classes: :py:class:`~ptypy.core.classes.Container`, :py:class:`~ptypy.core.classes.Storage`, and :py:class:`~ptypy.core.classes.View` First a short reminder from the :py:mod:`~ptypy.core.classes` Module about the classes this tutorial covers. .. parsed-literal:: **Container class** A high-level container that keeps track of sub-containers (Storage) and Views onto them. A container can copy itself to produce a buffer needed for calculations. Some basic Mathematical operations are implemented at this level as in-place operations. In general, operations on arrays should be done using the Views, which simply return numpy arrays. **Storage class** The sub-container, wrapping a numpy array buffer. A Storage defines a system of coordinate (for now only a scaled translation of the pixel coordinates, but more complicated affine transformation could be implemented if needed). The sub-class DynamicStorage can adapt the size of its buffer (cropping and/or padding) depending on the Views. **View class** A low-weight class that contains all information to access a 2D piece of a Storage within a Container. The basic idea is that the View access is controlled by a physical position and its frame, such that one is not bothered by memory/array addresses when accessing data. Import some modules :: >>> import matplotlib as mpl >>> import numpy as np >>> import ptypy >>> from ptypy import utils as u >>> from ptypy.core import View, Container, Storage, Base >>> plt = mpl.pyplot A single Storage in one Container --------------------------------- The governing object is a :any:`Container` instance, which we need to construct as first object. It does not need much for an input but the data type. :: >>> C1 = Container(data_type='real') This class itself holds nothing at first. In order to store data in ``C1`` we have to add a :any:`Storage` to that container. :: >>> S1 = C1.new_storage(shape=(1, 7, 7)) Similarly we can contruct the Storage from a data buffer. ``S1=C1.new_storage(data=np.ones((1,7,7)))`` would have also worked. All of |ptypy|'s special classes carry a uniue ID, which is needed to communicate these classes across nodes and for saving and loading. As we haven't specified an ID the Container class picks one for ``S1`` In this case that will be ``S0000`` where the *S* refers to the class type. :: >>> print(S1.ID) S0000 Let's have a look at what kind of data Storage holds. :: >>> print(S1.formatted_report()[0]) (C)ontnr : Memory : Shape : Pixel size : Dimensions : Views (S)torgs : (MB) : (Pixel) : (meters) : (meters) : act. -------------------------------------------------------------------------------- S0000 : 0.0 : 1 * 7 * 7 : 1.0*1.0e0 : 7.0*7.0e0 : 0 Apart from the ID on the left we discover a few other entries, for example the quantity ``psize`` which refers to the physical pixel size for the last two dimensions of the stored data. :: >>> print(S1.psize) [1. 1.] Many attributes of a :any:`Storage` are in fact *properties*. Changing their value may have an impact on other methods or attributes of the class. For example, one convenient method is Storage.\ :py:meth:`~ptypy.core.classes.Storage.grids` that creates coordinate grids for the last two dimensions (see also :py:func:`ptypy.utils.array_utils.grids`) :: >>> y, x = S1.grids() >>> print(y) [[[-3. -3. -3. -3. -3. -3. -3.] [-2. -2. -2. -2. -2. -2. -2.] [-1. -1. -1. -1. -1. -1. -1.] [ 0. 0. 0. 0. 0. 0. 0.] [ 1. 1. 1. 1. 1. 1. 1.] [ 2. 2. 2. 2. 2. 2. 2.] [ 3. 3. 3. 3. 3. 3. 3.]]] >>> print(x) [[[-3. -2. -1. 0. 1. 2. 3.] [-3. -2. -1. 0. 1. 2. 3.] [-3. -2. -1. 0. 1. 2. 3.] [-3. -2. -1. 0. 1. 2. 3.] [-3. -2. -1. 0. 1. 2. 3.] [-3. -2. -1. 0. 1. 2. 3.] [-3. -2. -1. 0. 1. 2. 3.]]] which are coordinate grids for vertical and horizontal axes respectively. We note that these coordinates have their center at :: >>> print(S1.center) [3 3] Now we change a few properties. For example, :: >>> S1.center = (2, 2) >>> S1.psize = 0.1 >>> g = S1.grids() >>> print(g[0]) [[[-0.2 -0.2 -0.2 -0.2 -0.2 -0.2 -0.2] [-0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1] [ 0. 0. 0. 0. 0. 0. 0. ] [ 0.1 0.1 0.1 0.1 0.1 0.1 0.1] [ 0.2 0.2 0.2 0.2 0.2 0.2 0.2] [ 0.3 0.3 0.3 0.3 0.3 0.3 0.3] [ 0.4 0.4 0.4 0.4 0.4 0.4 0.4]]] >>> print(g[1]) [[[-0.2 -0.1 0. 0.1 0.2 0.3 0.4] [-0.2 -0.1 0. 0.1 0.2 0.3 0.4] [-0.2 -0.1 0. 0.1 0.2 0.3 0.4] [-0.2 -0.1 0. 0.1 0.2 0.3 0.4] [-0.2 -0.1 0. 0.1 0.2 0.3 0.4] [-0.2 -0.1 0. 0.1 0.2 0.3 0.4] [-0.2 -0.1 0. 0.1 0.2 0.3 0.4]]] We observe that the center has in fact moved one pixel up and one left. The :py:func:`~ptypy.core.classes.Storage.center` property uses pixel units. If we want to use a physical quantity to shift the center, we may instead set the top left pixel to a new value, which shifts the center to a new position. :: >>> S1.origin -= 0.12 >>> y, x = S1.grids() >>> print(y) [[[-0.32 -0.32 -0.32 -0.32 -0.32 -0.32 -0.32] [-0.22 -0.22 -0.22 -0.22 -0.22 -0.22 -0.22] [-0.12 -0.12 -0.12 -0.12 -0.12 -0.12 -0.12] [-0.02 -0.02 -0.02 -0.02 -0.02 -0.02 -0.02] [ 0.08 0.08 0.08 0.08 0.08 0.08 0.08] [ 0.18 0.18 0.18 0.18 0.18 0.18 0.18] [ 0.28 0.28 0.28 0.28 0.28 0.28 0.28]]] >>> print(x) [[[-0.32 -0.22 -0.12 -0.02 0.08 0.18 0.28] [-0.32 -0.22 -0.12 -0.02 0.08 0.18 0.28] [-0.32 -0.22 -0.12 -0.02 0.08 0.18 0.28] [-0.32 -0.22 -0.12 -0.02 0.08 0.18 0.28] [-0.32 -0.22 -0.12 -0.02 0.08 0.18 0.28] [-0.32 -0.22 -0.12 -0.02 0.08 0.18 0.28] [-0.32 -0.22 -0.12 -0.02 0.08 0.18 0.28]]] >>> print(S1.center) [3.2 3.2] Up until now our actual *data* numpy array located at ``S1.data`` is still filled with ones, i.e. flat. We can use :any:`Storage.fill` to fill that container with an array. :: >>> S1.fill(x+y) >>> print(S1.data) [[[-0.64 -0.54 -0.44 -0.34 -0.24 -0.14 -0.04] [-0.54 -0.44 -0.34 -0.24 -0.14 -0.04 0.06] [-0.44 -0.34 -0.24 -0.14 -0.04 0.06 0.16] [-0.34 -0.24 -0.14 -0.04 0.06 0.16 0.26] [-0.24 -0.14 -0.04 0.06 0.16 0.26 0.36] [-0.14 -0.04 0.06 0.16 0.26 0.36 0.46] [-0.04 0.06 0.16 0.26 0.36 0.46 0.56]]] We can have visual check on the data using :py:func:`~ptypy.utils.plot_utils.plot_storage` :: >>> fig = u.plot_storage(S1, 0) See :numref:`ptypyclasses_00` for the plotted image. .. figure:: ../_script2rst/ptypyclasses_00.png :width: 72 % :figclass: highlights :name: ptypyclasses_00 A plot of the data stored in ``S1`` Adding Views as a way to access data ------------------------------------ Besides being able to access the data directly through its attribute and the corresponding *numpy* syntax, ptypy offers acces through a :any:`View` instance. The View invocation is a bit more complex. :: >>> from ptypy.core.classes import DEFAULT_ACCESSRULE >>> ar = DEFAULT_ACCESSRULE.copy() >>> print(ar) * id3V1AJ3O73G : ptypy.utils.parameters.Param(6) * storageID : None * shape : None * coord : None * psize : 1.0 * layer : 0 * active : True Let's say we want a 4x4 view on Storage ``S1`` around the origin. We set :: >>> ar.shape = (4, 4) # ar.shape = 4 would have been also valid. >>> ar.coord = 0. # ar.coord = (0.,0.) would have been accepted, too. >>> ar.storageID = S1.ID >>> ar.psize = None Now we can construct the View. Upon construction, the View will access information from Storage ``S1`` to initialize other attributes for access and geometry. :: >>> V1 = View(C1, ID=None, accessrule=ar) We see that a number of the accessrule items appear in the View now. :: >>> print(V1.shape) [4 4] >>> print(V1.coord) [0. 0.] >>> print(V1.storageID) S0000 A few others were set by the automatic update of Storage ``S1``. :: >>> print(V1.psize) [0.1 0.1] >>> print(V1.storage) S0000 : 0.00 MB :: data=(1, 7, 7) @float64 psize=[0.1 0.1] center=[3.2 3.2] The update also set new attributes of the View which all start with a lower ``d`` and are locally cached information about data access. :: >>> print(V1.dlayer, V1.dlow, V1.dhigh) 0 [1 1] [5 5] Finally, we can retrieve the data subset by applying the View to the storage. :: >>> data = S1[V1] >>> print(data) [[-0.44 -0.34 -0.24 -0.14] [-0.34 -0.24 -0.14 -0.04] [-0.24 -0.14 -0.04 0.06] [-0.14 -0.04 0.06 0.16]] It does not matter if we apply the View to Storage ``S1`` or the container ``C1``, or use the View internal View.\ :py:meth:`~ptypy.core.classes.View.data` property. :: >>> print(np.allclose(data, C1[V1])) True >>> print(np.allclose(data, V1.data)) True The first access yielded a similar result because the :py:attr:`~ptypy.core.classes.View.storageID` ``S0000`` is in ``C1`` and the second acces method worked because it uses the View's :py:attr:`~ptypy.core.classes.View.storage` attribute. :: >>> print(V1.storage is S1) True >>> print(V1.storageID in C1.S.keys()) True We observe that the coordinate [0.0,0.0] is not part of the grid in S1 anymore. Consequently, the View was put as close to [0.0,0.0] as possible. The coordinate in data space that the View would have as center is the attribute :py:meth:`~ptypy.core.classes.View.pcoord` while :py:meth:`~ptypy.core.classes.View.dcoord` is the closest data coordinate. The difference is held by :py:meth:`~ptypy.core.classes.View.sp` such that a subpixel correction may be applied if needed (future release) :: >>> print(V1.dcoord, V1.pcoord, V1.sp) [3 3] [3.2 3.2] [0.2 0.2] .. note:: Please note that we cannot guarantee any API stability for other attributes / properties besides *.data*, *.shape* and *.coord* If we set the coordinate to some other value in the grid, we can eliminate the subpixel misfit. By changing the *.coord* property, we need to update the View manually, as the View-Storage interaction is non-automatic apart from the moment the View is constructed - a measure to prevent unwanted feedback loops. :: >>> V1.coord = (0.08, 0.08) >>> S1.update_views(V1) >>> print(V1.dcoord, V1.pcoord, V1.sp) [4 4] [4. 4.] [0. 0.] We observe that the high range limit of the View is close to the border of the data buffer. :: >>> print(V1.dhigh) [6 6] What happens if we push the coordinate further? :: >>> V1.coord = (0.28, 0.28) >>> S1.update_views(V1) >>> print(V1.dhigh) [8 8] Now the higher range limit of the View is off bounds for sure. Applying this View to the Storage may lead to undesired behavior, i.e. array concatenation or data access errors. :: >>> print(S1[V1]) [[0.16 0.26 0.36] [0.26 0.36 0.46] [0.36 0.46 0.56]] >>> print(S1[V1].shape, V1.shape) (3, 3) [4 4] One important feature of the :any:`Storage` class is that it can detect all out-of-bounds accesses and reformat the data buffer accordingly. A simple call to *Storage*.\ :py:meth:`~ptypy.core.classes.Storage.reformat` should do. :: >>> print(S1.shape) (1, 7, 7) >>> mn = S1[V1].mean() >>> S1.fill_value = mn >>> S1.reformat() >>> print(S1.shape) (1, 4, 4) Oh no, the Storage data buffer has shrunk! But don't worry, that is intended behavior. A call to *.reformat()* crops and pads the data buffer around all **active** Views. We would need to set ``S1.padonly = True`` if we wanted to avoid that the data buffer is cropped. We leave this as an exercise for now. Instead, we add a new View at a different location to verify that the buffer will adapt to reach both Views. :: >>> ar2 = ar.copy() >>> ar2.coord = (-0.82, -0.82) >>> V2 = View(C1, ID=None, accessrule=ar2) >>> S1.fill_value = 0. >>> S1.reformat() >>> print(S1.shape) (1, 15, 15) Ok, we note that the buffer has grown in size. Now, we give the new View some copied values of the other view to make the View appear in a plot. :: >>> V2.data = V1.data.copy() >>> fig = u.plot_storage(S1, 2) See :numref:`ptypyclasses_02` for the plotted image. .. figure:: ../_script2rst/ptypyclasses_02.png :width: 72 % :figclass: highlights :name: ptypyclasses_02 Storage with 4x4 views of the same content. We observe that the data buffer spans both views. Now let us add more and more Views! :: >>> for i in range(1, 11): >>> ar2 = ar.copy() >>> ar2.coord = (-0.82+i*0.1, -0.82+i*0.1) >>> View(C1, ID=None, accessrule=ar2) >>> A handy method of the :any:`Storage` class is to determine the *coverage*, which maps, for every pixel of the storage, the number of views having access to this pixel. :: >>> S1.data[:] = S1.get_view_coverage() >>> fig = u.plot_storage(S1, 3) See :numref:`ptypyclasses_03` for the plotted image. .. figure:: ../_script2rst/ptypyclasses_03.png :width: 72 % :figclass: highlights :name: ptypyclasses_03 View coverage of data buffer of ``S1``. Another handy feature of the :any:`View` class is that it automatically creates a Storage instance to a ``storageID`` if that ID does not already exist. :: >>> ar = DEFAULT_ACCESSRULE.copy() >>> ar.shape = 200 >>> ar.coord = 0. >>> ar.storageID = 'S100' >>> ar.psize = 1.0 >>> V3=View(C1, ID=None, accessrule=ar) Finally we have a look at the mischief we managed so far. :: >>> print(C1.formatted_report()) (C)ontnr : Memory : Shape : Pixel size : Dimensions : Views (S)torgs : (MB) : (Pixel) : (meters) : (meters) : act. -------------------------------------------------------------------------------- None : 0.3 : float64 S0000 : 0.0 : 1 * 15 * 15 : 1.0*1.0e-1 : 1.5*1.5e0 : 12 S100 : 0.3 : 1 * 200 * 200 : 1.0*1.0e0 : 2.0*2.0e2 : 1 .. _simupod: Tutorial: Modeling the experiment - Pod, Geometry ================================================= .. note:: This tutorial was generated from the python source :file:`[ptypy_root]/tutorial/simupod.py` using :file:`ptypy/doc/script2rst.py`. You are encouraged to modify the parameters and rerun the tutorial with:: $ python [ptypy_root]/tutorial/simupod.py In the :ref:`ptypyclasses` we have learned to deal with the basic storage-and-access classes on small toy arrays. In this tutorial we will learn how to create :any:`POD` and :any:`Geo` instances to imitate a ptychography experiment and to use larger arrays. We would like to point out that the "data" created her is not actual data. There is neither light or other wave-like particle involved, nor actual diffraction. You will also not find an actual movement of motors or stages, nor is there an actual detector. Everything should be understood as a test for this software. The selected physical quantities only simulate an experimental setup. We start again with importing some modules. :: >>> import matplotlib as mpl >>> import numpy as np >>> import ptypy >>> from ptypy import utils as u >>> from ptypy.core import View, Container, Storage, Base, POD >>> plt = mpl.pyplot >>> import sys >>> scriptname = sys.argv[0].split('.')[0] We create a managing top-level class instance. We will not use the the :any:`Ptycho` class for now, as its rich set of methods may be a bit overwhelming to start with. Instead we start of a plain :py:class:`~ptypy.core.classes.Base` instance. As the ``Base`` class is slotted starting with version 0.3.0, we can't attach attributes after initialisation as we could with a normal python class. To reenable ths feature we have to subclasse once. :: >>> class Base2(Base): >>> pass >>> :: >>> P = Base2() Now we can start attaching attributes :: >>> P.CType = np.complex128 >>> P.FType = np.float64 Set experimental setup geometry and create propagator ----------------------------------------------------- Here, we accept a little help from the :any:`Geo` class to provide a propagator and pixel sizes for sample and detector space. :: >>> from ptypy.core import geometry >>> g = u.Param() >>> g.energy = None # u.keV2m(1.0)/6.32e-7 >>> g.lam = 5.32e-7 >>> g.distance = 15e-2 >>> g.psize = 24e-6 >>> g.shape = 256 >>> g.propagation = "farfield" >>> G = geometry.Geo(owner=P, pars=g) The Geo instance ``G`` has done a lot already at this moment. For example, we find forward and backward propagators at ``G.propagator.fw`` and ``G.propagator.bw``. It has also calculated the appropriate pixel size in the sample plane (aka resolution), :: >>> print(G.resolution) [1.29882812e-05 1.29882812e-05] which sets the shifting frame to be of the following size: :: >>> fsize = G.shape * G.resolution >>> print("%.2fx%.2fmm" % tuple(fsize*1e3)) 3.32x3.32mm Create probing illumination --------------------------- Next, we need to create a probing illumination. We start with a suited container that we call *probe* :: >>> P.probe = Container(P, 'Cprobe', data_type='complex') For convenience, there is a test probing illumination in |ptypy|'s resources. :: >>> from ptypy.resources import moon_pr >>> moon_probe = -moon_pr(G.shape) As of version 0.3.0 the Storage classes are no longer limited to 3d arrays. This means, however, that the shape is better explicitly specified upon creation to align the coordinate system. :: >>> pr_shape = (1,)+ tuple(G.shape) >>> pr = P.probe.new_storage(shape=pr_shape, psize=G.resolution) Probe is loaded through the :py:meth:`~ptypy.core.classes.Storage.fill` method. :: >>> pr.fill(moon_probe) >>> fig = u.plot_storage(pr, 0) See :numref:`simupod_00` for the plotted image. .. figure:: ../_script2rst/simupod_00.png :width: 72 % :figclass: highlights :name: simupod_00 Ptypy's default testing illumination, an image of the moon. Of course, we could have also used the coordinate grids from the propagator to model a probe, :: >>> y, x = G.propagator.grids_sam >>> apert = u.smooth_step(fsize[0]/5-np.sqrt(x**2+y**2), 1e-6) >>> pr2 = P.probe.new_storage(shape=pr_shape, psize=G.resolution) >>> pr2.fill(apert) >>> fig = u.plot_storage(pr2, 1, channel='c') See :numref:`simupod_01` for the plotted image. .. figure:: ../_script2rst/simupod_01.png :width: 72 % :figclass: highlights :name: simupod_01 Round test illumination (not necessarily smoothly varying). or the coordinate grids from the Storage itself. :: >>> pr3 = P.probe.new_storage(shape=pr_shape, psize=G.resolution) >>> y, x = pr3.grids() >>> apert = u.smooth_step(fsize[0]/5-np.abs(x), 3e-5)*u.smooth_step(fsize[1]/5-np.abs(y), 3e-5) >>> pr3.fill(apert) >>> fig = u.plot_storage(pr3, 2) See :numref:`simupod_02` for the plotted image. .. figure:: ../_script2rst/simupod_02.png :width: 72 % :figclass: highlights :name: simupod_02 Square test illumination. In order to put some physics in the illumination we set the number of photons to 1 billion :: >>> for pp in [pr, pr2, pr3]: >>> pp.data *= np.sqrt(1e9/np.sum(pp.data*pp.data.conj())) :: >>> print(u.norm2(pr.data)) 1000000000.0000001 and we quickly check if the propagation works. :: >>> ill = pr.data[0] >>> propagated_ill = G.propagator.fw(ill) >>> fig = plt.figure(3) >>> ax = fig.add_subplot(111) >>> im = ax.imshow(np.log10(np.abs(propagated_ill)+1)) >>> plt.colorbar(im) See :numref:`simupod_03` for the plotted image. .. figure:: ../_script2rst/simupod_03.png :width: 80 % :figclass: highlights :name: simupod_03 Logarithmic intensity of propagated illumination. Create scan pattern and object ------------------------------ We use the :py:mod:`ptypy.core.xy` module to create a scan pattern. :: >>> pos = u.Param() >>> pos.model = "round" >>> pos.spacing = fsize[0]/8 >>> pos.steps = None >>> pos.extent = fsize*1.5 >>> from ptypy.core import xy >>> positions = xy.from_pars(pos) >>> fig = plt.figure(4) >>> ax = fig.add_subplot(111) >>> ax.plot(positions[:, 1], positions[:, 0], 'o-') See :numref:`simupod_04` for the plotted image. .. figure:: ../_script2rst/simupod_04.png :width: 80 % :figclass: highlights :name: simupod_04 Created scan pattern. Next, we need to model a sample through an object transmisson function. We start of with a suited container which we call *obj*. :: >>> P.obj = Container(P, 'Cobj', data_type='complex') As we have learned from the previous :ref:`ptypyclasses`\ , we can use :any:`View`\ 's to create a Storage data buffer of the right size. :: >>> oar = View.DEFAULT_ACCESSRULE.copy() >>> oar.storageID = 'S00' >>> oar.psize = G.resolution >>> oar.layer = 0 >>> oar.shape = G.shape >>> oar.active = True :: >>> for pos in positions: >>> # the rule >>> r = oar.copy() >>> r.coord = pos >>> V = View(P.obj, None, r) >>> We let the Storages in ``P.obj`` reformat in order to include all Views. Conveniently, this can be initiated from the top with Container.\ :py:meth:`~ptypy.core.classes.Container.reformat` :: >>> P.obj.reformat() >>> print(P.obj.formatted_report()) (C)ontnr : Memory : Shape : Pixel size : Dimensions : Views (S)torgs : (MB) : (Pixel) : (meters) : (meters) : act. -------------------------------------------------------------------------------- Cobj : 6.5 : complex128 S00 : 6.5 : 1 * 638 * 640 : 1.3*1.3e-5 : 8.3*8.3e-3 : 116 At last we fill the object Storage ``S00`` with a complex transmission. Again there is a convenience transmission function in the resources :: >>> from ptypy.resources import flower_obj >>> storage = P.obj.storages['S00'] >>> storage.fill(flower_obj(storage.shape[-2:])) >>> fig = u.plot_storage(storage, 5) See :numref:`simupod_05` for the plotted image. .. figure:: ../_script2rst/simupod_05.png :width: 72 % :figclass: highlights :name: simupod_05 Complex transmission function of the object. Creating additional Views and the PODs -------------------------------------- A single coherent propagation in |ptypy| is represented by an instance of the :py:class:`~ptypy.core.classes.POD` class. :: >>> print(POD.__doc__) POD : Ptychographic Object Descriptor A POD brings together probe view, object view and diff view. It also gives access to "exit", a (coherent) exit wave, and to propagation objects to go from exit to diff space. >>> print(POD.__init__.__doc__) Parameters ---------- ptycho : Ptycho The instance of Ptycho associated with this pod. model : ~ptypy.core.manager.ScanModel The instance of ScanModel (or it subclasses) which describes this pod. ID : str or int The pod ID, If None it is managed by the ptycho. views : dict or Param The views. See :py:attr:`DEFAULT_VIEWS`. geometry : Geo Geometry class instance and attached propagator For creating a single POD we need a :py:class:`~ptypy.core.classes.View` to *probe*, *object*, *exit* wave and *diff*\ raction containers as well as the :any:`Geo` class instance which represents the experimental setup. First we create the missing :py:class:`~ptypy.core.classes.Container`'s. :: >>> P.exit = Container(P, 'Cexit', data_type='complex') >>> P.diff = Container(P, 'Cdiff', data_type='real') >>> P.mask = Container(P, 'Cmask', data_type='real') We start with one POD and its views. :: >>> objviews = list(P.obj.views.values()) >>> obview = objviews[0] We construct the probe View. :: >>> probe_ar = View.DEFAULT_ACCESSRULE.copy() >>> probe_ar.psize = G.resolution >>> probe_ar.shape = G.shape >>> probe_ar.active = True >>> probe_ar.storageID = pr.ID >>> prview = View(P.probe, None, probe_ar) We construct the exit wave View. This construction is shorter as we only change a few bits in the access rule. :: >>> exit_ar = probe_ar.copy() >>> exit_ar.layer = 0 >>> exit_ar.active = True >>> exview = View(P.exit, None, exit_ar) We construct diffraction and mask Views. Even shorter is the construction of the mask View as, for the mask, we are essentially using the same access as for the diffraction data. :: >>> diff_ar = probe_ar.copy() >>> diff_ar.layer = 0 >>> diff_ar.active = True >>> diff_ar.psize = G.psize >>> mask_ar = diff_ar.copy() >>> maview = View(P.mask, None, mask_ar) >>> diview = View(P.diff, None, diff_ar) Now we can create the POD. :: >>> pods = [] >>> views = {'probe': prview, 'obj': obview, 'exit': exview, 'diff': diview, 'mask': maview} >>> pod = POD(P, ID=None, views=views, geometry=G) >>> pods.append(pod) The :any:`POD` is the most important class in |ptypy|. Its instances are used to write the reconstruction algorithms using its attributes as local references. For example we can create and store an exit wave in the following convenient fashion: :: >>> pod.exit = pod.probe * pod.object The result of the calculation above is stored in the appropriate storage of ``P.exit``. Therefore we can use this command to plot the result. :: >>> exit_storage = list(P.exit.storages.values())[0] >>> fig = u.plot_storage(exit_storage, 6) See :numref:`simupod_06` for the plotted image. .. figure:: ../_script2rst/simupod_06.png :width: 72 % :figclass: highlights :name: simupod_06 Simulated exit wave using a pod. The diffraction plane is also conveniently accessible with :: >>> pod.diff = np.abs(pod.fw(pod.exit))**2 The result is stored in the diffraction container. :: >>> diff_storage = list(P.diff.storages.values())[0] >>> fig = u.plot_storage(diff_storage, 7, modulus='log') See :numref:`simupod_07` for the plotted image. .. figure:: ../_script2rst/simupod_07.png :width: 72 % :figclass: highlights :name: simupod_07 First simulated diffraction image (without noise) Creating the rest of the pods is now straight-forward since the data accesses are similar. :: >>> for obview in objviews[1:]: >>> # we keep the same probe access >>> prview = View(P.probe, None, probe_ar) >>> # For diffraction and exit wave we need to increase the >>> # layer index as there is a new exit wave and diffraction pattern for each >>> # scan position >>> exit_ar.layer += 1 >>> diff_ar.layer += 1 >>> exview = View(P.exit, None, exit_ar) >>> maview = View(P.mask, None, mask_ar) >>> diview = View(P.diff, None, diff_ar) >>> views = {'probe': prview, 'obj': obview, 'exit': exview, 'diff': diview, 'mask': maview} >>> pod = POD(P, ID=None, views=views, geometry=G) >>> pods.append(pod) >>> We let the storage arrays adapt to the new Views. :: >>> for C in [P.mask, P.exit, P.diff, P.probe]: >>> C.reformat() >>> And the rest of the simulation fits in three lines of code! :: >>> for pod in pods: >>> pod.exit = pod.probe * pod.object >>> # we use Poisson statistics for a tiny bit of realism in the >>> # diffraction images >>> pod.diff = np.random.poisson(np.abs(pod.fw(pod.exit))**2) >>> pod.mask = np.ones_like(pod.diff) >>> >>> We make a quick check on the diffraction patterns :: >>> fig = u.plot_storage(diff_storage, 8, slices=':2,:,:', modulus='log') See :numref:`simupod_08` for the plotted image. .. figure:: ../_script2rst/simupod_08.png :width: 100 % :figclass: highlights :name: simupod_08 Diffraction patterns with poisson statistics. **Well done!** We can move forward to create and run a reconstruction engine as in section :ref:`basic_algorithm` in :ref:`ownengine` or store the generated diffraction patterns as in the next section. .. _store: Storing the simulation ---------------------- On unix system we choose the ``/tmp`` folder :: >>> save_path = '/tmp/ptypy/sim/' >>> import os :: >>> if not os.path.exists(save_path): >>> os.makedirs(save_path) >>> First, we save the geometric info in a text file. :: >>> with open(save_path+'geometry.txt', 'w') as f: >>> f.write('distance %.4e\n' % G.p.distance) >>> f.write('energy %.4e\n' % G.energy) >>> f.write('psize %.4e\n' % G.psize[0]) >>> f.write('shape %d\n' % G.shape[0]) >>> f.close() >>> Next, we save positions and the diffraction images. We don't burden ourselves for now with converting to an image file format such as .tiff or a data format like .hdf5 but instead we use numpy's binary storage format. :: >>> with open(save_path+'positions.txt', 'w') as f: >>> if not os.path.exists(save_path+'ccd/'): >>> os.mkdir(save_path+'ccd/') >>> for pod in pods: >>> diff_frame = 'ccd/diffraction_%04d.npy' % pod.di_view.layer >>> f.write(diff_frame+' %.4e %.4e\n' % tuple(pod.ob_view.coord)) >>> frame = pod.diff.astype(np.int32) >>> np.save(save_path+diff_frame, frame) >>> If you want to learn how to convert this "experiment" into ptypy data file (``.ptyd``), see to :ref:`subclassptyscan` .. _ownengine: Tutorial: A reconstruction engine from scratch ============================================== .. note:: This tutorial was generated from the python source :file:`[ptypy_root]/tutorial/ownengine.py` using :file:`ptypy/doc/script2rst.py`. You are encouraged to modify the parameters and rerun the tutorial with:: $ python [ptypy_root]/tutorial/ownengine.py In this tutorial, we want to provide the information needed to create an engine compatible with the state mixture expansion of ptychogrpahy as described in Thibault et. al 2013 [#modes]_ . First we import ptypy and the utility module :: >>> import ptypy >>> from ptypy import utils as u >>> import numpy as np Preparing a managing Ptycho instance ------------------------------------ We need to prepare a managing :any:`Ptycho` instance. It requires a parameter tree, as specified in :ref:`parameters` First, we create a most basic input paramater tree. There are many default values, but we specify manually only a more verbose output and single precision for the data type. :: >>> p = u.Param() >>> p.verbose_level = 3 >>> p.data_type = "single" Now, we need to create a set of scans that we wish to reconstruct in the reconstruction run. We will use a single scan that we call 'MF' and mark the data source as 'test' to use the |ptypy| internal :py:class:`MoonFlowerScan` :: >>> p.scans = u.Param() >>> p.scans.MF = u.Param() >>> p.scans.MF.name = 'Full' >>> p.scans.MF.data = u.Param() >>> p.scans.MF.data.name = 'MoonFlowerScan' >>> p.scans.MF.data.shape = 128 >>> p.scans.MF.data.num_frames = 400 This bare parameter tree will be the input for the :any:`Ptycho` class which is constructed at ``level=2``. It means that it creates all necessary basic :any:`Container` instances like *probe*, *object* *diff* , etc. It also loads the first chunk of data and creates all :any:`View` and :any:`POD` instances, as the verbose output will tell. :: >>> P = ptypy.core.Ptycho(p, level=2) A quick look at the diffraction data :: >>> diff_storage = list(P.diff.storages.values())[0] >>> fig = u.plot_storage(diff_storage, 0, slices=(slice(2), slice(None), slice(None)), modulus='log') See :numref:`ownengine_00` for the plotted image. .. figure:: ../_script2rst/ownengine_00.png :width: 100 % :figclass: highlights :name: ownengine_00 Plot of simulated diffraction data for the first two positions. We don't need to use |ptypy|'s :any:`Ptycho` class to arrive at this point. The structure ``P`` that we arrive with at the end of :ref:`simupod` suffices completely. Probe and object are not so exciting to look at for now. As default, probes are initialized with an aperture like support. :: >>> probe_storage = list(P.probe.storages.values())[0] >>> fig = u.plot_storage(list(P.probe.S.values())[0], 1) See :numref:`ownengine_01` for the plotted image. .. figure:: ../_script2rst/ownengine_01.png :width: 72 % :figclass: highlights :name: ownengine_01 Plot of the starting guess for the probe. .. _basic_algorithm: A basic Difference-Map implementation ------------------------------------- Now we can start implementing a simple DM algorithm. We need three basic functions, one is the ``fourier_update`` that implements the Fourier modulus constraint. .. math:: \psi_{d,\lambda,k} = \mathcal{D}_{\lambda,z}^{-1}\left\{\sqrt{I_{d}}\frac{\mathcal{D}_{\lambda,z} \{\psi_{d,\lambda,k}\}}{\sum_{\lambda,k} |\mathcal{D}_{\lambda,z} \{\psi_{d,\lambda,k}\} |^2}\right\} :: >>> def fourier_update(pods): >>> import numpy as np >>> pod = list(pods.values())[0] >>> # Get Magnitude and Mask >>> mask = pod.mask >>> modulus = np.sqrt(np.abs(pod.diff)) >>> # Create temporary buffers >>> Imodel = np.zeros_like(pod.diff) >>> err = 0. >>> Dphi = {} >>> # Propagate the exit waves >>> for gamma, pod in pods.items(): >>> Dphi[gamma] = pod.fw(2*pod.probe*pod.object - pod.exit) >>> Imodel += np.abs(Dphi[gamma] * Dphi[gamma].conj()) >>> # Calculate common correction factor >>> factor = (1-mask) + mask * modulus / (np.sqrt(Imodel) + 1e-10) >>> # Apply correction and propagate back >>> for gamma, pod in pods.items(): >>> df = pod.bw(factor*Dphi[gamma]) - pod.probe*pod.object >>> pod.exit += df >>> err += np.mean(np.abs(df*df.conj())) >>> # Return difference map error on exit waves >>> return err >>> :: >>> def probe_update(probe, norm, pods, fill=0.): >>> """ >>> Updates `probe`. A portion `fill` of the probe is kept from >>> iteration to iteration. Requires `norm` buffer and pod dictionary >>> """ >>> probe *= fill >>> norm << fill + 1e-10 >>> for name, pod in pods.items(): >>> if not pod.active: continue >>> probe[pod.pr_view] += pod.object.conj() * pod.exit >>> norm[pod.pr_view] += pod.object * pod.object.conj() >>> # For parallel usage (MPI) we have to communicate the buffer arrays >>> probe.allreduce() >>> norm.allreduce() >>> probe /= norm >>> :: >>> def object_update(obj, norm, pods, fill=0.): >>> """ >>> Updates `object`. A portion `fill` of the object is kept from >>> iteration to iteration. Requires `norm` buffer and pod dictionary >>> """ >>> obj *= fill >>> norm << fill + 1e-10 >>> for pod in pods.values(): >>> if not pod.active: continue >>> pod.object += pod.probe.conj() * pod.exit >>> norm[pod.ob_view] += pod.probe * pod.probe.conj() >>> obj.allreduce() >>> norm.allreduce() >>> obj /= norm >>> :: >>> def iterate(Ptycho, num): >>> # generate container copies >>> obj_norm = P.obj.copy(fill=0.) >>> probe_norm = P.probe.copy(fill=0.) >>> errors = [] >>> for i in range(num): >>> err = 0 >>> # fourier update >>> for di_view in Ptycho.diff.V.values(): >>> if not di_view.active: continue >>> err += fourier_update(di_view.pods) >>> # probe update >>> probe_update(Ptycho.probe, probe_norm, Ptycho.pods) >>> # object update >>> object_update(Ptycho.obj, obj_norm, Ptycho.pods) >>> # print error >>> errors.append(err) >>> if i % 3==0: print(err) >>> # cleanup >>> P.obj.delete_copy() >>> P.probe.delete_copy() >>> #return error >>> return errors >>> We start off with a small number of iterations. :: >>> iterate(P, 9) 438026.50735607144 315508.66172933887 237978.26703255248 We note that the error (here only displayed for 3 iterations) is already declining. That is a good sign. Let us have a look how the probe has developed. :: >>> fig = u.plot_storage(list(P.probe.S.values())[0], 2) See :numref:`ownengine_02` for the plotted image. .. figure:: ../_script2rst/ownengine_02.png :width: 72 % :figclass: highlights :name: ownengine_02 Plot of the reconstructed probe after 9 iterations. We observe that the actaul illumination of the sample must be larger than the initial guess. Looks like the probe is on a good way. How about the object? :: >>> fig = u.plot_storage(list(P.obj.S.values())[0], 3, slices='0,120:-120,120:-120') See :numref:`ownengine_03` for the plotted image. .. figure:: ../_script2rst/ownengine_03.png :width: 72 % :figclass: highlights :name: ownengine_03 Plot of the reconstructed object after 9 iterations. It is not quite clear what object is reconstructed Ok, let us do some more iterations. 36 will do. :: >>> iterate(P, 36) 158603.9740395377 91072.73033678501 61211.753241738574 46684.904276998306 37483.80987198136 28430.613557863566 22350.477334251234 19980.413638760696 18558.9012019789 17776.334601063198 17781.86195633913 17806.30507956264 Error is still on a steady descent. Let us look at the final reconstructed probe and object. :: >>> fig = u.plot_storage(list(P.probe.S.values())[0], 4) See :numref:`ownengine_04` for the plotted image. .. figure:: ../_script2rst/ownengine_04.png :width: 72 % :figclass: highlights :name: ownengine_04 Plot of the reconstructed probe after a total of 45 iterations. It's a moon ! >>> fig = u.plot_storage(list(P.obj.S.values())[0], 5, slices='0,120:-120,120:-120') See :numref:`ownengine_05` for the plotted image. .. figure:: ../_script2rst/ownengine_05.png :width: 72 % :figclass: highlights :name: ownengine_05 Plot of the reconstructed object after a total of 45 iterations. It's a bunch of flowers ! .. [#modes] P. Thibault and A. Menzel, **Nature** 494, 68 (2013)