diff --git a/.travis.yml b/.travis.yml index 5e18e4424..fac0049c0 100644 --- a/.travis.yml +++ b/.travis.yml @@ -2,8 +2,7 @@ cache: apt sudo: true language: python python: - - "2.7" - - "3.3" + - 2.7 before_install: - wget https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh -O miniconda.sh; # grab miniconda - bash miniconda.sh -b -p $HOME/miniconda # install miniconda @@ -13,12 +12,21 @@ before_install: - conda install pyyaml - conda info -a # and print the info -install: - - conda env create --file environment.yml #set up the environment thats in the file - - source activate test-environment # activate it - - python setup.py install # install ptypy +env: + - TEST_ENV_NAME=core_dependencies + - TEST_ENV_NAME=full_dependencies + +install: True script: + - conda env create --file ${TEST_ENV_NAME}.yml #set up the environment thats in the file + - source activate ${TEST_ENV_NAME} # activate it + - conda install pytest # additional dependencies for the tests + - pip install pytest-cov + - pip install coveralls + - echo $PYTHONPATH + - conda list + - python setup.py install # install ptypy - py.test ptypy/test -v --cov ptypy --cov-report term-missing # now run the tests after_script: diff --git a/benchmark/model_speed.py b/benchmark/model_speed.py new file mode 100644 index 000000000..1983abc77 --- /dev/null +++ b/benchmark/model_speed.py @@ -0,0 +1,16 @@ + +from ptypy import utils as u +from ptypy.core import Ptycho, Vanilla + +p = u.Param() +p.verbose_level = 3 +p.io = u.Param() +p.io.home = "/tmp/ptypy/" +p.scans = u.Param() +p.scans.MF = u.Param() +p.scans.MF.data= u.Param() +p.scans.MF.name = 'Vanilla' +p.scans.MF.data.name = 'QuickScan' +p.scans.MF.data.num_frames = 50000 +p.scans.MF.data.shape = 32 +Ptycho(p,level=2) diff --git a/core_dependencies.yml b/core_dependencies.yml new file mode 100644 index 000000000..957f3c8b0 --- /dev/null +++ b/core_dependencies.yml @@ -0,0 +1,8 @@ +name: core_dependencies +channels: + - conda-forge +dependencies: + - python=2.7 + - numpy + - scipy + - h5py diff --git a/doc/Makefile b/doc/Makefile index decf44118..9de9237e2 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -3,8 +3,13 @@ # You can set these variables from the command line. SPHINXOPTS = -SPHINXBUILD = ~/.local/bin/sphinx-build + +ifeq ($(SPHINXBUILD),"") + SPHINXBUILD = ~/.local/bin/sphinx-build +endif + #SPHINXBUILD = /usr/bin/sphinx-build + PAPER = BUILDDIR = build diff --git a/doc/conf.py b/doc/conf.py index 422944800..5a7ec654e 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -21,11 +21,11 @@ #sys.path.insert(0, os.path.abspath('.')) # generate paramters.rst and other rst -execfile('parameters2rst.py') -execfile('tmp2rst.py') -execfile('version.py') import subprocess -subprocess.call(['python', 'script2rst.py']) # We need this to have a clean sys.argv +subprocess.check_call(['python', 'script2rst.py']) # We need this to have a clean sys.argv +subprocess.check_call(['python','parameters2rst.py']) +subprocess.check_call(['python','tmp2rst.py']) +execfile('version.py') # -- General configuration ------------------------------------------------ @@ -42,73 +42,63 @@ ] +def truncate_docstring(app, what, name, obj, options, lines): + """ + Remove the Default parameter entries. + """ + if not hasattr(obj, 'DEFAULT'): + return + if any(l.strip().startswith('Defaults:') for l in lines): + while True: + if lines.pop(-1).strip().startswith('Defaults:'): + break + + def remove_mod_docstring(app, what, name, obj, options, lines): from ptypy import utils as u - import numpy as np + from ptypy import defaults_tree u.verbose.report.headernewline='\n\n' searchstr = ':py:data:' - def get_refs(dct,pd,depth=2, indent = ''): - if depth<0: + def get_refs(dct, pd, depth=2, indent=''): + if depth < 0: return for k, value in dct.iteritems(): ref = ', see :py:data:`~%s`' % pd.children[k].entry_point if pd.children.has_key(k) else '' - if hasattr(value,'items'): + if hasattr(value, 'items'): v = str(value.__class__.__name__) - elif str(value)==value: - v='"%s"' % value + elif str(value) == value: + v = '"%s"' % value else: - v=str(value) + v = str(value) - lines.append(indent+'* *' +k+'* = ``'+v+'``' +ref)#+'\n') + lines.append(indent + '* *' + k + '* = ``' + v + '``' + ref) - if hasattr(value,'items'): - #lines.append('\n\n') + if hasattr(value, 'items'): lines.append("") - get_refs(value,pd.children[k],depth=depth-1, indent = indent+' ') + get_refs(value, pd.children[k], depth=depth-1, indent=indent+' ') lines.append("") - #lines.append('\n\n') - #if name.find('DEFAULT')>=0: if isinstance(obj, u.Param) or isinstance(obj, dict): - keys = obj.keys() pd = None - """ - # auto_matching - for entry,pdesc in u.validator.entry_points_Param.iteritems(): - chkeys = ':'.join([k.split('.')[-1] for k in pdesc.children.keys()]) - #print chkeys - #print keys - matches = [key in chkeys for key in keys] - #print matches - print np.mean(matches) - if np.mean(matches)>0.8: - print 'Param match' - e=entry - print e - pd = pdesc - break - """ for l in lines: start = l.find(searchstr) if start > -1: newstr = l[start:] newstr = newstr.split('`')[1] newstr = newstr.replace('~', '') - #print newstr, what, name, options - pd = u.validator.entry_points_dct.get(newstr,None) + pd = defaults_tree.get(newstr) break if pd is not None: - #lines.append('Match with :py:data:`.%s` \n\n' %pd.entry_point) get_refs(obj, pd, depth=2, indent='') - #print lines def setup(app): app.connect('autodoc-process-docstring', remove_mod_docstring) + app.connect('autodoc-process-docstring', truncate_docstring) napoleon_use_ivar = True diff --git a/doc/html_templates/ptypysphinx/download.html b/doc/html_templates/ptypysphinx/download.html index c9a67a666..59a1f3960 100644 --- a/doc/html_templates/ptypysphinx/download.html +++ b/doc/html_templates/ptypysphinx/download.html @@ -32,5 +32,5 @@ diff --git a/doc/index.rst b/doc/index.rst index a59d037ed..bc755bc01 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -44,7 +44,7 @@ Highlights overcoming partial coherence or related phenomena. * **On-the-fly** reconstructions (while data is being acquired) using the - the :any:`PtyScan` class in the linking mode :ref:`linking mode` + the :py:class:`PtyScan` class in the linking mode :ref:`linking mode` Quicklinks @@ -68,10 +68,10 @@ Quicklinks .. rubric:: Footnotes -.. [#Enders2016] B.Enders and P.Thibault, **Proc. R. Soc. A** 472, 20160640 (2016), `doi `_ +.. [#Enders2016] B.Enders and P.Thibault, **Proc. R. Soc. A** 472, 20160640 (2016), `doi `__ -.. [#Thi2013] P.Thibault and A.Menzel, **Nature** 494, 68 (2013), `doi `_ +.. [#Thi2013] P.Thibault and A.Menzel, **Nature** 494, 68 (2013), `doi `__ -.. [#ml] P.Thibault and M.Guizar-Sicairos, **New J. of Phys.** 14, 6 (2012), `doi `_ +.. [#ml] P.Thibault and M.Guizar-Sicairos, **New J. of Phys.** 14, 6 (2012), `doi `__ -.. [#dm] P.Thibault, M.Dierolf *et al.*, **New J. of Phys. 14**, 6 (2012), `doi `_ +.. [#dm] P.Thibault, M.Dierolf *et al.*, **New J. of Phys. 14**, 6 (2012), `doi `__ diff --git a/doc/parameters2rst.py b/doc/parameters2rst.py index 419e3537f..37f93563a 100644 --- a/doc/parameters2rst.py +++ b/doc/parameters2rst.py @@ -1,5 +1,5 @@ -from ptypy import utils as u +from ptypy import defaults_tree prst = open('rst/parameters.rst','w') @@ -9,32 +9,52 @@ Header+= '************************\n\n' prst.write(Header) -names = u.validator.parameter_descriptions -for name,desc in u.validator.parameter_descriptions.iteritems(): - if name=='': +for path, desc in defaults_tree.descendants: + + types = desc.type + default = desc.default + lowlim, uplim = desc.limits + is_wildcard = (desc.name == '*') + + if is_wildcard: + path = path.replace('*', desc.parent.name[:-1] + '_00') + + if path == '': continue - if hasattr(desc,'children') and desc.parent is u.validator.pdroot: - prst.write('\n'+name+'\n') - prst.write('='*len(name)+'\n\n') - if hasattr(desc,'children') and desc.parent.parent is u.validator.pdroot: - prst.write('\n'+name+'\n') - prst.write('-'*len(name)+'\n\n') - - prst.write('.. py:data:: '+name) - prst.write('('+', '.join([t for t in desc.type])+')') + if desc.children and desc.parent is desc.root: + prst.write('\n'+path+'\n') + prst.write('='*len(path)+'\n\n') + if desc.children and desc.parent.parent is desc.root: + prst.write('\n'+path+'\n') + prst.write('-'*len(path)+'\n\n') + + prst.write('.. py:data:: '+path) + + if desc.is_symlink: + tp = 'Param' + else: + tp = ', '.join([str(t) for t in types]) + prst.write(' ('+tp+')') prst.write('\n\n') - num = str(desc.ID) if hasattr(desc,'ID') else "None" - prst.write(' *('+num+')* '+desc.shortdoc+'\n\n') - prst.write(' '+desc.longdoc.replace('\n','\n ')+'\n\n') - prst.write(' *default* = ``'+str(desc.default)) - if desc.lowlim is not None and desc.uplim is not None: - prst.write(' (>'+str(desc.lowlim)+', <'+str(desc.uplim)+')``\n') - elif desc.lowlim is not None and desc.uplim is None: - prst.write(' (>'+str(desc.lowlim)+')``\n') - elif desc.lowlim is None and desc.uplim is not None: - prst.write(' (<'+str(desc.uplim)+')``\n') + + if is_wildcard: + prst.write(' *Wildcard*: multiple entries with arbitrary names are accepted.\n\n') + + prst.write(' '+desc.help+'\n\n') + prst.write(' '+desc.doc.replace('','\n').replace('\n', '\n ')+'\n\n') + + if desc.is_symlink: + prst.write(' *default* = '+':py:data:`'+desc.path+'`\n') else: - prst.write('``\n') + prst.write(' *default* = ``'+repr(default)) + if lowlim is not None and uplim is not None: + prst.write(' (>'+str(lowlim)+', <'+str(uplim)+')``\n') + elif lowlim is not None and uplim is None: + prst.write(' (>'+str(lowlim)+')``\n') + elif lowlim is None and uplim is not None: + prst.write(' (<'+str(uplim)+')``\n') + else: + prst.write('``\n') prst.write('\n') prst.close() diff --git a/doc/rst/ptypy.core.rst b/doc/rst/ptypy.core.rst index b2fa0cba0..93837fdad 100644 --- a/doc/rst/ptypy.core.rst +++ b/doc/rst/ptypy.core.rst @@ -44,14 +44,6 @@ ptypy.core.manager module :undoc-members: :show-inheritance: -ptypy.core.model module ------------------------ - -.. automodule:: ptypy.core.model - :members: - :undoc-members: - :show-inheritance: - ptypy.core.paths module ----------------------- diff --git a/doc/rst/ptypy.debug.rst b/doc/rst/ptypy.debug.rst new file mode 100644 index 000000000..edc921576 --- /dev/null +++ b/doc/rst/ptypy.debug.rst @@ -0,0 +1,29 @@ +ptypy.debug package +=================== + +Submodules +---------- + +ptypy.debug.embedded_shell module +--------------------------------- + +.. automodule:: ptypy.debug.embedded_shell + :members: + :undoc-members: + :show-inheritance: + +ptypy.debug.ipython_kernel module +--------------------------------- + +.. automodule:: ptypy.debug.ipython_kernel + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: ptypy.debug + :members: + :undoc-members: + :show-inheritance: diff --git a/doc/rst/ptypy.rst b/doc/rst/ptypy.rst index d462b9cf7..46e7c1267 100644 --- a/doc/rst/ptypy.rst +++ b/doc/rst/ptypy.rst @@ -13,6 +13,7 @@ Subpackages ptypy.resources ptypy.simulations ptypy.utils + ptypy.debug Module contents --------------- diff --git a/doc/rst/ptypy.utils.rst b/doc/rst/ptypy.utils.rst index c4bda94ab..fb59aeea1 100644 --- a/doc/rst/ptypy.utils.rst +++ b/doc/rst/ptypy.utils.rst @@ -12,14 +12,6 @@ ptypy.utils.array_utils module :undoc-members: :show-inheritance: -ptypy.utils.embedded_shell module ---------------------------------- - -.. automodule:: ptypy.utils.embedded_shell - :members: - :undoc-members: - :show-inheritance: - ptypy.utils.math_utils module ----------------------------- @@ -76,10 +68,10 @@ ptypy.utils.scripts module :undoc-members: :show-inheritance: -ptypy.utils.validator module ----------------------------- +ptypy.utils.descriptor module +----------------------------- -.. automodule:: ptypy.utils.validator +.. automodule:: ptypy.utils.descriptor :members: :undoc-members: :show-inheritance: diff --git a/doc/rst_templates/data_management.tmp b/doc/rst_templates/data_management.tmp index fc15caaff..cb3fe7b4a 100644 --- a/doc/rst_templates/data_management.tmp +++ b/doc/rst_templates/data_management.tmp @@ -42,9 +42,9 @@ prior to a successful image reconstruction. reconstruction software naturally **cannot** provide actual experimental data. Nevertheless, the treatment of raw data is usually very similar for every experiment. Consequently, |ptypy| provides an abstract base class, -called :any:`PtyScan`, which aims to help with steps (B) and (C). In order +called :py:class:`PtyScan`, which aims to help with steps (B) and (C). In order to adapt |ptypy| for a specific experimental setup, we simply -subclass :any:`PtyScan` and reimplement only that subset of its methods which are +subclass :py:class:`PtyScan` and reimplement only that subset of its methods which are affected by the specifics of the experiemental setup (see :ref:`subclassptyscan`). @@ -53,7 +53,7 @@ affected by the specifics of the experiemental setup The PtyScan class ================= -:any:`PtyScan` is the abstract base class in |ptypy| that manages raw input +:py:class:`PtyScan` is the abstract base class in |ptypy| that manages raw input data. A PtyScan instance is constructed from a set of generic parameters, @@ -67,7 +67,7 @@ It provides the following features: process only loads the data it will later use in the reconstruction. Hence, the load on the network is not affected by the number of processes. - The parallel behavior of :any:`PtyScan`, is controlled by the parameter + The parallel behavior of :py:class:`PtyScan`, is controlled by the parameter :py:data:`.scan.data.load_parallel`. It uses the :py:class:`~ptypy.utils.parallel.LoadManager` **Preparation** @@ -135,7 +135,7 @@ The PtyScan class of |ptypy| provides support for three use cases. In this use case, the researcher has integrated |ptypy| into the beamline end-station or experimental setup - with the help of a custom subclass of :any:`PtyScan` that we call + with the help of a custom subclass of :py:class:`PtyScan` that we call ``UserScan``. This subclass has its own methods to extract many of the of the generic parameters of :py:data:`.scan.data` and also defaults for specific custom parameters, for instance file paths or file name @@ -149,7 +149,7 @@ The PtyScan class of |ptypy| provides support for three use cases. :figclass: highlights :name: case_integrated - Intergrated use case of :any:`PtyScan`. + Integrated use case of :py:class:`PtyScan`. A custom subclass ``UserScan`` serves as a translator between |ptypy|'s generic parameters and @@ -189,11 +189,11 @@ The PtyScan class of |ptypy| provides support for three use cases. :figclass: highlights :name: case_prepared - Standard supported use case of :any:`PtyScan`. + Standard supported use case of :py:class:`PtyScan`. If a structure-compatible (see :ref:`ptyd_file`) ``*.hdf5``-file is available, |ptypy| can be used without customizing a subclass of - :any:`PtyScan`. It will use the shipped subclass :any:`PtydScan` + :py:class:`PtyScan`. It will use the shipped subclass :py:class:`PtydScan` to read in the (prepared) raw data. **Preparation and reconstruction on-the-fly with data acquisition.** @@ -218,7 +218,7 @@ The PtyScan class of |ptypy| provides support for three use cases. :figclass: highlights :name: case_flyscan - On-the-fly or demon-like use case of :any:`PtyScan`. + On-the-fly or demon-like use case of :py:class:`PtyScan`. A separate process prepares the data *chunks* and saves them in separate files which are @@ -299,7 +299,7 @@ are general (optional) parameters. expected in the dataset (see :py:data:`~.scan.data.num_frames`) It is important to set this parameter when the data acquisition is not finished but the reconstruction has already started. If the dataset - is complete, the loading class :any:`PtydScan` retrieves the + is complete, the loading class :py:class:`PtydScan` retrieves the total number of frames from the payload ``chunks/`` The next set of optional parameters are diff --git a/doc/rst_templates/getting_started.tmp b/doc/rst_templates/getting_started.tmp index cdcfcdccd..0b5d7850e 100644 --- a/doc/rst_templates/getting_started.tmp +++ b/doc/rst_templates/getting_started.tmp @@ -378,4 +378,3 @@ branch :py:data:`.engine` for detailed parameter descriptions. .. literalinclude:: ../../templates/minimal_load_and_run.py :language: python :linenos: - :emphasize-lines: 25-56 diff --git a/doc/script2rst.py b/doc/script2rst.py index 352679301..0a19c45ac 100644 --- a/doc/script2rst.py +++ b/doc/script2rst.py @@ -12,12 +12,15 @@ if len(sys.argv) == 1: import pkg_resources + import subprocess + for script in scripts: scr = pkg_resources.resource_filename('ptypy', tutorial_dir+script) if not os.path.exists(scr): print('Using backup tutorial for %s' % script) scr = '../tutorial/'+script - os.system('python '+sys.argv[0]+' '+scr) + #subprocess.call(['python',sys.argv[0]+' '+scr]) # doesn't work + os.system('python ' + sys.argv[0]+' '+scr) sys.exit() indent_keys = ['for', 'if', 'with', 'def', 'class'] @@ -126,12 +129,17 @@ def check_for_fig(wline): frst.write(' ' + line2) continue + decorator = False indent = False for key in indent_keys: if line.startswith(key): indent = True break - + + if line.startswith('@'): + indent = True + decorator = True + if indent: frst.write('\n::\n\n >>> '+line) func = line @@ -139,9 +147,12 @@ def check_for_fig(wline): while True: line2 = fpy.readline() if line2.strip() and not line2.startswith(' '): - frst.write('\n') - fpy.seek(pt) - break + if decorator: + decorator = False + else: + frst.write('\n') + fpy.seek(pt) + break func += line2 frst.write(' >>> '+line2) pt = fpy.tell() diff --git a/doc/version.py b/doc/version.py index 6fa1f4553..e36958bf6 100644 --- a/doc/version.py +++ b/doc/version.py @@ -1,8 +1,8 @@ # THIS FILE IS GENERATED FROM ptypy/setup.py -short_version='0.2.1' -version='0.2.1' -release=True +short_version='0.3.0' +version='0.3.0' +release=False if not release: version += '.dev' diff --git a/extra/ipynb/plotclient.ipynb b/extra/ipynb/plotclient.ipynb index 88133f5ce..e9767d994 100644 --- a/extra/ipynb/plotclient.ipynb +++ b/extra/ipynb/plotclient.ipynb @@ -11,7 +11,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false, "extensions": { "jupyter_dashboards": { "version": 1, @@ -30,7 +29,8 @@ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import time\n", - "from ptypy.utils import PlotClient, validator, ortho\n", + "from ptypy.utils import PlotClient, ortho\n", + "from ptypy.utils.descriptor import defaults_tree\n", "from ptypy.utils.plot_utils import PtyAxis\n", "from ptypy.io import h5read" ] @@ -46,7 +46,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false, "extensions": { "jupyter_dashboards": { "version": 1, @@ -61,7 +60,8 @@ }, "outputs": [], "source": [ - "pc = PlotClient(validator.make_sub_default('.io.interaction'),in_thread=False)\n", + "#pc = PlotClient(defaults_tree['io.interaction'].make_default(depth=5), in_thread=False)\n", + "pc = PlotClient(in_thread=False)\n", "data = lambda: pc.get_data()\n", "pc.start()" ] @@ -76,9 +76,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "pc.stop()" @@ -94,9 +92,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "base_path = '...'\n", @@ -115,9 +111,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "ptya = PtyAxis(channel='a',cmap='viridis', vmin=None, vmax=None)\n", @@ -136,9 +130,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "ptya = PtyAxis(channel='c',cmap='magma', vmax=None)\n", @@ -157,9 +149,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "fig = plt.figure()\n", @@ -216,7 +206,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", - "version": "2.7.10" + "version": "2.7.12" } }, "nbformat": 4, diff --git a/environment.yml b/full_dependencies.yml similarity index 80% rename from environment.yml rename to full_dependencies.yml index 6c48b7ddd..6e7ab2da2 100644 --- a/environment.yml +++ b/full_dependencies.yml @@ -1,4 +1,4 @@ -name: test-environment +name: full_dependencies channels: - conda-forge dependencies: @@ -6,14 +6,13 @@ dependencies: - numpy - scipy - matplotlib - - pytest - h5py - pyzmq - pep8 - - pytest - mpi4py - pil - pyfftw - pip: - pytest-cov - coveralls + - fabio diff --git a/ptypy/__init__.py b/ptypy/__init__.py index 0c59f29fb..d6cee956f 100644 --- a/ptypy/__init__.py +++ b/ptypy/__init__.py @@ -1,73 +1,79 @@ #!/usr/bin/env python2 # -*- coding: utf-8 -*- -from version import short_version, version, release +from .version import short_version, version, release + +__doc__ = \ """ PTYPY(v%(short)s): A ptychography reconstruction package. To cite PTYPY in publications, use @article{ ... } - :copyright: Copyright 2014 by the PTYPY team, see AUTHORS. + :copyright: Copyright 2018 by the PTYPY team, see AUTHORS. :license: GPLv2, see LICENSE for details. :version: %(version)s :status: %(devrel)s -""" % {'version':version,'devrel':'release' if release else 'development','short':short_version} +""" % {'version': version, 'devrel': 'release' if release else 'development', 'short': short_version} del short_version, release try: import zmq - __has_zmq__= True +except ImportError as ie: + __has_zmq__ = False +else: + __has_zmq__ = True del zmq -except ImportError('ZeroMQ not found.\nInteraction server & client disabled.\n\ -Install python-zmq via the package repositories or with `pip install --user pyzmq`'): - __has_zmq__= False try: import mpi4py - __has_mpi4py__= True +except ImportError as ie: + __has_mpi4py__ = False +else: + __has_mpi4py__ = True del mpi4py -except ImportError('Message Passaging for Python (mpi4py) not found.\n\ -CPU-parallelization disabled.\n\ -Install python-mpi4py via the package repositories or with `pip install --user mpi4py`'): - __has_mpi4py__= False try: import matplotlib - __has_matplotlib__= True +except ImportError as ie: + __has_matplotlib__ = False +else: + __has_matplotlib__ = True del matplotlib -except ImportError('Plotting for Python (matplotlib) not found.\n\ -Plotting disabled.\n\ -Install python-matplotlib via the package repositories or with `pip install --user matplotlib`'): - __has_matplotlib__= False # Initialize MPI (eventually GPU) -from utils import parallel +from .utils import parallel # Logging -from utils import verbose -#verbose.set_level(2) - -import utils -# Import core modules -import io -import experiment -import core +from .utils import verbose +# Log immediately dependency information +if not __has_zmq__: + __zmq_msg = 'ZeroMQ not found.\nInteraction server & client disabled.\n\ + Install python-zmq via the package repositories or with `pip install --user pyzmq`' + verbose.logger.warn(__zmq_msg) +if not __has_mpi4py__: + __mpi_msg = 'Message Passaging for Python (mpi4py) not found.\n\ + CPU-parallelization disabled.\n\ + Install python-mpi4py via the package repositories or with `pip install --user mpi4py`' + verbose.logger.warn(__mpi_msg) +if not __has_matplotlib__: + __mpl_msg = 'Plotting for Python (matplotlib) not found.\n\ + Plotting disabled.\n\ + Install python-matplotlib via the package repositories or with `pip install --user matplotlib`' + verbose.logger.warn(__mpl_msg) -#from core import * -import simulations -import resources +# Start a parameter tree +from .utils.descriptor import EvalDescriptor +defaults_tree = EvalDescriptor('root') +del EvalDescriptor -if __name__ == "__main__": - # TODO: parse command line arguments for extra options - import sys - # Get parameters from command line argument - param_filename = sys.argv[1] - p = parameters.load(param_filename) - # Initialize Ptycho object - pt = Ptycho(p) +# Import core modules +from . import utils +from . import io +from . import experiment +from . import core +from . import simulations +from . import resources - # Start reconstruction - pt.run() diff --git a/ptypy/core/__init__.py b/ptypy/core/__init__.py index ed7352e35..38d47f3d5 100644 --- a/ptypy/core/__init__.py +++ b/ptypy/core/__init__.py @@ -7,12 +7,13 @@ :copyright: Copyright 2014 by the PTYPY team, see AUTHORS. :license: GPLv2, see LICENSE for details. """ -from classes import * -from ptycho import * -from manager import * -import xy -import illumination -import sample -import geometry -import save_load -import data +from .classes import * +from .ptycho import * +from .manager import * +from . import xy +from . import illumination +from . import sample +from . import geometry +from . import geometry_bragg +from . import save_load +from . import data diff --git a/ptypy/core/classes.py b/ptypy/core/classes.py index 051ebfb5a..3e33fecdd 100644 --- a/ptypy/core/classes.py +++ b/ptypy/core/classes.py @@ -34,6 +34,8 @@ """ import numpy as np import weakref +from collections import OrderedDict + try: from pympler.asizeof import asizeof use_asizeof = True @@ -86,7 +88,10 @@ class Base(object): _CHILD_PREFIX = 'ID' _PREFIX = BASE_PREFIX - + + __slots__ = ['ID','numID','owner','_pool','_recs','_record'] + _fields = [('ID','= l: + nl = l + 8192 if idx > 10000 else 2*l + recs = np.resize(recs,(nl,)) + self._recs[prefix] = recs + rec = recs[idx] + obj._record = rec + rec['ID'] = nID + return - + @staticmethod def _num_to_id(num): """ @@ -179,12 +198,18 @@ def _num_to_id(num): @classmethod def _from_dict(cls, dct): """ - create new instance from dictionary dct + Create new instance from dictionary dct should be compatible with _to_dict() """ inst = cls.__new__(cls) - inst.__dict__.update(dct) - + for k in inst.__slots__: + if k not in dct: + continue + else: + setattr(inst,k ,dct[k]) + if hasattr(inst,'__dict__'): + inst.__dict__.update(dct) + # Calling post dictionary import routine (empty in base) inst._post_dict_import() @@ -203,7 +228,12 @@ def _to_dict(self): a dictionary. Overwrite in child for custom behavior. Default. Returns shallow copy of internal dict as default """ - return self.__dict__.copy() + res = OrderedDict() + for k in self.__slots__: + res[k] = getattr(self, k) + if hasattr(self, '__dict__'): + res.update(self.__dict__.copy()) + return res def calc_mem_usage(self): space = 64 # that is for the class itself @@ -217,22 +247,24 @@ def calc_mem_usage(self): space += asizeof(v, limit=0) for kk, vv in v.iteritems(): pool_space += vv.calc_mem_usage()[0] - - for k, v in self.__dict__.iteritems(): - if issubclass(type(v), Base): - continue - elif str(k) == '_pool' or str(k) == 'pods': - continue - else: - if use_asizeof: - s = asizeof(v) - space += s - if type(v) is np.ndarray: - npy_space += v.nbytes + + if hasattr(self, '__dict__'): + for k, v in self.__dict__.iteritems(): + if issubclass(type(v), Base): + continue + elif str(k) == '_pool' or str(k) == 'pods': + continue + else: + if use_asizeof: + s = asizeof(v) + space += s + if type(v) is np.ndarray: + npy_space += v.nbytes return space + pool_space + npy_space, pool_space, npy_space + def get_class(ID): """ Determine ptypy class from unique `ID` @@ -290,8 +322,8 @@ class Storage(Base): _PREFIX = STORAGE_PREFIX - def __init__(self, container, ID=None, data=None, shape=(1, 1, 1), fill=0., - psize=None, origin=None, layermap=None, padonly=False, + def __init__(self, container, ID=None, data=None, shape=DEFAULT_SHAPE, + fill=0., psize=1., origin=None, layermap=None, padonly=False, **kwargs): """ Parameters @@ -305,9 +337,11 @@ def __init__(self, container, ID=None, data=None, shape=(1, 1, 1), fill=0., data: ndarray or None A numpy array to use as initial buffer. + *Deprecated* in v0.3, use fill method instead. - shape : 3-tuple - The shape of the buffer + shape : tuple, int, or None + The shape of the buffer. If None or int, the dimensionality + is found from the owning Container instance. fill : float or complex The default value to fill storage with, will be converted to @@ -338,19 +372,28 @@ def __init__(self, container, ID=None, data=None, shape=(1, 1, 1), fill=0., self.fill_value = fill if fill is not None else 0. # For documentation - #: Three dimensional array as data buffer + #: Three/four or potentially N-dimensional array as data buffer self.data = None + # dimensionality suggestion from container + ndim = container.ndim if container.ndim is not None else 2 + if shape is None: - shape = DEFAULT_SHAPE + shape = (1,) + (1,) * ndim elif np.isscalar(shape): - shape = (1, int(shape), int(shape)) - elif len(shape) == 2: - shape = (1,) + tuple(shape) - elif len(shape) != 3: - raise ValueError( - '`shape` must be None or scalar or 2-tuple or 3-tuple of int') + shape = (1,) + (int(shape),) * ndim + else: + shape = tuple(shape) + if len(shape) not in [3, 4]: + logger.warn('Storage view access dimension %d is not in regular ' + 'scope (2,3). Behavior is untested.' % len(shape[1:])) + + self.shape = shape + self.data = np.empty(self.shape, self.dtype) + self.data.fill(self.fill_value) + + """ # Set data buffer if data is None: # Create data buffer @@ -360,15 +403,17 @@ def __init__(self, container, ID=None, data=None, shape=(1, 1, 1), fill=0., else: # Set initial buffer. Casting the type makes a copy data = np.asarray(data).astype(self.dtype) - if data.ndim < 2 or data.ndim > 3: + if data.ndim < self.ndim or data.ndim > (self.ndim + 1): raise ValueError( - 'Initial buffer must be 2D or 3D, this one is %dD.' - % data.ndim) - elif data.ndim == 2: + 'For `data_dims` = %d, initial buffer must be' + ' %dD or %dD, this one is %dD' + % (self.ndim, self.ndim, self.ndim+1, data.ndim)) + elif data.ndim == self.ndim: self.data = data.reshape((1,) + data.shape) else: self.data = data self.shape = self.data.shape + """ if layermap is None: layermap = range(len(self.data)) @@ -379,7 +424,7 @@ def __init__(self, container, ID=None, data=None, shape=(1, 1, 1), fill=0., # Need to bootstrap the parameters. # We set the initial center to the middle of the array - self._center = u.expect2(self.shape[-2:]) // 2 + self._center = u.expectN(self.shape[-self.ndim:], self.ndim) // 2 # Set pixel size (in physical units) self.psize = psize if psize is not None else DEFAULT_PSIZE @@ -406,40 +451,12 @@ def __init__(self, container, ID=None, data=None, shape=(1, 1, 1), fill=0., # solution required # self._origin = None - def _to_dict(self): + @property + def ndim(self): """ - We will have to recompute the datalist here + Number of dimensions for :any:`View` access """ - cp = self.__dict__.copy() - # Delete datalist reference - try: - del cp['_datalist'] - except: - pass - - return cp - # self._make_datalist() - - def _make_datalist(self): - pass - """ - # BE does not give the same result on all nodes - #self._datalist = [None] * (max(self.layermap)+1) - #u.parallel.barrier() - #print u.parallel.rank - self._datalist = [None] * max(self.nlayers,max(self.layermap)+1) - for k,i in enumerate(self.layermap): - self._datalist[i] = self.data[k] - """ - """ - @property - def datalist(self): - - if not hasattr(self,'datalist'): - self._make_datalist() - - return self._datalist - """ + return len(self.shape[1:]) @property def dtype(self): @@ -459,24 +476,28 @@ def copy(self, owner=None, ID=None, fill=None): If float, set the content to this value. If None, copy the current content. """ - if fill is None: - # Return a new Storage or sub-class object with a copy of the data. - return type(self)(owner, - ID, - data=self.data.copy(), - psize=self.psize, - origin=self.origin, - layermap=self.layermap) - else: - # Return a new Storage or sub-class object with an empty buffer - new_storage = type(self)(owner, - ID, - shape=self.shape, - psize=self.psize, - origin=self.origin, - layermap=self.layermap) + + # if fill is None: + # # Return a new Storage or sub-class object with a copy of the data. + # return type(self)(owner, + # ID, + # data=self.data.copy(), + # psize=self.psize, + # origin=self.origin, + # layermap=self.layermap) + # else: + # # Return a new Storage or sub-class object with an empty buffer + new_storage = type(self)(owner, + ID, + shape=self.shape, + psize=self.psize, + origin=self.origin, + layermap=self.layermap) + if fill is not None: new_storage.fill(fill) - return new_storage + else: + new_storage.fill(self.data.copy()) + return new_storage def fill(self, fill=None): """ @@ -501,12 +522,13 @@ def fill(self, fill=None): self.data.fill(fill) elif type(fill) is np.ndarray: # Replace the buffer - if fill.ndim < 2 or fill.ndim > 3: - raise ValueError('Numpy ndarray fill must be 2D or 3D.') - elif fill.ndim == 2: - fill = np.resize(fill, (self.shape[0], - fill.shape[0], - fill.shape[1])) + if fill.ndim < self.ndim or fill.ndim > (self.ndim + 1): + raise ValueError( + 'For `data_dims` = %d, initial buffer must be' + ' %dD or %dD, this one is %dD' + % (self.ndim, self.ndim, self.ndim+1, fill.ndim)) + elif fill.ndim == self.ndim: + fill = np.resize(fill, (self.shape[0],) + fill.shape) self.data = fill.astype(self.dtype) self.shape = self.data.shape @@ -535,33 +557,32 @@ def update_views(self, v=None): self.update_views(v) return + if not self.ndim == v.ndim: + raise ValueError( + 'Storage %s(ndim=%d) and View %s(ndim=%d) have conflicting ' + 'data dimensions' % (self.ID, self.ndim, v.ID, v.ndim)) + # Synchronize pixel size - v.psize = self.psize.copy() + v.psize = self.psize #.copy() - # v.shape can be None upon initialization - this means "full frame" - if v.shape is None: - v.shape = u.expect2(self.shape[-2:]) - v.pcoord = v.shape / 2. - v.coord = self._to_phys(v.pcoord) - else: - # Convert the physical coordinates of the view to pixel coordinates - v.pcoord = self._to_pix(v.coord) + # Convert the physical coordinates of the view to pixel coordinates + pcoord = self._to_pix(v.coord) # Integer part (note that np.round is not stable for odd arrays) - v.dcoord = np.round(v.pcoord + 0.00001).astype(int) + v.dcoord = np.round(pcoord + 0.00001).astype(int) # These are the important attributes used when accessing the data - v.dlow = v.dcoord - v.shape/2 - v.dhigh = v.dcoord + (v.shape + 1)/2 + v.dlow = v.dcoord - v.shape / 2 + v.dhigh = v.dcoord + (v.shape + 1) / 2 # v.roi = np.array([pix - v.shape/2, pix + (v.shape + 1)/2]) - v.sp = v.pcoord - v.dcoord + v.sp = pcoord - v.dcoord # if self.layermap is None: # v.slayer = 0 # else: # v.slayer = self.layermap.index(v.layer) - def reformat(self, newID=None): + def reformat(self, newID=None, update=True): """ Crop or pad if required. @@ -570,6 +591,11 @@ def reformat(self, newID=None): newID : str or int If None (default) act on self. Otherwise create a copy of self before doing the cropping and padding. + + update : bool + If True, updates all Views before reformatting. Not necessarily + needed, if Views have been recently instantiated. Roughly doubles + execution time. Returns ------- @@ -580,10 +606,11 @@ def reformat(self, newID=None): # If a new storage is requested, make a copy. if newID is not None: s = self.copy(newID) - s.reformat(newID=None) + s.reformat() return s # Make sure all views are up to date + # This call takes roughly half the time of .reformat() self.update() # List of views on this storage @@ -594,34 +621,39 @@ def reformat(self, newID=None): logger.debug('%s[%s] :: %d views for this storage' % (self.owner.ID, self.ID, len(views))) + sh = self.data.shape + # Loop through all active views to get individual boundaries - rows = [] - cols = [] + #axes = [[]] * self.ndim + + mn = [np.inf] * self.ndim + mx = [-np.inf] * self.ndim layers = [] + dims = range(self.ndim) for v in views: if not v.active: continue # Accumulate the regions of interest to # compute the full field of view - # rows += [v.roi[0, 0], v.roi[1, 0]] - # cols += [v.roi[0, 1], v.roi[1, 1]] - rows += [v.dlow[0], v.dhigh[0]] - cols += [v.dlow[1], v.dhigh[1]] - + for d in dims: + #axes[d] += [v.dlow[d], v.dhigh[d]] + mn[d] = min(mn[d],v.dlow[d]) + mx[d] = max(mx[d],v.dhigh[d]) + # Gather a (unique) list of layers if v.layer not in layers: layers.append(v.layer) sh = self.data.shape - # Compute 2d misfit (distance between the buffer boundaries and the + # Compute Nd misfit (distance between the buffer boundaries and the # region required to fit all the views) - misfit = np.array([[-np.min(rows), np.max(rows) - sh[-2]], - [-np.min(cols), np.max(cols) - sh[-1]]]) - - logger.debug('%s[%s] :: misfit = [%s, %s]' - % (self.owner.ID, self.ID, misfit[0], misfit[1])) + misfit = np.array([[-mn[d], mx[d] - sh[d+1]] for d in dims]) + + _misfit_str = ', '.join(['%s' % m for m in misfit]) + logger.debug('%s[%s] :: misfit = [%s]' + % (self.owner.ID, self.ID, _misfit_str)) posmisfit = (misfit > 0) negmisfit = (misfit < 0) @@ -631,19 +663,20 @@ def reformat(self, newID=None): if posmisfit.any() or negmisfit.any(): logger.debug( - 'Storage %s of container %s has a misfit of [%s, %s] between ' + 'Storage %s of container %s has a misfit of [%s] between ' 'its data and its views' - % (str(self.ID), str(self.owner.ID), misfit[0], misfit[1])) - + % (str(self.ID), str(self.owner.ID), _misfit_str)) + if needtocrop_or_pad: if self.padonly: misfit[negmisfit] = 0 # Recompute center and shape new_center = self.center + misfit[:, 0] - new_shape = (sh[0], - sh[1] + misfit[0].sum(), - sh[2] + misfit[1].sum()) + new_shape = (sh[0],) + for d in dims: + new_shape += (sh[d+1] + misfit[d].sum(),) + logger.debug('%s[%s] :: center: %s -> %s' % (self.owner.ID, self.ID, str(self.center), str(new_center))) @@ -655,7 +688,7 @@ def reformat(self, newID=None): raise RuntimeError('Arrays larger than 50M not supported. You ' 'requested %.2fM pixels.' % megapixels) - # Apply 2d misfit + # Apply Nd misfit if self.data is not None: new_data = u.crop_pad( self.data, @@ -669,9 +702,10 @@ def reformat(self, newID=None): new_data = self.data new_shape = sh new_center = self.center - + # Deal with layermap new_layermap = sorted(layers) + if self.layermap != new_layermap: relaid_data = [] for i in new_layermap: @@ -680,7 +714,7 @@ def reformat(self, newID=None): d = new_data[self.layermap.index(i)] else: # A new layer - d = np.empty(new_shape[-2:], self.dtype) + d = np.empty(new_shape[-self.ndim:], self.dtype) d.fill(self.fill_value) relaid_data.append(d) new_data = np.array(relaid_data) @@ -688,9 +722,8 @@ def reformat(self, newID=None): self.layermap = new_layermap self.nlayers = len(new_layermap) - - # BE: set a layer index in the view the datalist access has proven to - # be too slow. + + # set layer index in the view for v in views: v.dlayer = self.layermap.index(v.layer) @@ -705,9 +738,6 @@ def reformat(self, newID=None): # A storage is "distributed" if and only if layer maps are different across nodes. self._update_distributed() - # make datalist - #self._make_datalist() - def _update_distributed(self): self.distributed = False if u.parallel.MPIenabled: @@ -742,7 +772,7 @@ def _to_phys(self, pix): @property def psize(self): """ - :returns: The pixel size. + The pixel size. """ return self._psize @@ -751,7 +781,7 @@ def psize(self, v): """ Set the pixel size, and update all the internal variables. """ - self._psize = u.expect2(v) + self._psize = u.expectN(v, self.ndim) self._origin = - self._center * self._psize self.update() @@ -767,7 +797,7 @@ def origin(self, v): """ Set the origin and update all the internal variables. """ - self._origin = u.expect2(v) + self._origin = u.expectN(v, self.ndim) self._center = - self._origin / self._psize self.update() @@ -784,7 +814,7 @@ def center(self, v): """ Set the center and update all the internal variables. """ - self._center = u.expect2(v) + self._center = u.expectN(v, self.ndim) self._origin = - self._center * self._psize self.update() @@ -826,8 +856,8 @@ def zoom_to_psize(self, new_psize, **kwargs): new_psize : scalar or array_like new pixel size """ - new_psize = u.expect2(new_psize) - sh = np.asarray(self.shape[-2:]) + new_psize = u.expectN(new_psize, self.ndim) + sh = np.asarray(self.shape[1:]) # psize is quantized new_sh = np.round(self.psize / new_psize * sh) new_psize = self.psize / new_sh * sh @@ -839,7 +869,7 @@ def zoom_to_psize(self, new_psize, **kwargs): # Zoom data buffer. # Could be that it is faster and cleaner to loop over first axis zoom = new_sh / sh - self.fill(u.zoom(self.data, [1.0, zoom[0], zoom[1]], **kwargs)) + self.fill(u.zoom(self.data, [1.0] + [z for z in zoom], **kwargs)) self._psize = new_psize self.zoom_cycle += 1 # !!! BUG: Unresolved attribute reference @@ -855,15 +885,15 @@ def grids(self): """ Returns ------- - x, y: ndarray + x, y : ndarray grids in the shape of internal buffer """ sh = self.data.shape - nm = np.indices(sh)[-2:] - flat = nm.reshape((2, self.data.size)) + nm = np.indices(sh)[1:] + flat = nm.reshape((self.ndim, self.data.size)) c = self._to_phys(flat.T).T - c = c.reshape((2,) + sh) - return c[0], c[1] + c = c.reshape((self.ndim,) + sh) + return tuple(c) def get_view_coverage(self): """ @@ -949,16 +979,18 @@ def formatted_report(self, table_format=None, offset=8, align='right', for key, column in fr.table: if str(key) == 'shape': dct[key] = tuple(self.data.shape) - info = '%d * %d * %d' % dct[key] + info = ('%d' + ' * %d' * self.ndim) % dct[key] elif str(key) == 'psize': + dec = np.floor(np.log10(self.psize).min()) dct[key] = tuple(self.psize) - info = '%.2e * %.2e' % tuple(dct[key]) - info = info.split('e', 1)[0] + info.split('e', 1)[1][3:] + info = '*'.join(['%.1f' % p for p in self.psize * 10**(-dec)]) + info += 'e%d' % dec elif str(key) == 'dimension': - dct[key] = (self.psize[0] * self.data.shape[-2], - self.psize[1] * self.data.shape[-1]) - info = '%.2e*%.2e' % tuple(dct[key]) - info = info.split('e', 1)[0] + info.split('e', 1)[1][3:] + r = self.psize * self.data.shape[1:] + dec = np.floor(np.log10(r).min()) + dct[key] = tuple(r) + info = '*'.join(['%.1f' % p for p in r * 10**(-dec)]) + info += 'e%d' % dec elif str(key) == 'memory': dct[key] = float(self.data.nbytes) / 1e6 info = '%.1f' % dct[key] @@ -1001,9 +1033,14 @@ def __getitem__(self, v): # return shift(self.data[v.slayer, v.roi[0, 0]:v.roi[1, 0], # v.roi[0, 1]:v.roi[1, 1]], v.sp) if isinstance(v, View): - return shift(self.data[ - v.dlayer, v.dlow[0]:v.dhigh[0], v.dlow[1]:v.dhigh[1]], - v.sp) + if self.ndim == 2: + return shift(self.data[ + v.dlayer, v.dlow[0]:v.dhigh[0], v.dlow[1]:v.dhigh[1]], + v.sp) + elif self.ndim == 3: + return shift(self.data[ + v.dlayer, v.dlow[0]:v.dhigh[0], v.dlow[1]:v.dhigh[1], + v.dlow[2]:v.dhigh[2]], v.sp) elif v in self.layermap: return self.data[self.layermap.index(v)] else: @@ -1032,8 +1069,30 @@ def __setitem__(self, v, newdata): # self.data[v.slayer, v.roi[0, 0]:v.roi[1, 0], # v.roi[0, 1]:v.roi[1, 1]] = shift(newdata, -v.sp) if isinstance(v, View): - self.data[v.dlayer, v.dlow[0]:v.dhigh[0], v.dlow[1]:v.dhigh[1]] = ( - shift(newdata, -v.sp)) + # there must be a nicer way to do this, numpy.take is nearly + # right, but returns copies and not views. + if self.ndim == 2: + self.data[v.dlayer, + v.dlow[0]:v.dhigh[0], + v.dlow[1]:v.dhigh[1]] = (shift(newdata, -v.sp)) + elif self.ndim == 3: + self.data[v.dlayer, + v.dlow[0]:v.dhigh[0], + v.dlow[1]:v.dhigh[1], + v.dlow[2]:v.dhigh[2]] = (shift(newdata, -v.sp)) + elif self.ndim == 4: + self.data[v.dlayer, + v.dlow[0]:v.dhigh[0], + v.dlow[1]:v.dhigh[1], + v.dlow[2]:v.dhigh[2], + v.dlow[3]:v.dhigh[3]] = (shift(newdata, -v.sp)) + elif self.ndim == 5: + self.data[v.dlayer, + v.dlow[0]:v.dhigh[0], + v.dlow[1]:v.dhigh[1], + v.dlow[2]:v.dhigh[2], + v.dlow[3]:v.dhigh[3], + v.dlow[4]:v.dhigh[4]] = (shift(newdata, -v.sp)) elif v in self.layermap: self.data[self.layermap.index(v)] = newdata else: @@ -1058,7 +1117,7 @@ def shift(v, sp): class View(Base): """ - A "window" on a Container. + A 'window' on a Container. A view stores all the slicing information to extract a 2D piece of Container. @@ -1067,16 +1126,29 @@ class View(Base): ---- The final structure of this class is yet up to debate and the constructor signature may change. Especially since - "DEFAULT_ACCESSRULE" is yet so small, its contents could be + 'DEFAULT_ACCESSRULE' is yet so small, its contents could be incorporated in the constructor call. """ + _fields = Base._fields + \ + [('active', 'b1'), + ('dlayer',' %s[%s] : shape = %s layer = %s coord = %s' % (self.owner.ID, self.storage.ID, self.ID, @@ -1207,6 +1276,35 @@ def __str__(self): else: return first + '\n ACTIVE : slice = %s' % str(self.slice) + @property + def active(self): + return self._record['active'] + + @active.setter + def active(self,v): + self._record['active'] = v + + @property + def dlayer(self): + return self._record['dlayer'] + + @dlayer.setter + def dlayer(self,v): + self._record['dlayer'] = v + + + @property + def layer(self): + return self._record['layer'] + + @layer.setter + def layer(self, v): + self._record['layer'] = v + + @property + def ndim(self): + return self._ndim + @property def slice(self): """ @@ -1219,9 +1317,10 @@ def slice(self): # else: # slayer = self.storage.layermap.index(self.layer) - return (self.dlayer, - slice(self.dlow[0], self.dhigh[0]), - slice(self.dlow[1], self.dhigh[1])) + res = (self.dlayer,) + for d in range(self.ndim): + res += (slice(self.dlow[d], self.dhigh[d]),) + return res @property def pod(self): @@ -1230,8 +1329,22 @@ def pod(self): This is a common call in the code and has therefore found its way here. May return ``None`` if there is no pod connected. """ - return self._pod() # weak reference + if isinstance(self._pod, weakref.ref): + return self._pod() # weak reference + else: + return self._pod # return self.pods.values()[0] + + @property + def pods(self): + """ + Returns all :any:`POD`s still connected to this view as a dict. + """ + if self._pods is not None: + return self._pods + else: + pod = self.pod + return {} if pod is None else {pod.ID: pod} @property def data(self): @@ -1252,7 +1365,9 @@ def shape(self): """ Two dimensional shape of View. """ - return self._arint[0] if (self._arint[0] > 0).all() else None + + sh = self._record['shape'] + return None if (sh==0).all() else sh[:self._ndim] @shape.setter def shape(self, v): @@ -1260,75 +1375,81 @@ def shape(self, v): Set two dimensional shape of View. """ if v is None: - self._arint[0] = u.expect2(0) + self._record['shape'][:] = 0 + elif np.isscalar(v): + sh = (int(v),) * self.owner.ndim + self._ndim = len(sh) + self._record['shape'][:len(sh)] = sh else: - self._arint[0] = u.expect2(v) + self._ndim = len(v) + self._record['shape'][:len(v)] = v @property def dlow(self): """ Low side of the View's data range. """ - return self._arint[1] + return self._record['dlow'][:self._ndim] @dlow.setter def dlow(self, v): """ Set low side of the View's data range. """ - self._arint[1] = v + self._record['dlow'][:self._ndim] = v @property def dhigh(self): """ High side of the View's data range. """ - return self._arint[2] + return self._record['dhigh'][:self._ndim] @dhigh.setter def dhigh(self, v): """ Set high side of the View's data range. """ - self._arint[2] = v + self._record['dhigh'][:self._ndim] = v @property def dcoord(self): """ Center coordinate (index) in data buffer. """ - return self._arint[3] + return self._record['dcoord'][:self._ndim] @dcoord.setter def dcoord(self, v): """ Set high side of the View's data range. """ - self._arint[3] = v + self._record['dcoord'][:self._ndim] = v @property def psize(self): """ Pixel size of the View. """ - return self._arfloat[0] if (self._arfloat[0] > 0.).all() else None + ps = self._record['psize'][:self._ndim] + return ps if (ps > 0.).all() else None @psize.setter def psize(self, v): """ - Set pixel size + Set pixel size. """ if v is None: - self._arfloat[0] = u.expect2(0.) + self._record['psize'][:] = 0. else: - self._arfloat[0] = u.expect2(v) + self._record['psize'][:self._ndim] = u.expectN(v, self._ndim) @property def coord(self): """ The View's physical coordinate (meters) """ - return self._arfloat[1] + return self._record['coord'][:self._ndim] @coord.setter def coord(self, v): @@ -1336,11 +1457,11 @@ def coord(self, v): Set the View's physical coordinate (meters) """ if v is None: - self._arfloat[1] = u.expect2(0.) + self._record['coord'][:] = 0. elif type(v) is not np.ndarray: - self._arfloat[1] = u.expect2(v) + self._record['coord'][:self._ndim] = u.expectN(v, self._ndim) else: - self._arfloat[1] = v + self._record['coord'][:self._ndim] = v @property def sp(self): @@ -1348,7 +1469,7 @@ def sp(self): The subpixel difference (meters) between physical coordinate and data coordinate. """ - return self._arfloat[2] + return self._record['sp'][:self._ndim] @sp.setter def sp(self, v): @@ -1357,11 +1478,18 @@ def sp(self, v): and data coordinate. """ if v is None: - self._arfloat[2] = u.expect2(0.) + self._record['sp'][:] = 0. elif type(v) is not np.ndarray: - self._arfloat[2] = u.expect2(v) + self._record['sp'][:self._ndim] = u.expectN(v, self._ndim) else: - self._arfloat[2] = v + self._record['sp'][:self._ndim] = v + + @property + def pcoord(self): + """ + The physical coordinate in pixel space + """ + return self.dcoord + self.sp class Container(Base): @@ -1402,32 +1530,31 @@ class Container(Base): """ _PREFIX = CONTAINER_PREFIX - def __init__(self, ptycho=None, ID=None, data_type='complex', **kwargs): + def __init__(self, owner=None, ID=None, data_type='complex', data_dims=2): """ Parameters ---------- ID : str or int - A unique ID, managed by the parent + A unique ID, managed by the parent. - ptycho : Ptycho - The instance of Ptycho associated with this pod. + owner : Base + A possible subclass of :any:`Base` that holds this Container. data_type : str or numpy.dtype data type - either a numpy.dtype object or 'complex' or 'real' (precision is taken from ptycho.FType or ptycho.CType) """ - # if ptycho is None: - # ptycho = ptypy.currentPtycho - - super(Container, self).__init__(ptycho, ID) - # if len(kwargs) > 0: - # self._initialize(**kwargs) - # def _initialize(self, original=None, data_type='complex'): + super(Container, self).__init__(owner, ID) self.data_type = data_type + #: Default data dimensionality for Views and Storages. + self.ndim = data_dims + if self.ndim not in (2, 3): + logger.warn('Container untested for `data_dim` other than 2 or 3') + # Prepare for copy # self.original = original if original is not None else self self.original = self @@ -1543,7 +1670,7 @@ def nbytes(self): sz += s.data.nbytes return sz - def views_in_storage(self, s, active=True): + def views_in_storage(self, s, active_only=True): """ Return a list of views on :any:`Storage` `s`. @@ -1551,10 +1678,10 @@ def views_in_storage(self, s, active=True): ---------- s : Storage The storage to look for. - active : True or False + active_only : True or False If True (default), return only active views. """ - if active: + if active_only: return [v for v in self.original.V.values() if v.active and (v.storageID == s.ID)] else: @@ -1582,7 +1709,7 @@ def copy(self, ID=None, fill=None, dtype=None): # Create new container data_type = self.data_type if dtype is None else dtype - new_cont = type(self)(ptycho=self.owner, + new_cont = type(self)(self.owner, ID=ID, data_type=data_type) new_cont.original = self @@ -1607,7 +1734,6 @@ def fill(self, fill=0.0): self += fill for s in self.storages.itervalues(): s.fill(fill) - s._make_datalist() def allreduce(self, op=None): """ @@ -1907,14 +2033,18 @@ class POD(Base): _PREFIX = POD_PREFIX - def __init__(self, ptycho=None, ID=None, views=None, geometry=None, - **kwargs): + def __init__(self, ptycho=None, model=None, ID=None, views=None, + geometry=None, **kwargs): """ 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. @@ -1926,10 +2056,7 @@ def __init__(self, ptycho=None, ID=None, views=None, geometry=None, """ super(POD, self).__init__(ptycho, ID, False) - # if len(kwargs) > 0: - # self._initialize(**kwargs) - - # def _initialize(self, views=None, geometry=None): #,meta=None): + self.model = model # other defaults: self.is_empty = False @@ -1943,8 +2070,16 @@ def __init__(self, ptycho=None, ID=None, views=None, geometry=None, for v in self.V.values(): if v is None: continue - v.pods[self.ID] = self - v._pod = weakref.ref(self) + if v._pod is None: + # you are first + v._pod = weakref.ref(self) + else: + # View has at least one POD connected + if v._pods is None: + v._pods = weakref.WeakValueDictionary() + # register the older view + v._pods[v.pod.ID] = v.pod + v._pods[self.ID] = self #: :any:`Geo` instance with propagators self.geometry = geometry @@ -2008,11 +2143,13 @@ def object(self): if not self.is_empty: return self.ob_view.data else: + # Empty probe means no object (perfect transmission) return np.ones(self.geometry.shape, dtype=self.owner.CType) @object.setter def object(self, v): - self.ob_view.data = v + if not self.is_empty: + self.ob_view.data = v @property def probe(self): diff --git a/ptypy/core/data.py b/ptypy/core/data.py index c593d6e12..ecaf476d5 100644 --- a/ptypy/core/data.py +++ b/ptypy/core/data.py @@ -17,23 +17,17 @@ :copyright: Copyright 2014 by the PTYPY team, see AUTHORS. :license: GPLv2, see LICENSE for details. """ -from ..utils import parallel -from ptypy import resources -from ptypy.core import xy import numpy as np import os import h5py -if __name__ == "__main__": - from ptypy import utils as u - from ptypy import io - from ptypy.core import geometry - from ptypy.utils.verbose import logger, log, headerline -else: - from .. import utils as u - from .. import io - from .. import resources - from ..utils.verbose import logger, log, headerline - import geometry +from . import geometry +from . import xy +from .. import utils as u +from .. import io +from .. import resources +from ..utils import parallel +from ..utils.verbose import logger, log, headerline +from .. import defaults_tree PTYD = dict( # frames, positions @@ -45,60 +39,17 @@ ) """ Basic Structure of a .ptyd datafile """ -META = dict( - # Label will be set internally - label=None, - # A unique label of user choice - experimentID=None, - version='0.1', - shape=None, - psize=None, - # lam=None, - energy=None, - center=None, - distance=None, -) - -GENERIC = u.Param( - # Filename (e.g. 'foo.ptyd') - dfile=None, - # Format for chunk file appendix. - chunk_format='.chunk%02d', - # 2-tuple or int for the desired fina frame size - # roi=None, - # Saving option: None, 'merge', 'append', 'extlink' - save=None, - # Auto center: if False, no automatic center, None only - # if center is None, True it will be enforced - auto_center=None, - # Parallel loading: None, 'data', 'common', 'all' - load_parallel='data', - # Rebin diffraction data - rebin=None, - # Switching orientation : None, int or 3-tuple switch - # Actions are (transpose, invert rows, invert cols) - orientation=None, - # Minimum number of frames of one chunk if not at end of scan - min_frames=1, - # Theoretical position list (This input parameter may get deprecated) - positions_theory=None, - # Total number of frames to be prepared - num_frames=None, - recipe={}, -) -""" Default data parameters. See :py:data:`.scan.data` - and a short listing below """ - -GENERIC.update(META) WAIT = 'msg1' EOS = 'msgEOS' CODES = {WAIT: 'Scan unfinished. More frames available after a pause', EOS: 'End of scan reached'} -__all__ = ['GENERIC', 'PtyScan', 'PTYD', 'PtydScan', 'MoonFlowerScan'] +__all__ = ['PtyScan', 'PTYD', 'PtydScan', + 'MoonFlowerScan'] +@defaults_tree.parse_doc('scandata.PtyScan') class PtyScan(object): """ PtyScan: A single ptychography scan, created on the fly or read from file. @@ -112,44 +63,208 @@ class PtyScan(object): - On-the-fly support in form of chunked data. - mpi capable, child classes should not worry about mpi + Default data parameters. See :py:data:`.scan.data` + + Defaults: + + [name] + default = PtyScan + type = str + help = + + [dfile] + type = file + default = None + help = Prepared data file path + doc = If source was ``None`` or ``'file'``, data will be loaded from this file and processing as + well as saving is deactivated. If source is the name of an experiment recipe or path to a + file, data will be saved to this file + userlevel = 0 + + [chunk_format] + type = str + default = .chunk%02d + help = Appendix to saved files if save == 'link' + doc = + userlevel = 2 + + [save] + type = str + default = None + help = Saving mode + doc = Mode to use to save data to file. + + - ``None``: No saving + - ``'merge'``: attemts to merge data in single chunk **[not implemented]** + - ``'append'``: appends each chunk in master \*.ptyd file + - ``'link'``: appends external links in master \*.ptyd file and stores chunks separately + + in the path given by the link. Links file paths are relative to master file. + userlevel = 1 + + [auto_center] + type = bool + default = None + help = Determine if center in data is calculated automatically + doc = + - ``False``, no automatic centering + - ``None``, only if :py:data:`center` is ``None`` + - ``True``, it will be enforced + userlevel = 0 + + [load_parallel] + type = str + default = data + help = Determines what will be loaded in parallel + doc = Choose from ``None``, ``'data'``, ``'common'``, ``'all'`` + + [rebin] + type = int + default = None + help = Rebinning factor + doc = Rebinning factor for the raw data frames. ``'None'`` or ``1`` both mean *no binning* + userlevel = 1 + lowlim = 1 + uplim = 8 + + [orientation] + type = int, tuple, list + default = None + help = Data frame orientation + doc = Choose + + - ``None`` or ``0``: correct orientation + - ``1``: invert columns (numpy.flip_lr) + - ``2``: invert rows (numpy.flip_ud) + - ``3``: invert columns, invert rows + - ``4``: transpose (numpy.transpose) + - ``4+i``: tranpose + other operations from above + + Alternatively, a 3-tuple of booleans may be provided ``(do_transpose, + do_flipud, do_fliplr)`` + userlevel = 1 + + [min_frames] + type = int + default = 1 + help = Minimum number of frames loaded by each node + doc = + userlevel = 2 + + [positions_theory] + type = ndarray + default = None + help = Theoretical positions for this scan + doc = If provided, experimental positions from :py:class:`PtyScan` subclass will be ignored. If data + preparation is called from Ptycho instance, the calculated positions from the + :py:func:`ptypy.core.xy.from_pars` dict will be inserted here + userlevel = 2 + + [num_frames] + type = int + default = None + help = Maximum number of frames to be prepared + doc = If `positions_theory` are provided, num_frames will be ovverriden with the number of + positions available + userlevel = 1 + + [label] + type = str + default = None + help = The scan label + doc = Unique string identifying the scan + userlevel = 1 + + [experimentID] + type = str + default = None + help = Name of the experiment + doc = If None, a default value will be provided by the recipe. **unused** + userlevel = 2 + + [version] + type = float + default = 0.1 + help = TODO: Explain this and decide if it is a user parameter. + doc = + userlevel = 2 + + [shape] + type = int, tuple + default = 256 + help = Shape of the region of interest cropped from the raw data. + doc = Cropping dimension of the diffraction frame + Can be None, (dimx, dimy), or dim. In the latter case shape will be (dim, dim). + userlevel = 1 + + [center] + type = tuple, str + default = 'fftshift' + help = Center (pixel) of the optical axes in raw data + doc = If ``None``, this parameter will be set by :py:data:`~.scan.data.auto_center` or elsewhere + userlevel = 1 + + [psize] + type = float, tuple + default = 0.000172 + help = Detector pixel size + doc = Dimensions of the detector pixels (in meters) + userlevel = 0 + lowlim = 0 + + [distance] + type = float + default = 7.19 + help = Sample to detector distance + doc = In meters. + userlevel = 0 + lowlim = 0 + + [energy] + type = float + default = 7.2 + help = Photon energy of the incident radiation in keV + doc = + userlevel = 0 + lowlim = 0 + + [add_poisson_noise] + default = True + type = bool + help = Decides whether the scan should have poisson noise or not """ - DEFAULT = GENERIC.copy() WAIT = WAIT EOS = EOS CODES = CODES + METAKEYS = ['version', 'num_frames', 'label', 'shape', 'psize', 'energy', 'center', 'distance'] + """ Keys to store in meta param """ + def __init__(self, pars=None, **kwargs): # filename='./foo.ptyd', shape=None, save=True): """ - Class creation with minimum set of parameters, see :py:data:`GENERIC` + Class creation with minimum set of parameters, see :py:data:`PtyScan.DEFAULT` Please note that class creation is not meant to load data. Call :py:data:`initialize` to begin loading and data file creation. """ # Load default parameter structure - info = u.Param(self.DEFAULT.copy()) - - # FIXME this overwrites the child's recipe defaults - info.update(pars, in_place_depth=1) - info.update(kwargs) - - # validate(pars, '.scan.preparation') - - # Prepare meta data - self.meta = u.Param(META.copy()) + p = self.DEFAULT.copy(99) + p.update(pars) + p.update(kwargs) # Attempt to get number of frames. - self.num_frames = info.num_frames + self.num_frames = p.num_frames """ Total number of frames to prepare / load. Set by :py:data:`~.scan.data.num_frames` """ - self.min_frames = info.min_frames * parallel.size + self.min_frames = p.min_frames * parallel.size """ Minimum number of frames to prepare / load with call of :py:meth:`auto` """ - if info.positions_theory is not None: - num = len(info.positions_theory) + if p.positions_theory is not None: + num = len(p.positions_theory) logger.info('Theoretical positions are available. ' 'There will be %d frames.' % num) logger.info( @@ -158,33 +273,17 @@ def __init__(self, pars=None, **kwargs): 'Former input value of frame number `num_frames` %s is ' 'overridden to %d.' % (str(self.num_frames), num)) self.num_frames = num - """ - # check if we got information on geometry from ptycho - if info.geometry is not None: - for k, v in info.geometry.items(): - # FIXME: This is a bit ugly - - # some parameters are added to info without documentation. - info[k] = v if info.get(k) is None else None - # FIXME: This should probably be done more transparently: - # it is not clear for the user that info.roi - # has precedence over geometry.N - if info.roi is None: - info.roi = u.expect2(info.geometry.N) - """ + # None for rebin should be allowed, as in "don't rebin". - if info.rebin is None: - info.rebin = 1 + if p.rebin is None: + p.rebin = 1 - self.info = info - """:any:`Param` container that stores all input parameters.""" + self.info = p + """:py:class:`Param` container that stores all input parameters.""" # Print a report log(4, 'Ptypy Scan instance got the following parameters:') - log(4, u.verbose.report(info)) - - # Dump all input parameters as class attributes. - # FIXME: This duplication of parameters can lead to much confusion... - # self.__dict__.update(info) + log(4, u.verbose.report(p)) # Check MPI settings lp = str(self.info.load_parallel) @@ -205,12 +304,9 @@ def __init__(self, pars=None, **kwargs): self.dfile = None self.save = self.info.save - # Copy all values for meta - for k in self.meta.keys(): - self.meta[k] = self.info[k] - # self.center = None # Center will be set later - # self.roi = self.info.roi #None # ROI will be set later - # self.shape = None + # Construct meta + self.meta = u.Param({k: self.info[k] for k in self.METAKEYS}) + self.orientation = self.info.orientation self.rebin = self.info.rebin @@ -218,9 +314,6 @@ def __init__(self, pars=None, **kwargs): self._flags = np.array([0, 0, 0], dtype=int) self.is_initialized = False - # post init method call - self.post_init() - def initialize(self): """ Begins the Data preparation and intended as the first method @@ -234,7 +327,8 @@ def initialize(self): * Sets :py:attr:`num_frames` if needed * Calls :py:meth:`post_initialize` """ - logger.info(headerline('Enter PtyScan.initialize()', 'l')) + logger.info(headerline('Enter %s.initialize()' + % self.__class__.__name__, 'l')) # Prepare writing to file if self.info.save is not None: @@ -282,6 +376,9 @@ def initialize(self): if self.has_weight2d: logger.info('shape = '.rjust(29) + str(self.weight2d.shape)) + # FIXME: Saving weight to info. This is not ideal, info is optional + self.info.weight2d = self.has_weight2d + logger.info('All experimental positions : ' + str(self.has_positions)) if self.has_positions: logger.info('shape = '.rjust(29) + str(positions.shape)) @@ -479,13 +576,6 @@ def post_initialize(self): """ pass - def post_init(self): - """ - Placeholder. Called at the end of construction by all - processes. - """ - pass - def _mpi_check(self, chunksize, start=None): """ Executes the check() function on master node and communicates @@ -567,240 +657,228 @@ def get_data_chunk(self, chunksize, start=None): if msg in [EOS, WAIT]: logger.info(CODES[msg]) return msg - else: - start, step = msg - # Get scan point index distribution - indices = self._mpi_indices(start, step) + start, step = msg - # What did we get? - data, positions, weights = self._mpi_pipeline_with_dictionaries( - indices) - # All these dictionaries could be empty - # Fill weights dictionary with references to the weights in common + # Get scan point index distribution + indices = self._mpi_indices(start, step) - has_data = (len(data) > 0) - has_weights = (len(weights) > 0) and len(weights.values()[0]) > 0 + # What did we get? + data, positions, weights = self._mpi_pipeline_with_dictionaries( + indices) + # All these dictionaries could be empty + # Fill weights dictionary with references to the weights in common - if has_data: - dsh = np.array(data.values()[0].shape[-2:]) - else: - dsh = np.array([0, 0]) + has_data = (len(data) > 0) + has_weights = (len(weights) > 0) and len(weights.values()[0]) > 0 - # Communicate result - dsh[0] = parallel.MPImax([dsh[0]]) - dsh[1] = parallel.MPImax([dsh[1]]) + if has_data: + dsh = np.array(data.values()[0].shape[-2:]) + else: + dsh = np.array([0, 0]) - if not has_weights: - # Peak at first item - if self.has_weight2d: - altweight = self.weight2d - else: - try: - altweight = self.meta.weight2d - except: - altweight = np.ones(dsh) - weights = dict.fromkeys(data.keys(), altweight) - - assert len(weights) == len(data), ( - 'Data and Weight frames unbalanced %d vs %d' - % (len(data), len(weights))) - - sh = self.info.shape - # Adapt roi if not set - if sh is None: - logger.info('ROI not set. Using full frame shape of (%d, %d).' - % tuple(dsh)) - sh = dsh - else: - sh = u.expect2(sh) - - # Only allow square slices in data - if sh[0] != sh[1]: - roi = u.expect2(sh.min()) - logger.warning('Asymmetric data ROI not allowed. Setting ROI ' - 'from (%d, %d) to (%d, %d).' - % (sh[0], sh[1], roi[0], roi[1])) - sh = roi - - self.info.shape = sh - - cen = self.info.center - if str(cen) == cen: - cen = geometry.translate_to_pix(sh, cen) - - auto = self.info.auto_center - # Get center in diffraction image - if auto is None or auto is True: - auto_cen = self._mpi_autocenter(data, weights) - else: - auto_cen = None + # Communicate result + dsh[0] = parallel.MPImax([dsh[0]]) + dsh[1] = parallel.MPImax([dsh[1]]) - if cen is None and auto_cen is not None: - logger.info('Setting center for ROI from %s to %s.' - % (str(cen), str(auto_cen))) - cen = auto_cen - elif cen is None and auto is False: - cen = dsh // 2 + if not has_weights: + # Peak at first item + if self.has_weight2d: + altweight = self.weight2d else: - # Center is number or tuple - cen = u.expect2(cen[-2:]) - if auto_cen is not None: - logger.info('ROI center is %s, automatic guess is %s.' - % (str(cen), str(auto_cen))) - - # It is important to set center again in order to NOT have - # different centers for each chunk, the downside is that the center - # may be off if only a few diffraction patterns are used for the - # analysis. In such case it is beneficial to set the center in the - # parameters - self.info.center = cen - - # Make sure center is in the image frame - assert (cen > 0).all() and (dsh - cen > 0).all(), ( - 'Optical axes (center = (%.1f, %.1f) outside diffraction image ' - 'frame (%d, %d).' % tuple(cen) + tuple(dsh)) - - # Determine if the arrays require further processing - do_flip = (self.orientation is not None - and np.array(self.orientation).any()) - do_crop = (np.abs(sh - dsh) > 0.5).any() - do_rebin = self.rebin is not None and (self.rebin != 1) - - if do_flip or do_crop or do_rebin: - logger.info( - 'Enter preprocessing ' - '(crop/pad %s, rebin %s, flip/rotate %s) ... \n' - % (str(do_crop), str(do_rebin), str(do_flip))) - - # We proceed with numpy arrays.That is probably now more memory - # intensive but shorter in writing - if has_data: - d = np.array([data[ind] for ind in indices.node]) - w = np.array([weights[ind] for ind in indices.node]) - else: - d = np.ones((1,) + tuple(dsh)) - w = np.ones((1,) + tuple(dsh)) - - # Crop data - d, tmp = u.crop_pad_symmetric_2d(d, sh, cen) - - # Check if provided mask has the same shape as data, if not, - # use the mask's center for cropping the mask. The latter is - # also needed if no mask is provided, as weights is then - # created using the requested cropping shape and thus might - # have a different center than the raw data. - # NOTE: Maybe distinguish between no mask provided and mask - # with wrong size in warning - if (dsh == np.array(w[0].shape)).all(): - w, cen = u.crop_pad_symmetric_2d(w, sh, cen) - else: - logger.warning('Mask does not have the same shape as data. ' - 'Will use mask center for cropping mask.') - cen = np.array(w[0].shape) // 2 - w, cen = u.crop_pad_symmetric_2d(w, sh, cen) - - # Flip, rotate etc. - d, tmp = u.switch_orientation(d, self.orientation, cen) - w, cen = u.switch_orientation(w, self.orientation, cen) - - # Rebin, check if rebinning is neither to strong nor impossible - rebin = self.rebin - if rebin <= 1: - pass - elif (rebin in range(2, 6) - and (((sh / float(rebin)) % 1) == 0.0).all()): - mask = w > 0 - d = u.rebin_2d(d, rebin) - w = u.rebin_2d(w, rebin) - mask = u.rebin_2d(mask, rebin) - # We keep only the pixels that do not include a masked pixel - # w[mask < mask.max()] = 0 - # TODO: apply this operation when weights actually are weights - w = (mask == mask.max()) - else: - raise RuntimeError( - 'Binning (%d) is to large or incompatible with array ' - 'shape (%s).' % (rebin, str(tuple(sh)))) + try: + altweight = self.info.weight2d + except: + altweight = np.ones(dsh) + weights = dict.fromkeys(data.keys(), altweight) + + assert len(weights) == len(data), ( + 'Data and Weight frames unbalanced %d vs %d' + % (len(data), len(weights))) + + sh = self.info.shape + # Adapt roi if not set + if sh is None: + logger.info('ROI not set. Using full frame shape of (%d, %d).' + % tuple(dsh)) + sh = dsh + else: + sh = u.expect2(sh) + + # Only allow square slices in data + if sh[0] != sh[1]: + roi = u.expect2(sh.min()) + logger.warning('Asymmetric data ROI not allowed. Setting ROI ' + 'from (%d, %d) to (%d, %d).' + % (sh[0], sh[1], roi[0], roi[1])) + sh = roi + + self.info.shape = sh + + cen = self.info.center + if str(cen) == cen: + cen = geometry.translate_to_pix(dsh, cen) + + auto = self.info.auto_center + # Get center in diffraction image + if auto is None or auto is True: + auto_cen = self._mpi_autocenter(data, weights) + else: + auto_cen = None + + if cen is None and auto_cen is not None: + logger.info('Setting center for ROI from %s to %s.' + % (str(cen), str(auto_cen))) + cen = auto_cen + elif cen is None and auto is False: + cen = dsh // 2 + else: + # Center is number or tuple + cen = u.expect2(cen[-2:]) + if auto_cen is not None: + logger.info('ROI center is %s, automatic guess is %s.' + % (str(cen), str(auto_cen))) - # restore contiguity of the cropped/padded/rotated/flipped array - d = np.ascontiguousarray(d) + # It is important to set center again in order to NOT have + # different centers for each chunk, the downside is that the center + # may be off if only a few diffraction patterns are used for the + # analysis. In such case it is beneficial to set the center in the + # parameters + self.info.center = cen + + # Make sure center is in the image frame + assert (cen > 0).all() and (dsh - cen > 0).all(), ( + 'Optical axes (center = (%.1f, %.1f) outside diffraction image ' + 'frame (%d, %d).' % (tuple(cen) + tuple(dsh))) + + # Determine if the arrays require further processing + do_flip = (self.orientation is not None + and np.array(self.orientation).any()) + do_crop = (np.abs(sh - dsh) > 0.5).any() + do_rebin = self.rebin is not None and (self.rebin != 1) + + if do_flip or do_crop or do_rebin: + logger.info( + 'Enter preprocessing ' + '(crop/pad %s, rebin %s, flip/rotate %s) ... \n' + % (str(do_crop), str(do_rebin), str(do_flip))) - if has_data: - # Translate back to dictionaries - data = dict(zip(indices.node, d)) - weights = dict(zip(indices.node, w)) + # We proceed with numpy arrays.That is probably now more memory + # intensive but shorter in writing + if has_data: + d = np.array([data[ind] for ind in indices.node]) + w = np.array([weights[ind] for ind in indices.node]) + else: + d = np.ones((1,) + tuple(dsh)) + w = np.ones((1,) + tuple(dsh)) + + # Crop data + d, tmp = u.crop_pad_symmetric_2d(d, sh, cen) + + # Check if provided mask has the same shape as data, if not, + # use the mask's center for cropping the mask. The latter is + # also needed if no mask is provided, as weights is then + # created using the requested cropping shape and thus might + # have a different center than the raw data. + # NOTE: Maybe distinguish between no mask provided and mask + # with wrong size in warning + + if (dsh == np.array(w[0].shape)).all(): + w, cen = u.crop_pad_symmetric_2d(w, sh, cen) + else: + logger.warning('Mask does not have the same shape as data. ' + 'Will use mask center for cropping mask.') + cen = np.array(w[0].shape) // 2 + w, cen = u.crop_pad_symmetric_2d(w, sh, cen) + + # Flip, rotate etc. + d, tmp = u.switch_orientation(d, self.orientation, cen) + w, cen = u.switch_orientation(w, self.orientation, cen) + + # Rebin, check if rebinning is neither to strong nor impossible + rebin = self.rebin + if rebin <= 1: + pass + elif (rebin in range(2, 6) + and (((sh / float(rebin)) % 1) == 0.0).all()): + mask = w > 0 + d = u.rebin_2d(d, rebin) + w = u.rebin_2d(w, rebin) + mask = u.rebin_2d(mask, rebin) + # We keep only the pixels that do not include a masked pixel + # w[mask < mask.max()] = 0 + # TODO: apply this operation when weights actually are weights + w = (mask == mask.max()) + else: + raise RuntimeError( + 'Binning (%d) is to large or incompatible with array ' + 'shape (%s).' % (rebin, str(tuple(sh)))) - # Adapt meta info - self.meta.center = cen / float(self.rebin) - self.meta.shape = u.expect2(sh) / self.rebin + # restore contiguity of the cropped/padded/rotated/flipped array + d = np.ascontiguousarray(d) + w = np.ascontiguousarray(w) - if self.info.psize is not None: - self.meta.psize = u.expect2(self.info.psize) * self.rebin - else: - self.meta.psize = None - - # Prepare chunk of data - chunk = u.Param() - chunk.indices = indices.chunk - chunk.indices_node = indices.node - chunk.num = self.chunknum - chunk.data = data - - # If there are weights we add them to chunk, - # otherwise we push it into meta - if has_weights: - chunk.weights = weights - elif has_data: - chunk.weights = {} - self.meta.weight2d = weights.values()[0] - - # Slice positions from common if they are empty too - if positions is None or len(positions) == 0: - pt = self.info.positions_theory - if pt is not None: - chunk.positions = pt[indices.chunk] - else: - try: - chunk.positions = ( - self.info.positions_scan[indices.chunk]) - except: - logger.info('Unable to slice position information from ' - 'experimental or theoretical resource.') - chunk.positions = [None] * len(indices.chunk) + if has_data: + # Translate back to dictionaries + data = dict(zip(indices.node, d)) + weights = dict(zip(indices.node, w)) + + # Adapt geometric info + self.meta.center = cen / float(self.rebin) + self.meta.shape = u.expect2(sh) / self.rebin + + if self.info.psize is not None: + self.meta.psize = u.expect2(self.info.psize) * self.rebin + + # Prepare chunk of data + chunk = u.Param() + chunk.indices = indices.chunk + chunk.indices_node = indices.node + chunk.num = self.chunknum + chunk.data = data + + # If there are weights we add them to chunk, + # otherwise we push it into meta + if has_weights: + chunk.weights = weights + elif has_data: + chunk.weights = {} + self.weight2d = weights.values()[0] + + # Slice positions from common if they are empty too + if positions is None or len(positions) == 0: + pt = self.info.positions_theory + if pt is not None: + chunk.positions = pt[indices.chunk] else: - # A dict : sort positions to indices.chunk - # This may fail if there are less positions than scan points - # (unlikely) - chunk.positions = np.asarray( - [positions[k] for k in indices.chunk]) - # Positions complete - - # With first chunk we update meta - if self.chunknum < 1: - """ - for k, v in self.meta.items(): - # FIXME: I would like to avoid this "blind copying" - # BE: This is not a blind copy as only keys - # in META above are used - if v is None: - self.meta[k] = self.__dict__.get(k, self.info.get(k)) - else: - self.meta[k] = v - self.meta['center'] = cen - """ - - if self.info.save is not None and parallel.master: - io.h5append(self.dfile, meta=dict(self.meta)) + try: + chunk.positions = ( + self.info.positions_scan[indices.chunk]) + except: + logger.info('Unable to slice position information from ' + 'experimental or theoretical resource.') + chunk.positions = [None] * len(indices.chunk) + else: + # A dict : sort positions to indices.chunk + # This may fail if there are less positions than scan points + # (unlikely) + chunk.positions = np.asarray( + [positions[k] for k in indices.chunk]) + # Positions complete + + # With first chunk we update info + if self.chunknum < 1: + if self.info.save is not None and parallel.master: + io.h5append(self.dfile, meta=dict(self.meta)) - parallel.barrier() + parallel.barrier() - self.chunk = chunk - self.chunknum += 1 + self.chunk = chunk + self.chunknum += 1 - return chunk + return chunk - def auto(self, frames, chunk_form='dp'): + def auto(self, frames): """ Repeated calls to this function will process the data. @@ -809,16 +887,13 @@ def auto(self, frames, chunk_form='dp'): frames : int Number of frames to process. - chunk_form : str - Currently only type data package 'dp' implemented - Returns ------- variable one of the following - - None, if scan's end is not reached, + - WAIT, if scan's end is not reached, but no data could be prepared yet - - False, if scan's end is reached + - EOS, if scan's end is reached - a data package otherwise """ # attempt to get data: @@ -832,7 +907,7 @@ def auto(self, frames, chunk_form='dp'): # del self.chunk return msg else: - out = self.return_chunk_as(msg, chunk_form) + out = self._make_data_package(msg) # save chunk if self.info.save is not None: self._mpi_save_chunk(self.info.save, msg) @@ -840,15 +915,10 @@ def auto(self, frames, chunk_form='dp'): del self.chunk return out - def return_chunk_as(self, chunk, kind='dp'): + def _make_data_package(self, chunk): """ - Returns the loaded data chunk `chunk` in the format `kind`. - - For now only kind=='dp' (data package) is valid. + Returns the loaded data chunk `chunk` as a data package. """ - # This is a bit ugly now - if kind != 'dp': - raise RuntimeError('Unknown kind of chunck format: %s' % str(kind)) # The "common" part out = {'common': self.meta} @@ -867,9 +937,11 @@ def return_chunk_as(self, chunk, kind='dp'): # First look in chunk for a weight to this index, then # look for a 2d-weight in meta, then arbitrarily set # weight to ones. - w = chunk.weights.get( - index, self.meta.get('weight2d', - np.ones_like(frame['data']))) + try: + fallback = self.weight2d + except AttributeError: + fallback = np.ones_like(frame['data']) + w = chunk.weights.get(index, fallback) frame['mask'] = (w > 0) iterables.append(frame) @@ -926,7 +998,7 @@ def check(self, frames=None, start=None): for preparation and should determine if data acquisition for this scan is finished. Its main purpose is to allow for a data acquisition scheme, where the number of frames is not known - when :any:`PtyScan` is constructed, i.e. a data stream or an + when :py:class:`PtyScan` is constructed, i.e. a data stream or an on-the-fly reconstructions. Note @@ -1004,7 +1076,7 @@ def load(self, indices): Note ---- This is the *most* important method to change when subclassing - :any:`PtyScan`. Most often it suffices to override the constructor + :py:class:`PtyScan`. Most often it suffices to override the constructor and this method to create a subclass suited for a specific experiment. """ @@ -1050,17 +1122,14 @@ def _mpi_autocenter(self, data, weights): cen[k] = u.mass_center(d * (weights[k] > 0)) # For some nodes, cen may still be empty. - # Therefore we use gather_dict to be save + # Therefore we use gather_dict to be safe cen = parallel.gather_dict(cen) parallel.barrier() # Now master possesses all calculated centers if parallel.master: cen = np.array(cen.values()).mean(0) - else: - cen = np.array([0., 0.]) - - parallel.allreduce(cen) + cen = parallel.bcast(cen) return cen @@ -1145,27 +1214,36 @@ def _mpi_save_chunk(self, kind='link', chunk=None): parallel.barrier() +@defaults_tree.parse_doc('scandata.PtydScan') class PtydScan(PtyScan): """ PtyScan provided by native "ptyd" file format. + + Defaults: + + [name] + default = PtydScan + type = str + help = + doc = + + [source] + default = 'file' + type = str, None + help = Alternate source file path if data is meant to be reprocessed. + doc = `None` for input shall be deprecated in future + """ - DEFAULT = GENERIC.copy() - def __init__(self, pars=None, source=None, **kwargs): + def __init__(self, pars=None, **kwargs): """ PtyScan provided by native "ptyd" file format. - - :param source: Explicit source file. If not None or 'file', - the data may get processed depending on user input - - :param pars: Input like PtyScan """ # Create parameter set - p = u.Param(self.DEFAULT.copy()) - - # Copy the label - # if pars is not None: - # p.label = pars.get('label') + p = self.DEFAULT.copy(99) + p.update(pars) + p.update(kwargs) + source = p.source if source is None or str(source) == 'file': # This is the case of absolutely no additional work @@ -1183,7 +1261,7 @@ def __init__(self, pars=None, source=None, **kwargs): dfile = pars['dfile'] # Check for conflict - if str(u.unique_path(source)) == str(u.unique_path(dfile)): + if dfile and (str(u.unique_path(source)) == str(u.unique_path(dfile))): logger.info('Source and Sink files are the same.') dfile = os.path.splitext(dfile) dfile = dfile[0] + '_n.' + dfile[1] @@ -1203,18 +1281,20 @@ def __init__(self, pars=None, source=None, **kwargs): # At least ONE chunk must exist to ensure everything works with h5py.File(source, 'r') as f: check = f.get('chunks/0') + f.close() # Get number of frames supposedly in the file - # FIXME: try/except clause only for backward compatibilty + # FIXME: try/except clause only for backward compatibilty # for .ptyd files created priot to commit 2e626ff #try: # source_frames = f.get('info/num_frames_actual')[...].item() #except TypeError: # source_frames = len(f.get('info/positions_scan')[...]) - f.close() + #f.close() if check is None: raise IOError('Ptyd source %s contains no data. Load aborted' % source) + """ if source_frames is None: logger.warning('Ptyd source is not aware of the total' @@ -1224,22 +1304,32 @@ def __init__(self, pars=None, source=None, **kwargs): # Get meta information meta = u.Param(io.h5read(self.source, 'meta')['meta']) + if meta.get('num_frames') is None: + logger.warning('Ptyd source is not aware of the total' + 'number of diffraction frames expected') + if len(meta) == 0: logger.warning('There should be meta information in ' '%s. Something is odd here.' % source) # Update given parameters when they are None if not manipulate: - super(PtydScan, self).__init__(meta, **kwargs) + p.update(meta) else: - # Overwrite only those set to None + # Replace only None entries in p + # FIXME: + # BE: This was the former right way when the defaults + # were mostly None, now this no longer applies, unless + # defaults are overwritten to None. I guess t would be + # canonical now to overwrite the defaults in the + # docstring. But since reprocessing is rare for k, v in meta.items(): - if p.get(k) is None: # should be replace by 'unset' + if p.get(k) is None: p[k] = v - # Initialize parent class and fill self - super(PtydScan, self).__init__(p, **kwargs) - """ + super(PtydScan, self).__init__(p) + + """ if source_frames is not None: if self.num_frames is None: self.num_frames = source_frames @@ -1350,29 +1440,65 @@ def load(self, indices): return (out.get(key, {}) for key in ['data', 'positions', 'weights']) - +@defaults_tree.parse_doc('scandata.MoonFlowerScan') class MoonFlowerScan(PtyScan): """ Test PtyScan class producing a romantic ptychographic data set of a moon illuminating flowers. - """ - DEFAULT = GENERIC.copy() - DEFAULT.update(geometry.DEFAULT.copy()) - RECIPE = u.Param( - # Position distance in fraction of illumination frame - density=0.2, - photons=1e8, - psf=0. - ) + Override parent class default: + + Defaults: + + [name] + default = MoonFlowerScan + type = str + help = + doc = + + [num_frames] + default = 100 + type = int + help = Number of frames to simulate + doc = + + [shape] + type = int, tuple + default = 128 + help = Shape of the region of interest cropped from the raw data. + doc = Cropping dimension of the diffraction frame + Can be None, (dimx, dimy), or dim. In the latter case shape will be (dim, dim). + userlevel = 1 + + [density] + default = 0.2 + type = float + help = Position distance in fraction of illumination frame + + [model] + default = 'round' + type = str + help = The scan pattern + + [photons] + default = 1e8 + type = float + help = Total number of photons for Poisson noise + + [psf] + default = 0. + type = float + help = Point spread function of the detector + + """ def __init__(self, pars=None, **kwargs): """ Parent pars are for the """ - p = geometry.DEFAULT.copy() - if pars is not None: - p.update(pars) + + p = self.DEFAULT.copy(depth=99) + p.update(pars) # Initialize parent class super(MoonFlowerScan, self).__init__(p, **kwargs) @@ -1380,18 +1506,24 @@ def __init__(self, pars=None, **kwargs): # Derive geometry from input geo = geometry.Geo(pars=self.meta) - # Recipe specific things - r = self.RECIPE.copy() - r.update(self.info.recipe) - # Derive scan pattern - pos = u.Param() - pos.spacing = geo.resolution * geo.shape * r.density - pos.steps = np.int(np.round(np.sqrt(self.num_frames))) + 1 - pos.extent = pos.steps * pos.spacing - pos.model = 'round' - pos.count = self.num_frames - self.pos = xy.from_pars(pos) + if p.model is 'raster': + pos = u.Param() + pos.spacing = geo.resolution * geo.shape * p.density + pos.steps = np.int(np.round(np.sqrt(self.num_frames))) + 1 + pos.extent = pos.steps * pos.spacing + pos.model = p.model + self.num_frames = pos.steps**2 + pos.count = self.num_frames + self.pos = xy.raster_scan(pos.spacing[0], pos.spacing[1], pos.steps, pos.steps) + else: + pos = u.Param() + pos.spacing = geo.resolution * geo.shape * p.density + pos.steps = np.int(np.round(np.sqrt(self.num_frames))) + 1 + pos.extent = pos.steps * pos.spacing + pos.model = p.model + pos.count = self.num_frames + self.pos = xy.from_pars(pos) # Calculate pixel positions pixel = self.pos / geo.resolution @@ -1399,14 +1531,25 @@ def __init__(self, pars=None, **kwargs): self.pixel = np.round(pixel).astype(int) + 10 frame = self.pixel.max(0) + 10 + geo.shape self.geo = geo - self.obj = resources.flower_obj(frame) - + + try: + self.obj = resources.flower_obj(frame) + except: + # matplotlib failsafe + self.obj = u.gf_2d(u.parallel.MPInoise2d(frame),1.) + # Get probe - moon = resources.moon_pr(self.geo.shape) - moon /= np.sqrt(u.abs2(moon).sum() / r.photons) + try: + moon = resources.moon_pr(self.geo.shape) + except: + # matplotlib failsafe + moon = u.ellipsis(u.grids(self.geo.shape)).astype(complex) + + moon /= np.sqrt(u.abs2(moon).sum() / p.photons) self.pr = moon self.load_common_in_parallel = True - self.r = r + + self.p = p def load_positions(self): return self.pos @@ -1419,149 +1562,103 @@ def load(self, indices): s = self.geo.shape raw = {} + if self.p.add_poisson_noise: + logger.info("Generating data with poisson noise.") + else: + logger.info("Generating data without poisson noise.") + for k in indices: intensity_j = u.abs2(self.geo.propagator.fw( self.pr * self.obj[p[k][0]:p[k][0] + s[0], p[k][1]:p[k][1] + s[1]])) - if self.r.psf > 0.: - intensity_j = u.gf(intensity_j, self.r.psf) + if self.p.psf > 0.: + intensity_j = u.gf(intensity_j, self.p.psf) - raw[k] = np.random.poisson(intensity_j).astype(np.int32) + if self.p.add_poisson_noise: + raw[k] = np.random.poisson(intensity_j).astype(np.int32) + else: + raw[k] = intensity_j.astype(np.int32) return raw, {}, {} - -class DataSource(object): +@defaults_tree.parse_doc('scandata.QuickScan') +class QuickScan(PtyScan): """ - A source of data for ptychographic reconstructions. - - The important method is "feed_data", which returns packets of diffraction - patterns with their meta-data. + Test PtyScan to benchmark graph creation further down the line. + + Override parent class default: + + Defaults: + + [name] + default = MoonFlowerScan + type = str + help = + doc = + + [num_frames] + default = 100 + type = int + help = Number of frames to simulate + doc = + + [shape] + type = int, tuple + default = 64 + help = Shape of the region of interest cropped from the raw data. + doc = Cropping dimension of the diffraction frame + Can be None, (dimx, dimy), or dim. In the latter case shape will be (dim, dim). + userlevel = 1 + + [density] + default = 0.05 + type = float + help = Position distance in fraction of illumination frame """ - def __init__(self, scans, frames_per_call=10000000, feed_format='dp'): - """ - DataSource initialization. - - Parameters - ---------- - scans : - a dictionary of scan structures. - - frames_per_call : (optional) - number of frames to load in one call. - By default, load as many as possible. - feed_format : - the format in with the data is packaged. - For now only 'dp' is implemented. - """ - from ..experiment import PtyScanTypes - # FIXME: SC: when moved to top, import fails - - self.frames_per_call = frames_per_call - self.feed_format = feed_format - self.scans = scans - - # Sort after given keys - self.labels = sorted(scans.keys()) - - # Empty list for the scans - self.pty_scans = [] - - for label in self.labels: - # We are making a copy of the root as we want to fill it - scan = scans[label] - s = scan['pars'] - - # Copy other relevant information - prep = s.data.copy() - - # Assign label - prep.label = label - - # Assign source, recipe, and positions_theory - source = prep.source - recipe = prep.get('recipe', {}) - if prep.get('positions_theory') is None: - prep.positions_theory = scan['pos_theory'] - - # prep.dfile = s.data_file - # prep.geometry = s.geometry.copy() - # prep.xy = s.xy.copy() - - # if source is not None: - # source = source.lower() - if source is None or source.lower() == 'empty': - prep.recipe = None - logger.warning('Generating dummy PtyScan for scan `%s` - This ' - 'label will source only zeros as data' % label) - self.pty_scans.append(PtyScan(prep)) - elif source.lower() in PtyScanTypes: - pst = PtyScanTypes[source.lower()] - logger.info('Scan %s will be prepared with the recipe "%s"' - % (label, source)) - self.pty_scans.append(pst(prep, recipe=recipe)) - elif (source.endswith('.ptyd') or source.endswith('.pty') or - str(source) == 'file'): - self.pty_scans.append(PtydScan(prep, source=source)) - elif source.lower() == 'test': - self.pty_scans.append(MoonFlowerScan(prep)) - elif source.lower() == 'sim': - from ..simulations import SimScan - logger.info('Scan %s will simulated' % label) - self.pty_scans.append(SimScan(prep, s.copy())) - else: - raise RuntimeError('Could not manage source "%s" for scan `%s`.' - % (str(source), label)) - - # Initialize flags - self.scan_current = -1 - self.data_available = True - self.scan_total = len(self.pty_scans) - - @property - def scan_available(self): - return self.scan_current < (self.scan_total - 1) - - def feed_data(self): + def __init__(self, pars=None, **kwargs): """ - Yield data packages. + Parent pars are for the """ - # Get PtyScan instance to scan_number - cur_scan = self.pty_scans[self.scan_current] - label = self.labels[self.scan_current] - - # Initialize if that has not been done yet - if not cur_scan.is_initialized: - cur_scan.initialize() - msg = cur_scan.auto(self.frames_per_call, self.feed_format) + p = self.DEFAULT.copy(depth=99) + p.update(pars) - # If we catch a scan that has ended, look for an unfinished scan - while msg == EOS and self.scan_available: - self.scan_current += 1 - cur_scan = self.pty_scans[self.scan_current] - label = self.labels[self.scan_current] - - if not cur_scan.is_initialized: - cur_scan.initialize() + # Initialize parent class + super(QuickScan, self).__init__(p, **kwargs) - msg = cur_scan.auto(self.frames_per_call, self.feed_format) + # Derive geometry from input + geo = geometry.Geo(pars=self.meta) - self.data_available = (msg != EOS or self.scan_available) + # Derive scan pattern + pos = u.Param() + pos.spacing = geo.resolution * geo.shape * p.density + pos.steps = np.int(np.round(np.sqrt(self.num_frames))) + 1 + pos.extent = pos.steps * pos.spacing + pos.model = 'round' + pos.count = self.num_frames + self.pos = xy.from_pars(pos) + self.geo = geo + # dumm diff pattern + self.diff = np.zeros(geo.shape) + self.diff[0,0] = 1e7 + self.diff = np.fft.fftshift(self.diff) + + self.p = p - logger.debug(u.verbose.report(msg)) + def load_positions(self): + return self.pos - if msg != WAIT and msg != EOS: - # Ok that would be a data package - # Attach inner label - msg['common']['ptylabel'] = label - logger.info('Feeding data chunk') - return msg - else: - return None + def load_weight(self): + return np.ones(self.diff.shape) + def load(self, indices): + raw = {} + for k in indices: + raw[k] = self.diff.copy().astype(np.int32) + return raw, {}, {} + if __name__ == "__main__": u.verbose.set_level(3) MS = MoonFlowerScan(num_frames=100) diff --git a/ptypy/core/data_old.py b/ptypy/core/data_old.py deleted file mode 100644 index 630a4a3c6..000000000 --- a/ptypy/core/data_old.py +++ /dev/null @@ -1,1218 +0,0 @@ -""" -data - Diffraction data access - -This module defines a PtyScan, a container to hold the experimental -data of a ptychography scan. Instrument-specific reduction routines should -inherit PtyScan to prepare data for the Ptycho Instance in a uniform format. - -The experiment specific child class only needs to overwrite 2 functions -of the base class: - -For the moment the module contains two main objects: -PtyScan, which holds a single ptychography scan, and DataSource, which -holds a collection of datascans and feeds the data as required. - -This file is part of the PTYPY package. - - :copyright: Copyright 2014 by the PTYPY team, see AUTHORS. - :license: GPLv2, see LICENSE for details. -""" -import numpy as np -import os -import h5py -if __name__ == "__main__": - from ptypy import utils as u - from ptypy import io - from ptypy import resources - from ptypy.core import geometry - from ptypy.utils.verbose import logger, log, headerline -else: - from .. import utils as u - from .. import io - from .. import resources - from ..utils.verbose import logger, log, headerline - import geometry - -parallel = u.parallel - - -PTYD = dict( - chunks={}, # frames, positions - meta={}, # important to understand data. loaded by every process - info={}, # this dictionary is not loaded from a ptyd. Mainly for documentation - common={}, -) -""" Basic Structure of a .ptyd datafile """ - -META = dict( - label=None, # label will be set internally - experimentID=None, # a unique label of user choice - version='0.1', - shape=None, - psize=None, - #lam=None, - energy=None, - center=None, - distance = None, -) -GENERIC = dict( - dfile = None, # filename (e.g. 'foo.ptyd') - chunk_format='.chunk%02d', # Format for chunk file appendix. - #roi = None, # 2-tuple or int for the desired fina frame size - save = None, # None, 'merge', 'append', 'extlink' - auto_center = None, # False: no automatic center,None only if center is None, True it will be enforced - load_parallel = 'data', # None, 'data', 'common', 'all' - rebin = None, # rebin diffraction data - orientation = None, # None,int or 3-tuple switch, actions are (transpose, invert rows, invert cols) - min_frames = 1, # minimum number of frames of one chunk if not at end of scan - positions_theory = None, # Theoretical position list (This input parameter may get deprecated) - num_frames = None, # Total number of frames to be prepared - recipe = {}, -) -"""Default data parameters. See :py:data:`.scan.data` and a short listing below""" -GENERIC.update(META) - -WAIT = 'msg1' -EOS = 'msgEOS' -CODES = {WAIT: 'Scan unfinished. More frames available after a pause', - EOS: 'End of scan reached'} - -__all__ = ['GENERIC','PtyScan','PTYD','PtydScan','MoonFlowerScan'] - -import warnings -warnings.simplefilter('always', DeprecationWarning) -warnings.warn('This module is deprecated and will be removed from the package on 30/11/16',DeprecationWarning) - -class PtyScan(object): - """\ - PtyScan: A single ptychography scan, created on the fly or read from file. - BASECLASS - - Objectives: - - Stand alone functionality - - Can produce .ptyd data formats - - Child instances should be able to prepare from raw data - - On-the-fly support in form of chunked data. - - mpi capable, child classes should not worry about mpi - - """ - - DEFAULTS = GENERIC - - def __init__(self, pars=None, **kwargs): # filename='./foo.ptyd',shape=None, save=True): - """ - Class creation with minimum set of parameters, see DEFAULT dict in core/data.py - Please note the the class creation does not necessarily load data. - Call .initialize() to begin - """ - - # Load default parameter structure - info = u.Param(self.DEFAULTS.copy()) - info.update(pars) - info.update(kwargs) - - # validate(pars, '.scan.preparation') - - # Prepare meta data - self.meta = u.Param(META.copy()) - - # Attempt to get number of frames. - self.num_frames = info.num_frames - self.min_frames = info.min_frames * parallel.size - #logger.info('Looking for position information input parameter structure ....\n') - #if (info.positions_theory is None) and (info.xy is not None): - # from ptypy.core import xy - # info.positions_theory = xy.from_pars(info.xy) - - if info.positions_theory is not None: - num = len(info.positions_theory ) - logger.info('Theoretical positions are available. There will be %d frames.' % num) - logger.info('Any experimental position information will be ignored.') - logger.info('Former input value of frame number `num_frames` %s is overriden to %d' %(str(self.num_frames),num)) - self.num_frames = num - """ - # check if we got information on geometry from ptycho - if info.geometry is not None: - for k, v in info.geometry.items(): - # FIXME: This is a bit ugly - some parameters are added to info without documentation. - info[k] = v if info.get(k) is None else None - # FIXME: This should probably be done more transparently: it is not clear for the user that info.roi has precedence over geometry.N - if info.roi is None: - info.roi = u.expect2(info.geometry.N) - """ - # None for rebin should be allowed, as in "don't rebin". - if info.rebin is None: - info.rebin = 1 - - self.info = info - - # Print a report - log(4,'Ptypy Scan instance got the following parameters:') - log(4,u.verbose.report(info)) - - # Dump all input parameters as class attributes. - # FIXME: This duplication of parameters can lead to much confusion... - # self.__dict__.update(info) - - # Check MPI settings - lp = str(self.info.load_parallel) - self.load_common_in_parallel = (lp == 'all' or lp == 'common') - self.load_in_parallel = (lp == 'all' or lp == 'data') - - # set data chunk and frame counters to zero - self.framestart = 0 - self.chunknum = 0 - self.start = 0 - self.chunk = None - - # Initialize other instance attributes - self.common = {} - self.has_weight2d = None - self.has_positions = None - self.dfile = None - self.save = self.info.save - - # copy all values for meta - for k in self.meta.keys(): - self.meta[k] = self.info[k] - #self.center = None # Center will be set later - #self.roi = self.info.roi #None # ROI will be set later - #self.shape = None - self.orientation = self.info.orientation - self.rebin = self.info.rebin - - # initialize flags - self._flags = np.array([0, 0, 0], dtype=int) - self.is_initialized = False - - def initialize(self): - """ - Time for some read /write access - """ - - # Prepare writing to file - if self.info.save is not None: - # We will create a .ptyd - self.dfile = self.info.dfile - if parallel.master: - if os.path.exists(self.dfile): - backup = self.dfile + '.old' - logger.warning('File %s already exist. Renamed to %s' % (self.dfile, backup)) - os.rename(self.dfile, backup) - # Prepare an empty file with the appropriate structure - io.h5write(self.dfile, PTYD.copy()) - # Wait for master - parallel.barrier() - - if parallel.master or self.load_common_in_parallel: - self.common = self.load_common() - # FIXME - # Replace Nones, because we cannot broadcast them later - # I don't get it, why didn't we allow for a missing key? - # Also I disagree that we should enforce a 'weight2d' or - # positions. The user can very well add them later. - self.common = dict([(k,v) if v is not None else (k,np.array([])) for k,v in self.common.items() ]) - - # broadcast - if not self.load_common_in_parallel: - parallel.bcast_dict(self.common) - - self.common = u.Param(self.common) - assert 'weight2d' in self.common and 'positions_scan' in self.common - - logger.info('\n'+headerline('Analysis of the "common" arrays','l')) - # Check if weights (or mask) have been loaded by load_common. - weight2d = self.common.weight2d - self.has_weight2d = weight2d is not None and len(weight2d)>0 - logger.info('Check for weight or mask, "weight2d" .... %s : shape = %s' % (str(self.has_weight2d),str(weight2d.shape))) - - - # Check if positions have been loaded by load_common - positions = self.common.positions_scan - self.has_positions = positions is not None and len(positions)>0 - logger.info('Check for positions, "positions_scan" .... %s : shape = %s' % (str(self.has_positions),str(positions.shape))) - - - if self.info.positions_theory is not None: - logger.info('Skipping experimental positions `positions_scan`') - elif self.has_positions: - # Store positions in the info dictionary - self.info.positions_scan = positions - num_pos = len(positions) - if self.num_frames is None: - # Frame number was not known. We just set it now. - logger.info('Scanning positions found. There will be %d frames' % num_pos) - self.num_frames = num_pos - else: - # Frame number was already specified. Maybe we didn't want to use everything? - if num_pos > self.num_frames: - #logger.info('Scanning positions have the same number of points as the theoretical ones (%d).' % num_pos) - logger.info('Scanning positions (%d) exceed the desired number of scan points (%d).' % (num_pos,self.num_frames)) - logger.info('Set `num_frames` to None or to a larger value for more scan points') - elif num_pos < self.num_frames: - logger.info('Scanning positions (%d) are fewer than the desired number of scan points (%d).' % (num_pos,self.num_frames)) - logger.info('Resetting `num_frames` to lower value') - self.num_frames = num_pos - #raise RuntimeError('Scanning positions have a number of points (%d) inconsistent with what was previously deduced (%d).' - # % (num_pos, self.info.num_frames)) - else: - logger.info('No scanning position have been provided at this stage.') - - if self.num_frames is None: - logger.warning('Number of frames `num_frames` not specified at this stage\n.') - - parallel.barrier() - - #logger.info('####### MPI Report: ########\n') - log(4,u.verbose.report(self.common),True) - parallel.barrier() - logger.info(headerline('Analysis done','l')+'\n') - - if self.info.save is not None and parallel.master: - logger.info('Appending common dict to file %s\n' % self.info.dfile) - io.h5append(self.info.dfile, common=dict(self.common), info=dict(self.info)) - # wait for master - parallel.barrier() - - self.is_initialized = True - - def _finalize(self): - """ - Last actions when Eon-of-Scan is reached - """ - # maybe do this at end of everything - if self.info.save is not None and parallel.master: - io.h5append(self.info.dfile, common=dict(self.common), info=dict(self.info)) - - def load_common(self): - """ - **Overwrite in child class** - - Loads arrays that are common and needed for preparation of all - other data frames coming in. Any array loaded here and returned - in a dict will be distributed afterwards through a broadcoast. - It is not intended to load diffraction data in this step. - - The main purpose is that there may be common data to all processes, that - is slow to retrieve or large files, such that if all processes - attempt to get the data, perfomance will decline or RAM usage - may be too large. - - It is a good idea to load dark field, flat field, etc - Also positions may be handy to load here - - If `load_parallel` is set to `all` or common`, this function is - executed by all nodes, otherwise the master node executes this - function and braodcasts the results to other nodes. - - Returns - ------- - common : dict of numpy arrays - At least two keys, `weight2d` and `positions_scan` must be - given (they can be None). The return dictionary is available - throught - - Note - ---- - The return signature of this function is not yet fixed and may - get altered in near future. One Option would be to include - `weight2d` and `positions_scan` as part of the return signature - and thus fixing them in the Base class. - """ - weight2d = None if self.info.shape is None else np.ones(u.expect2(self.info.shape), dtype='bool') - positions_scan = None if self.num_frames is None else np.indices((self.num_frames, 2)).sum(0) - return {'weight2d': weight2d, - 'positions_scan': positions_scan} - - def _mpi_check(self, chunksize, start=None): - """ - Executes the check() function on master node and communicates - the result with the other nodes. - This function determines if the end of the scan is reached - or if there is more data after a pause. - - returns: - - codes WAIT or EOS - - or (start, frames) if data can be loaded - """ - # take internal counter if not specified - s = self.framestart if start is None else int(start) - - # check may contain data system access, so we keep it to one process - if parallel.master: - self.frames_accessible, eos = self.check(chunksize, start=s) - if self.num_frames is None and eos is None: - logger.warning('Number of frames not specified and .check() cannot determine end-of-scan. Aborting..') - self.abort = True - self.end_of_scan = (s + self.frames_accessible >= self.num_frames) if eos is None else eos - # wait for master - parallel.barrier() - # communicate result - self._flags = parallel.bcast(self._flags) - - # abort here if the flag was set - if self.abort: - raise RuntimeError('Load routine incapable to determine end-of-scan.') - - N = self.frames_accessible - # wait if too few frames are available and we are not at the end - if N < self.min_frames and not self.end_of_scan: - return WAIT - elif self.end_of_scan and N <= 0: - return EOS - else: - # move forward,set new starting point - self.framestart += N - return s, N - - def _mpi_indices(self, start, step): - """ - Funtion to generate the diffraction data index lists that - determine which node contains which data. - """ - indices = u.Param() - - # all indices in this chunk of data - indices.chunk = range(start, start + step) - - # let parallel.loadmanager take care of assigning these indices to nodes - indices.lm = parallel.loadmanager.assign(indices.chunk) - # this one contains now a list of indices listed after ranke - - # index list (node specific) - indices.node = [indices.chunk[i] for i in indices.lm[parallel.rank]] - - # store internally - self.indices = indices - return indices - - def get_data_chunk(self, chunksize, start=None): - """ - This function prepares a container that is compatible to data package - This function is called from the auto() function. - """ - msg = self._mpi_check(chunksize, start) - if msg in [EOS, WAIT]: - logger.info(CODES[msg]) - return msg - else: - start, step = msg - - # get scan point index distribution - indices = self._mpi_indices(start, step) - - # what did we get? - data, positions, weights = self._mpi_pipeline_with_dictionaries(indices) - # all these dictionaries could be empty - # fill weights dictionary with references to the weights in common - - has_data = (len(data) > 0) - has_weights = (len(weights) > 0) and len(weights.values()[0])>0 - - if has_data: - dsh = np.array(data.values()[0].shape[-2:]) - else: - dsh = np.zeros([0, 0]) - - # communicate result - parallel.MPImax(dsh) - - if not has_weights: - # peak at first item - if self.has_weight2d: - altweight = self.common.weight2d - else: - try: - altweight = self.meta.weight2d - except: - altweight = np.ones(dsh) - weights = dict.fromkeys(data.keys(), altweight) - - assert len(weights) == len(data), 'Data and Weight frames unbalanced %d vs %d' % (len(data), len(weights)) - - sh = self.info.shape - # adapt roi if not set - if sh is None: - logger.info('ROI not set. Using full frame shape of (%d,%d).' % tuple(dsh)) - sh = dsh - else: - sh = u.expect2(sh) - - # only allow square slices in data - if sh[0] != sh[1]: - roi = u.expect2(sh.min()) - logger.warning('Asymmetric data ROI not allowed. Setting ROI from (%d,%d) to (%d,%d)' % (sh[0],sh[1],roi[0],roi[1])) - sh = roi - - self.info.shape = sh - - cen = self.info.center - if str(cen)==cen: - cen=geometry.translate_to_pix(sh,cen) - - auto = self.info.auto_center - # get center in diffraction image - if auto is None or auto is True: - auto_cen = self._mpi_autocenter(data, weights) - else: - auto_cen = None - - if cen is None and auto_cen is not None: - logger.info('Setting center for ROI from %s to %s.' %(str(cen),str(auto_cen))) - cen = auto_cen - elif cen is None and auto is False: - cen = dsh // 2 - else: - # center is number or tuple - cen = u.expect2(cen[-2:]) - if auto_cen is not None: - logger.info('ROI center is %s, automatic guess is %s.' %(str(cen),str(auto_cen))) - - # It is important to set center again in order to NOT have different centers for each chunk - # the downside is that the center may be off if only a few diffraction patterns were - # used for the analysis. In such case it is beneficial to set the center in the parameters - self.info.center = cen - - # make sure center is in the image frame - assert (cen > 0).all() and ( - dsh - cen > 0).all(), 'Optical axes (center = (%.1f,%.1f) outside diffraction image frame (%d,%d)' % tuple(cen)+tuple(dsh) - - # determine if the arrays require further processing - do_flip = self.orientation is not None and np.array(self.orientation).any() - do_crop = (np.abs(sh - dsh) > 0.5).any() - do_rebin = self.rebin is not None and (self.rebin != 1) - - if do_flip or do_crop or do_rebin: - logger.info('Enter preprocessing (crop/pad %s, rebin %s, flip/rotate %s) ... \n' % - (str(do_crop), str(do_flip), str(do_rebin))) - # we proceed with numpy arrays. That is probably now - # more memory intensive but shorter in writing - if has_data: - d = np.array([data[ind] for ind in indices.node]) - w = np.array([weights[ind] for ind in indices.node]) - else: - d = np.ones((1,)+tuple(dsh)) - w = np.ones((1,)+tuple(dsh)) - - # flip, rotate etc. - d, tmp = u.switch_orientation(d, self.orientation, cen) - w, cen = u.switch_orientation(w, self.orientation, cen) - - # crop - d, tmp = u.crop_pad_symmetric_2d(d, sh, cen) - w, cen = u.crop_pad_symmetric_2d(w, sh, cen) - - # rebin, check if rebinning is neither to strong nor impossible - rebin = self.rebin - if rebin<=1: - pass - elif rebin in range(2, 6) and (((self.roi / float(rebin)) % 1) == 0.0).all(): - mask = w > 0 - d = u.rebin_2d(d, rebin) - w = u.rebin_2d(w, rebin) - mask = u.rebin_2d(mask, rebin) - w[mask < 1] = 0 - else: - raise RuntimeError('Binning (%d) is to large or incompatible with array shape (%s)' % (rebin,str(tuple(sh)))) - - if has_data: - # translate back to dictionaries - data = dict(zip(indices.node, d)) - weights = dict(zip(indices.node, w)) - - # adapt meta info - self.meta.center = cen / float(self.rebin) - self.meta.shape = u.expect2(sh) / self.rebin - self.meta.psize = u.expect2(self.info.psize) * self.rebin if self.info.psize is not None else None - - # prepare chunk of data - chunk = u.Param() - chunk.indices = indices.chunk - chunk.indices_node = indices.node - chunk.num = self.chunknum - chunk.data = data - - # if there were weights we add them to chunk, - # otherwise we push it into meta - if has_weights: - chunk.weights = weights - elif has_data: - chunk.weights = {} - self.meta.weight2d = weights.values()[0] - - # slice positions from common if they are empty too - if positions is None or len(positions) == 0: - pt = self.info.positions_theory - if pt is not None: - chunk.positions = pt[indices.chunk] - else: - try: - chunk.positions = self.info.positions_scan[indices.chunk] - except: - logger.info('Unable to slice position information from experimental or theoretical ressource') - chunk.positions = [None]*len(indices.chunk) - else: - # a dict : sort positions to indices.chunk - # this may fail if there a less positions than scan points (unlikely) - chunk.positions = np.asarray([positions[i] for i in indices.chunk]) - # positions complete - - # with first chunk we update meta - if self.chunknum < 1: - """ - for k, v in self.meta.items(): - # FIXME: I would like to avoid this "blind copying" - # BE: This is not a blind copy as only keys in META above are used - self.meta[k] = self.__dict__.get(k, self.info.get(k)) if v is None else v - self.meta['center'] = cen - """ - - if self.info.save is not None and parallel.master: - io.h5append(self.dfile, meta=dict(self.meta)) - parallel.barrier() - - self.chunk = chunk - self.chunknum += 1 - - return chunk - - def auto(self, frames, chunk_form='dp'): - """ - Repeated calls to this function will process the data - - Parameters - ---------- - frames : int - Number of frames to process. - - Returns - ------- - variable - one of the following - - None, if scan's end is not reached, but no data could be prepared yet - - False, if scan's end is reached - - a data package otherwise - """ - # attempt to get data: - msg = self.get_data_chunk(frames) - if msg == WAIT: - return msg - elif msg == EOS: - # cleanup maybe? - self._finalize() - # del self.common - # del self.chunk - return msg - else: - out = self.return_chunk_as(msg, chunk_form) - # save chunk - if self.info.save is not None: - self._mpi_save_chunk(self.info.save,msg) - # delete chunk - del self.chunk - return out - - def return_chunk_as(self, chunk, kind='dp'): - """ - Returns the loaded data chunk `chunk` in the format `kind` - For now only kind=='dp' (data package) is valid. - """ - # this is a bit ugly now - if kind != 'dp': - raise RuntimeError('Unknown kind of chunck format: %s' % str(kind)) - - out = {} - - # The "common" part - out['common'] = self.meta - - # The "iterable" part - iterables = [] - for pos, index in zip(chunk.positions, chunk.indices): - frame = {'index': index, - 'data': chunk.data.get(index), - 'position': pos} - if frame['data'] is None: - frame['mask'] = None - else: - # ok we now know that we need a mask since data is not None - # first look in chunk for a weight to this index, then - # look for a 2d-weight in meta, then arbitrarily set - # weight to ones. - w = chunk.weights.get(index, self.meta.get('weight2d', np.ones_like(frame['data']))) - frame['mask'] = (w > 0) - iterables.append(frame) - out['iterable'] = iterables - return out - - def _mpi_pipeline_with_dictionaries(self, indices): - """ - example processing pipeline using dictionaries. - - return : - positions, data, weights - -- Dictionaries. Keys are the respective scan point indices - `positions` and `weights` may be empty. If so, the information - is taken from the self.common dictionary - - """ - if self.load_in_parallel: - # all nodes load raw_data and slice according to indices - raw, pos, weights = self.load(indices=indices.node) - - # gather postion information as every node needs it later - pos = parallel.gather_dict(pos) - else: - if parallel.master: - raw, pos, weights = self.load(indices=indices.chunk) - else: - raw = {} - pos = {} - weights = {} - # distribute raw data across nodes according to indices - raw = parallel.bcast_dict(raw, indices.node) - weights = parallel.bcast_dict(weights, indices.node) - - # (re)distribute position information - every node should now be - # aware of all positions - parallel.bcast_dict(pos) - - # prepare data across nodes - data, weights = self.correct(raw, weights, self.common) - - return data, pos, weights - - def check(self, frames=None, start=None): - """ - **Overwrite in child class** - - This method checks how many frames the preparation routine may - process, starting from frame `start` at a request of `frames_requested`. - - This method is supposed to return the number of accessible frames - for preparation and should determine if data acquistion for this - scan is finished. - - Parameters - ---------- - frames : int or None - Number of frames requested. - start : int or None - Scanpoint index to start checking from. - - Returns - ------- - frames_accessible : int - Number of frames readable. - end_of_scan : int or None - is one of the following: - - 0, end of the scan is not reached - - 1, end of scan will be reached or is - - None, can't say - - """ - if start is None: - start = self.framestart - if frames is None: - frames = self.min_frames - - frames_accessible = min((frames, self.num_frames - start)) - - return frames_accessible, None - - @property - def end_of_scan(self): - return not (self._flags[1] == 0) - - @end_of_scan.setter - def end_of_scan(self, eos): - self._flags[1] = int(eos) - - @property - def frames_accessible(self): - return self._flags[0] - - @frames_accessible.setter - def frames_accessible(self, frames): - self._flags[0] = frames - - @property - def abort(self): - return not (self._flags[2] == 0) - - @abort.setter - def abort(self, abort): - self._flags[2] = int(abort) - - def load(self, indices): - """ - **Overwrite in child class** - - Loads data according to node specific scanpoint indeces that have - been determined by the loadmanager from utils.parallel or otherwise - - Returns - ------- - raw, positions, weight : dict - Dictionaries whose keys are the given scan point `indices` - and whose values are the respective frame / position according - to the scan point index. `weight` and `positions` may be empty - """ - # dummy fill - raw = dict((i, i * np.ones(u.expect2(self.info.shape))) for i in indices) - return raw, {}, {} - - def correct(self, raw, weights, common): - """ - **Overwrite in child class** - - Place holder for dark and flatfield correction. If :any:`load` - already provides data in the form of photon counts, and no frame - specific weight is needed, this method may be left as is - - May get *merged* with :any:`load` in future. - - Returns - ------- - data,weights : dict - Flat and dark-corrected data dictionaries. These dictionaries - must have the same keys as the input `raw` and contain - corrected frames (`data`) and statistical weights (`weights`) - which are zero for invalid or masked pixel other the number - of detector counts that correspond to one photon count - """ - # c=dict(indices=None,data=None,weight=None) - data = raw - return data, weights - - def _mpi_autocenter(self, data, weights): - """ - Calculates the frame center across all nodes. - Data and weights are dicts of the same length and different on each node - """ - cen = dict([(k, u.mass_center(d * (weights[k] > 0))) for k, d in data.iteritems()]) - # for some nodes, cen may still be empty. Therefore we use gather_dict to be save - cen = parallel.gather_dict(cen) - parallel.barrier() - # now master possesses all calculated centers - if parallel.master: - cen = np.array(cen.values()).mean(0) - else: - cen = np.array([0., 0.]) - # print cen - parallel.allreduce(cen) - return cen - - def report(self, what=None, shout=True): - """ - Make a report on internal structure - """ - what = what if what is not None else self.__dict__ - msg = u.verbose.report(what) - if shout: - logger.info(msg, extra={'allprocesses': True}) - else: - return msg - - def _mpi_save_chunk(self, kind='link', chunk=None): - """ - Saves data chunk to hdf5 file specified with `dfile` - It works by gathering weights and data to the master node. - Master node then writes to disk. - - In case you support parallel hdf5 writing, please modify this - function to suit your installation. - - 2 out of 3 modes currently supported - - kind : 'merge','append','link' - - 'append' : appends chunks of data in same file - 'link' : saves chunks in seperate files and adds ExternalLinks - - TODO: - * For the 'link case, saving still requires access to - main file `dfile` even so for just adding the link. - This may result in conflict if main file is polled often - by seperate read process. - Workaraound would be to write the links on startup in - initialise() - - """ - # gather all distributed dictionary data. - c = chunk if chunk is not None else self.chunk - - # shallow copy - todisk = dict(c) - num = todisk.pop('num') - ind = todisk.pop('indices_node') - for k in ['data', 'weights']: - if k in c.keys(): - if hasattr(c[k], 'iteritems'): - v = c[k] - else: - v = dict(zip(ind, np.asarray(c[k]))) - parallel.barrier() - # gather the content - newv = parallel.gather_dict(v) - todisk[k] = np.asarray([newv[i] for i in sorted(newv.keys())]) - - parallel.barrier() - # all information is at master node. - if parallel.master: - - # form a dictionary - if str(kind) == 'append': - h5address = 'chunks/%d' % num - io.h5append(self.dfile, {h5address: todisk}) - - elif str(kind) == 'link': - import h5py - h5address = 'chunks/%d' % num - hddaddress = self.dfile + '.part%03d' % num - with h5py.File(self.dfile) as f: - f[h5address] = h5py.ExternalLink(hddaddress, '/') - f.close() - io.h5write(hddaddress, todisk) - elif str(kind) == 'merge': - raise NotImplementedError('Merging all data into single chunk is not yet implemented') - parallel.barrier() - -class PtydScan(PtyScan): - """ - PtyScan provided by native "ptyd" file format. - """ - DEFAULT = GENERIC.copy() - - def __init__(self, pars=None, source=None,**kwargs): - """ - PtyScan provided by native "ptyd" file format. - - :param source: Explicit source file. If not None or 'file', - the data may get processed depending on user input - - :param pars: Input like PtyScan - """ - # create parameter set - p = u.Param(self.DEFAULT.copy()) - - # copy the label - #if pars is not None: - # p.label = pars.get('label') - - if source is None or str(source)=='file': - # this is the case of absolutely no additional work - logger.info('No explicit source file was given. Will continue read only') - source = pars['dfile'] - manipulate = False - elif pars is None or len(pars)==0: - logger.info('No parameters provided. Saving / modification disabled') - manipulate = False - else: - logger.info('Explicit source file given. Modification is possible\n') - dfile = pars['dfile'] - # check for conflict - if str(u.unique_path(source))==str(u.unique_path(dfile)): - logger.info('Source and Sink files are the same.') - dfile = os.path.splitext(dfile) - dfile = dfile[0] +'_n.'+ dfile[1] - logger.info('Will instead save to %s if necessary.' % os.path.split(dfile)[1]) - - pars['dfile']= dfile - manipulate = True - p.update(pars) - - - # make sure the source exists. - assert os.path.exists(u.unique_path(source)), 'Source File (%s) not found' % source - - self.source = source - - # get meta information - meta = u.Param(io.h5read(self.source, 'meta')['meta']) - - # update given parameters when they are None - if not manipulate: - super(PtydScan, self).__init__(meta, **kwargs) - else: - # overwrite only those set to None - for k, v in meta.items(): - if p.get(k) is None: - p[k] = v - # Initialize parent class and fill self - super(PtydScan, self).__init__(p, **kwargs) - - # enforce that orientations are correct - # Other instance attributes - self._checked = {} - self._ch_frame_ind = None - - - def check(self, frames=None, start=None): - """ - Implementation of the check routine for a .ptyd file format - - See also - -------- - Ptyscan.check - """ - - if start is None: - start = self.framestart - if frames is None: - frames = self.min_frames - - # Get info about size of currently available chunks. - # Dead external links will produce None and are excluded - with h5py.File(self.source, 'r') as f: - d = {} - chitems = sorted([(int(k), v) for k, v in f['chunks'].iteritems() if v is not None], key=lambda t: t[0]) - for chkey in chitems[0][1].keys(): - d[chkey] = np.array([((int(k),) + v[chkey].shape) for k, v in chitems if v is not None]) - f.close() - - self._checked = d - allframes = int(sum([ch[1] for ch in d['data']])) - self._ch_frame_ind = np.array([(dd[0], frame) for dd in d['data'] for frame in range(dd[1])]) - - return min((frames, allframes - start)), allframes < start+frames - - def _coord_to_h5_calls(self, key, coord): - return 'chunks/%d/%s' % (coord[0], key), slice(coord[1], coord[1] + 1) - - def load_common(self): - """ - In ptyd, 'common' must exist - """ - # this total buggy right now - return {'weight2d' : self.info.weight2d, 'positions_scan' : None} - - def load(self, indices): - """ - Load from ptyd. Due to possible chunked data, slicing frames is - non-trivial - """ - # ok we need to communicate the some internal info - parallel.barrier() - self._ch_frame_ind=parallel.bcast(self._ch_frame_ind) - parallel.barrier() - parallel.bcast_dict(self._checked) - - # get the coordinates in the chunks - coords = self._ch_frame_ind[indices] - calls = {} - for key in self._checked.keys(): - calls[key] = [self._coord_to_h5_calls(key, c) for c in coords] - - # get our data from the ptyd file - out = {} - with h5py.File(self.source, 'r') as f: - for array, call in calls.iteritems(): - out[array] = [np.squeeze(f[path][slce]) for path, slce in call] - f.close() - - # if the chunk provided indices, we use those instead of our own - # Dangerous and not yet implemented - # indices = out.get('indices',indices) - - # wrap in a dict - for k, v in out.iteritems(): - out[k] = dict(zip(indices, v)) - - return (out.get(key, {}) for key in ['data', 'positions', 'weights']) - - -class MoonFlowerScan(PtyScan): - """ - Test Ptyscan class producing a romantic ptychographic dataset of a moon - illuminating flowers. - """ - - DEFAULT = GENERIC.copy().update(geometry.DEFAULT.copy()) - - def __init__(self, pars = None, **kwargs): - """ - Parent pars are for the - """ - p = geometry.DEFAULT.copy() - if pars is not None: - p.update(pars) - # Initialize parent class - super(MoonFlowerScan, self).__init__(p, **kwargs) - - from ptypy import resources - from ptypy.core import xy - - # derive geometry from input - G = geometry.Geo(pars = self.meta) - #G._initialize(self.meta) - - # derive scan pattern - pos = u.Param() - pos.spacing = G.resolution * G.shape / 5. - pos.layers = np.int(np.round(np.sqrt(self.num_frames)))+1 - pos.extent = pos.layers * pos.spacing - pos.model = 'round' - pos.count = self.num_frames - self.pos = xy.from_pars(pos) - - # calculate pixel positions - pixel = self.pos / G.resolution - pixel -= pixel.min(0) - self.pixel = np.round(pixel).astype(int) + 10 - frame = self.pixel.max(0) + 10 + G.shape - self.G = G - #from matplotlib import pyplot as plt - #plt.figure(200);plt.imshow(u.imsave(G.propagator.pre_ifft)) - #plt.figure(101);plt.imshow(G.propagator.grids_det[0]**2+G.propagator.grids_det[1]**2) - # get object - self.obj = resources.flower_obj(frame) - - # get probe - moon = resources.moon_pr(self.G.shape) - moon /= np.sqrt(u.abs2(moon).sum() / 1e8) - self.pr = moon - self.load_common_in_parallel = True - - def load_common(self): - """ - Transmit positions - """ - return {'weight2d': np.ones(self.pr.shape), - 'positions_scan': self.pos} - - def load(self, indices): - """ - Forward propagation - """ - # dummy fill - p=self.pixel - s=self.G.shape - raw = {} - for i in indices: - raw[i]=np.random.poisson(u.abs2(self.G.propagator.fw(self.pr * self.obj[p[i][0]:p[i][0]+s[0],p[i][1]:p[i][1]+s[1]]))).astype(np.int32) - return raw, {}, {} - - - -class DataSource(object): - """ - A source of data for ptychographic reconstructions. The important method is "feed_data", which returns - packets of diffraction patterns with their meta-data. - """ - def __init__(self, scans, frames_per_call=10000000, feed_format='dp'): - """ - DataSource initialization. - - scans: a dictionnary of scan structures. - frames_per_call: (optional) number of frames to load in one call. By default, load as many as possible. - feed_format: the format in with the data is packaged. For now only 'dp' is implemented. - """ - - from ..experiment import PtyScanTypes - - self.frames_per_call = frames_per_call - self.feed_format = feed_format - self.scans = scans - - # sort after keys given - self.labels = sorted(scans.keys()) - - # empty list for the scans - self.PS = [] - - for label in self.labels: - # we are making a copy of the root as we want to fill it - scan = scans[label] - s = scan['pars'] - - # Copy other relevant information - prep = s.data.copy() - - # Assign label - prep.label = label - - source = prep.source - recipe = prep.get('recipe',{}) - if prep.get('positions_theory') is None: - prep.positions_theory = scan['pos_theory'] - - prep.dfile = s.data_file - #prep.geometry = s.geometry.copy() - #prep.xy = s.xy.copy() - - if source is not None: - source = source.lower() - - if source in PtyScanTypes: - PS = PtyScanTypes[source] - logger.info('Scan %s will be prepared with the recipe "%s"' % (label, source)) - self.PS.append(PS(prep, recipe= recipe)) - elif source.endswith('.ptyd') or source.endswith('.pty') or str(source)=='file': - self.PS.append(PtydScan(prep, source=source)) - elif source=='test': - self.PS.append(MoonFlowerScan(prep)) - elif source=='sim': - from ..simulations import SimScan - logger.info('Scan %s will simulated' % (label)) - self.PS.append(SimScan(prep,s.copy())) - elif source=='empty' or source is None: - prep.recipe = None - logger.warning('Generating dummy PtyScan for scan `%s` - This label will source only zeros as data' % label) - self.PS.append(PtyScan(prep)) - else: - raise RuntimeError('Could not manage source "%s" for scan `%s`' % (str(source),label)) - - # Initialize flags - self.scan_current = -1 - self.data_available = True - self.scan_total = len(self.PS) - - @property - def scan_available(self): - return self.scan_current < (self.scan_total - 1) - - def feed_data(self): - """ - Yield data packages. - """ - # get PtyScan instance to scan_number - PS = self.PS[self.scan_current] - label = self.labels[self.scan_current] - - # initialize if that has not been done yet - if not PS.is_initialized: - PS.initialize() - - msg = PS.auto(self.frames_per_call, self.feed_format) - # if we catch a scan that has ended look for an unfinished scan - while msg == EOS and self.scan_available: - self.scan_current += 1 - PS = self.PS[self.scan_current] - label = self.labels[self.scan_current] - if not PS.is_initialized: - PS.initialize() - msg = PS.auto(self.frames_per_call, self.feed_format) - - self.data_available = (msg != EOS or self.scan_available) - - logger.debug(u.verbose.report(msg)) - if msg != WAIT and msg != EOS: - # ok that would be a data package - # attach inner label - msg['common']['ptylabel'] = label - #print u.verbose.report(msg) - logger.info('Feeding data chunk') - return msg - else: - return None - -if __name__ == "__main__": - u.verbose.set_level(3) - MS = MoonFlowerScan(num_frames=100) - MS.initialize() - for i in range(50): - msg = MS.auto(10) - logger.info(u.verbose.report(msg), extra={'allprocesses': True}) - parallel.barrier() - diff --git a/ptypy/core/geometry.py b/ptypy/core/geometry.py index f751dffca..4c0ade50b 100644 --- a/ptypy/core/geometry.py +++ b/ptypy/core/geometry.py @@ -6,51 +6,23 @@ :copyright: Copyright 2014 by the PTYPY team, see AUTHORS. :license: GPLv2, see LICENSE for details. """ -# for solo use ########## -if __name__ == "__main__": - from ptypy import utils as u - from ptypy.utils.verbose import logger - from ptypy.core import Base - GEO_PREFIX = 'G' -# for in package use ##### -else: - from .. import utils as u - from ..utils.verbose import logger - from classes import Base,GEO_PREFIX - import numpy as np from scipy import fftpack + +from .. import utils as u +from ..utils.verbose import logger +from .classes import Base, GEO_PREFIX +from ..utils.descriptor import EvalDescriptor + try: import pyfftw import pyfftw.interfaces.numpy_fft as fftw_np except ImportError: - print "Warning: unable to import pyFFTW! Will use a slower FFT method." + pass + #logger.warning("Unable to import pyFFTW! Will use a slower FFT method.") -__all__ = ['DEFAULT', 'Geo', 'BasicNearfieldPropagator', - 'BasicFarfieldPropagator'] +__all__ = ['Geo', 'BasicNearfieldPropagator', 'BasicFarfieldPropagator'] -DEFAULT = u.Param( - # Incident photon energy (in keV) - energy=6.2, - # Wavelength (in meters) - lam=None, - # Distance from object to screen - distance=7.0, - # Pixel size (in meters) at detector plane - psize=172e-6, - # Pixel size (in meters) at sample plane - resolution=None, - # Number of detector pixels - shape=256, - # Propagation type - propagation='farfield', - misfit=0, - center='fftshift', - origin='fftshift', - precedence=None, -) -""" Default geometry parameters. See also :py:data:`.scan.geometry` - and an unflattened representation below """ _old2new = u.Param( # Distance from object to screen @@ -67,6 +39,8 @@ ) +local_tree = EvalDescriptor('') +@local_tree.parse_doc() class Geo(Base): """ Hold and keep consistent the information about experimental parameters. @@ -80,11 +54,89 @@ class Geo(Base): interact : bool (True) If set to True, changes to properties like :py:meth:`energy`, :py:meth:`lam`, :py:meth:`shape` or :py:meth:`psize` will cause - a call to :py:meth:`update` - + a call to :py:meth:`update`. + + + Default geometry parameters. See also :py:data:`.scan.geometry` + + Defaults: + + [energy] + type = float + default = 6.2 + help = Energy (in keV) + doc = If ``None``, uses `lam` instead. + userlevel = 0 + lowlim = 0 + + [lam] + type = float + default = None + help = Wavelength (in meters) + doc = Used only if `energy` is ``None`` + userlevel = 0 + lowlim = 0 + + [distance] + type = float + default = 7.19 + help = Distance from object to detector (in meters) + doc = + userlevel = 0 + lowlim = 0 + + [psize] + type = float + default = 0.000172 + help = Pixel size in the detector plane (in meters) + doc = + userlevel = 1 + lowlim = 0 + + [resolution] + type = float + default = None + help = Pixel size in the sample plane + doc = This parameter is used only for simulations + userlevel = 2 + lowlim = 0 + + [propagation] + type = str + default = farfield + help = Propagation type + doc = Either "farfield" or "nearfield" + userlevel = 1 + + [shape] + type = int, tuple + default = 256 + help = Number of pixels in detector frame + doc = Can be a 2-tuple of int (Nx, Ny) or an int N, in which case it is interpreted as (N, N). + userlevel = 1 + + [misfit] + type = float, tuple + default = 0. + help = TODO: WHAT IS MISFIT? + doc = + userlevel = 2 + + [center] + type = int, tuple, str + default = 'fftshift' + help = TODO: document this parameter. + doc = + userlevel = 2 + + [origin] + type = int, tuple, str + default = 'fftshift' + help = TODO: document this parameter. + doc = + userlevel = 2 """ - DEFAULT = DEFAULT _keV2m = 1.23984193e-09 _PREFIX = GEO_PREFIX @@ -105,13 +157,9 @@ def __init__(self, owner=None, ID=None, pars=None, **kwargs): if the key exists in `Geo.DEFAULT`. """ super(Geo, self).__init__(owner, ID) - # if len(kwargs)>0: - # self._initialize(**kwargs) - - # def _initialize(self, pars=None, **kwargs): # Starting parameters - p = u.Param(DEFAULT) + p = self.DEFAULT.copy(99) if pars is not None: p.update(pars) for k, v in p.items(): @@ -122,6 +170,13 @@ def __init__(self, owner=None, ID=None, pars=None, **kwargs): p[k] = v self.p = p + self._initialize(p) + + def _initialize(self, p): + """ + Parse input parameters, fill missing parameters and set up a + propagator. + """ self.interact = False # Set distance @@ -410,8 +465,8 @@ def assign_fft(self): elif str(self.ffttype) == 'numpy': self._numpy_fft() elif isinstance(self.ffttype, tuple): - self.fft = ffttype[0] - self.ifft = ffttype[1] + self.fft = self.ffttype[0] + self.ifft = self.ffttype[1] return (self.fft, self.ifft) @@ -425,7 +480,6 @@ class BasicFarfieldPropagator(object): Be aware though, that if the origin is not in the center of the frame, coordinates are rolled periodically, just like in the conventional fft case. """ - DEFAULT = DEFAULT def __init__(self, geo_pars=None, ffttype='fftw', **kwargs): """ @@ -456,7 +510,7 @@ def __init__(self, geo_pars=None, ffttype='fftw', **kwargs): self.post_ifft = None # Get default parameters and update - self.p = u.Param(DEFAULT) + self.p = u.Param(Geo.DEFAULT) if 'dtype' in kwargs: self.dtype = kwargs['dtype'] else: @@ -495,12 +549,12 @@ def update(self, geo_pars=None, **kwargs): lz /= (self.sh[0] + mis[0] - self.crop_pad[0]) / self.sh[0] # Calculate the grids - if str(p.origin) == p.origin: + if u.isstr(p.origin): c_sam = p.origin else: c_sam = p.origin + self.crop_pad / 2. - if str(p.center) == p.center: + if u.isstr(p.center): c_det = p.center else: c_det = p.center + self.crop_pad / 2. @@ -604,7 +658,6 @@ class BasicNearfieldPropagator(object): """ Basic two step (i.e. two ffts) Nearfield Propagator. """ - DEFAULT = DEFAULT def __init__(self, geo_pars=None, ffttype='fftw', **kwargs): """ @@ -630,7 +683,7 @@ def __init__(self, geo_pars=None, ffttype='fftw', **kwargs): self.ikernel = None # Get default parameters and update - self.p = u.Param(DEFAULT) + self.p = u.Param(Geo.DEFAULT) self.dtype = kwargs['dtype'] if 'dtype' in kwargs else np.complex128 self.update(geo_pars, **kwargs) self.FFTch = FFTchooser(ffttype) diff --git a/ptypy/core/geometry_bragg.py b/ptypy/core/geometry_bragg.py new file mode 100644 index 000000000..55f8c8982 --- /dev/null +++ b/ptypy/core/geometry_bragg.py @@ -0,0 +1,549 @@ +""" +Geometry management and propagation for Bragg geometry. +""" + +from .. import utils as u +from ..utils.verbose import logger +from geometry import Geo as _Geo +from ..utils.descriptor import EvalDescriptor +from classes import Container, Storage, View +import numpy as np +from scipy.ndimage.interpolation import map_coordinates + +__all__ = ['DEFAULT', 'Geo_Bragg'] + + +local_tree = EvalDescriptor('') +@local_tree.parse_doc() +class Geo_Bragg(_Geo): + """ + Class which presents a Geo analog valid for the 3d Bragg case. + + This class follows the naming convention of: + Berenguer et al., Phys. Rev. B 88 (2013) 144101. + + Indexing into all q-space arrays and storages follows (q3, q1, q2), + which corresponds to (r3, r1, r2) in the so-called natural real + space coordinate system. These coordinates are transformed to + (x, z, y) as described below. + + Defaults: + + [psize] + type = tuple + default = (.065, 172e-6, 172e-6) + help = Rocking curve step (in degrees) and pixel sizes (in meters) + doc = First element is the rocking curve step. + + [propagation] + doc = Only "farfield" is valid for Bragg + + [shape] + type = tuple + default = (31, 128, 128) + help = Number of rocking curve positions and detector pixels + doc = First element is the number of rocking curve positions. + + [theta_bragg] + type = float + default = 6.89 + help = Diffraction angle (theta, not two theta) in degrees + + [resolution] + type = tuple + default = None + help = 3D sample pixel size (in meters) + doc = Refers to the conjugate (natural) coordinate system as (r3, r1, r2). + """ + + def _initialize(self, p): + """ + Parse input parameters, fill missing parameters and set up a + propagator. + """ + self.interact = False + + # Set distance + if self.p.distance is None or self.p.distance == 0: + raise ValueError( + 'Distance (geometry.distance) must not be None or 0') + + # Set frame shape + if self.p.shape is None or (np.array(self.p.shape) == 0).any(): + raise ValueError( + 'Frame size (geometry.shape) must not be None or 0') + else: + self.p.shape = u.expect3(p.shape) + + # Set energy and wavelength + if p.energy is None: + if p.lam is None: + raise ValueError( + 'Wavelength (geometry.lam) and energy (geometry.energy)\n' + 'must not both be None') + else: + self.lam = p.lam # also sets energy through a property + else: + if p.lam is not None: + logger.debug('Energy and wavelength both specified. ' + 'Energy takes precedence over wavelength') + + self.energy = p.energy # also sets lam through a property + + # Pixel size + self.p.psize_is_fix = p.psize is not None + self.p.resolution_is_fix = p.resolution is not None + + if not self.p.psize_is_fix and not self.p.resolution_is_fix: + raise ValueError( + 'Pixel size in sample plane (geometry.resolution) and ' + 'detector plane \n(geometry.psize) must not both be None') + + # Fill pixel sizes + if self.p.resolution_is_fix: + self.p.resolution = u.expect3(p.resolution) + else: + self.p.resolution = u.expect3(1.0) + + if self.p.psize_is_fix: + self.p.psize = u.expect3(p.psize) + else: + self.p.psize = u.expect3(1.0) + + # Update other values + self.update(False) + + # Attach propagator + self._propagator = self._get_propagator() + self.interact = True + + def update(self, update_propagator=True): + """ + Update things which need a little computation + """ + + # Update the internal pixel sizes: 4 cases + if not self.p.resolution_is_fix and not self.p.psize_is_fix: + raise ValueError( + 'Neither pixel size nor sample resolution specified.') + elif not self.p.resolution_is_fix and self.p.psize_is_fix: + dq1 = self.psize[1] * 2 * np.pi / self.distance / self.lam + dq2 = self.psize[2] * 2 * np.pi / self.distance / self.lam + dq3 = np.deg2rad(self.psize[0]) * 4 * np.pi / self.lam * self.sintheta + self.p.resolution[1] = 2 * np.pi / \ + (self.shape[1] * dq1 * self.costheta) + self.p.resolution[2] = 2 * np.pi / (self.shape[2] * dq2) + self.p.resolution[0] = 2 * np.pi / \ + (self.shape[0] * dq3 * self.costheta) + elif self.p.resolution_is_fix and not self.p.psize_is_fix: + dq1 = 2 * np.pi / \ + (self.shape[1] * self.resolution[1] * self.costheta) + dq2 = 2 * np.pi / (self.shape[2] * self.resolution[2]) + dq3 = 2 * np.pi / \ + (self.shape[0] * self.resolution[0] * self.costheta) + self.p.psize[1] = dq1 * self.distance * self.lam / (2 * np.pi) + self.p.psize[2] = dq2 * self.distance * self.lam / (2 * np.pi) + self.p.psize[0] = np.rad2deg( + dq3 * self.lam / (4 * np.pi * self.sintheta)) + else: + raise ValueError( + 'Both pixel size and sample resolution specified.') + + # These are useful to have on hand + self.dq1 = dq1 + self.dq2 = dq2 + self.dq3 = dq3 + + # Establish transforms between coordinate systems + # ...from {x z y} to {r3 r1 r2} + self.A_r3r1r2 = [[1 / self.costheta, 0, 0], + [-self.sintheta / self.costheta, 1, 0], + [0, 0, 1]] + # ...from {r3 r1 r2} to {x z y} + self.A_xzy = [[self.costheta, 0, 0], + [self.sintheta, 1, 0], + [0, 0, 1]] + # ...from {qx qz qy} to {q3 q1 q2} + self.A_q3q1q2 = [[1, self.sintheta / self.costheta, 0], + [0, 1 / self.costheta, 0], + [0, 0, 1]] + # ...from {q3 q1 q2} to {qx qz qy} + self.A_qxqzqy = [[1, -self.sintheta, 0], + [0, self.costheta, 0], + [0, 0, 1]] + + # Update the propagator too + if update_propagator: + self.propagator.update() + + @property + def theta_bragg(self): + return self.p.theta_bragg + + @theta_bragg.setter + def theta_bragg(self, v): + self.p.theta_bragg = v + if self.interact: + self.update() + + @property + def sintheta(self): + return np.sin(np.deg2rad(self.theta_bragg)) + + @property + def costheta(self): + return np.cos(np.deg2rad(self.theta_bragg)) + + @property + def tantheta(self): + return np.tan(np.deg2rad(self.theta_bragg)) + + @property + def resolution(self): + return self.p.resolution + + @resolution.setter + def resolution(self, v): + self.p.resolution[:] = u.expect3(v) + if self.interact: + self.update() + + @property + def psize(self): + return self.p.psize + + @psize.setter + def psize(self, v): + self.p.psize[:] = u.expect3(v) + if self.interact: + self.update() + + @property + def shape(self): + return self.p.shape + + @shape.setter + def shape(self, v): + self.p.shape[:] = u.expect3(v).astype(int) + if self.interact: + self.update() + + def _get_propagator(self): + """ + The real space pixel size in the cartesian system. + """ + prop = BasicBragg3dPropagator(self) + return prop + + def _r3r1r2(self, p): + """ + Transforms a single point from [x z y] to [r3 r1 r2] + """ + return np.dot(self.A_r3r1r2, p) + + def _xzy(self, p): + """ + Transforms a single point from [r3 r1 r2] to [x z y] + """ + return np.dot(self.A_xzy, p) + + def _q3q1q2(self, p): + """ + Transforms a single point from [qx qz qy] to [q3 q1 q2] + """ + return np.dot(self.A_q3q1q2, p) + + def _qzqyqx(self, p): + """ + Transforms a single point from [q3 q1 q2] to [qx qz qy] + """ + return np.dot(self.A_qxqzqy, p) + + def transformed_grid(self, grids, input_space='real', input_system='natural'): + """ + + Transforms a coordinate grid between the cartesian and natural + coordinate systems in real or reciprocal space. + + Parameters + ---------- + grids : 3-tuple of 3-dimensional arrays: (x, z, y), + (r3, r1, r2), (qx, qz, qy), or (q3, q1, q2), + or a 3-dimensional Storage instance. + + input_space: `real` or `reciprocal` + + input_system: `cartesian` or `natural` + + """ + + if isinstance(grids, Storage): + grids = grids.grids() + + # choose transformation operator: 4 cases + if input_space == 'real' and input_system == 'natural': + r3, r1, r2 = grids + z = r1 + self.sintheta * r3 + y = r2 + x = self.costheta * r3 + return x, z, y + elif input_space == 'real' and input_system == 'cartesian': + x, z, y = grids + r1 = z - self.sintheta / self.costheta * x + r2 = y + r3 = 1 / self.costheta * x + return r3, r1, r2 + elif input_space == 'reciprocal' and input_system == 'natural': + q3, q1, q2 = grids + qz = self.costheta * q1 + qy = q2 + qx = q3 - self.sintheta * q1 + return qx, qz, qy + elif input_space == 'reciprocal' and input_system == 'cartesian': + qx, qz, qy = grids + q1 = 1 / self.costheta * qz + q2 = qy + q3 = qx + self.sintheta / self.costheta * qz + return q3, q1, q2 + else: + raise ValueError('invalid options') + + + def coordinate_shift(self, input_storage, input_space='real', + input_system='natural', keep_dims=True, + layer=0): + """ + Transforms a 3D storage between the cartesian and natural + coordinate systems in real or reciprocal space by simply rolling + the axes. It tries to do this symmetrically so that the center + is maintained. + + Note that this transform can be done in any way, and always + involves the choice of a new grid. This method (arbitrarily) + chooses the grid which results from skewing the along the + z/qx direction for real/reciprocal space. + + The shifting is identical to doing a nearest neighbor + interpolation, and it would not be difficult to use other + interpolation orders by instead shifting an index array and + using scipy.ndimage.interpolation.map_coordinates(). But then + you have to decide how to interpolate complex numbers. + + Parameters + ---------- + input_storage : The storage to operate on + + input_space: `real` or `reciprocal` + + input_system: `cartesian` or `natural` + + keep_dims : If True, maintain pixel size and number of pixels. + If False, keeps all the data of the input storage, which means + that the shape of the output storage will be larger than the + input. + + """ + + C_ = Container(data_type=input_storage.dtype, data_dims=3) + S = input_storage + + # Four cases. In real and reciprocal space, these skewing + # operations are done along different axes. For each space, the + # direction of the transform is taken care of. + + if input_space == 'real': + # create a padded copy of the data array + shape = S.shape[1:] + pad = int(np.ceil(self.sintheta * + shape[0] * S.psize[0] / S.psize[1])) + d = np.pad(S.data[layer], pad_width=( + (0, 0), (0, pad), (0, 0)), mode='constant') + # walk along the r3/x axis and roll the r1/z axis. the + # array is padded at the bottom (high indices) so the + # actual shifts have to be positive. + for i in range(shape[0]): + if input_system == 'cartesian': + # roll the z axis in the negative direction for more + # positive x + shift = int( + round((shape[0] - i) * S.psize[0] * self.sintheta / S.psize[1])) + elif input_system == 'natural': + # roll the r1 axis in the positive direction for more + # positive r3 + shift = int( + round(i * S.psize[0] * self.sintheta / S.psize[1])) + d_old = np.copy(d) + d[i, :, :] = np.roll(d[i, :, :], shift, axis=0) + + # optionally crop the new array + if keep_dims: + d = d[:, pad / 2:shape[1] + pad / 2, :] + # construct a new Storage + if input_system == 'cartesian': + new_psize = S.psize * np.array([1 / self.costheta, 1, 1]) + elif input_system == 'natural': + new_psize = S.psize * np.array([self.costheta, 1, 1]) + old_center = S.origin + S.psize * shape / 2 + S_out = C_.new_storage(ID='S0', psize=new_psize, + padonly=False, shape=None) + V = View(container=C_, storageID='S0', coord=old_center, + shape=d.shape, psize=new_psize) + S_out.reformat() + # should use the view here, but there is a bug (#74) + S_out.data[0] = d + + elif input_space == 'reciprocal': + # create a padded copy of the data array + shape = S.shape[1:] + pad = int(np.ceil(self.sintheta * shape[1])) + d = np.pad(S.data[layer], pad_width=( + (0, pad), (0, 0), (0, 0)), mode='constant') + # walk along the q1/qz axis and roll the q3/qx axis. the + # array is padded at the right (high indices) so the + # actual shifts have to be positive. + for i in range(shape[1]): + if input_system == 'cartesian': + # roll the qx axis in the positive direction for more + # positive qz + shift = int(round(i * self.sintheta)) + elif input_system == 'natural': + # roll the q3 axis in the positive direction for more + # negative q1 + shift = int(round((shape[1] - i) * self.sintheta)) + d[:, i, :] = np.roll(d[:, i, :], shift, axis=0) + # optionally crop the new array + if keep_dims: + d = d[pad / 2:shape[0] + pad / 2, :, :] + # construct a new Storage + if input_system == 'cartesian': + new_psize = S.psize * np.array([1, 1 / self.costheta, 1]) + elif input_system == 'natural': + new_psize = S.psize * np.array([1, self.costheta, 1]) + old_center = S.origin + S.psize * shape / 2 + S_out = C_.new_storage(ID='S0', psize=new_psize, + padonly=False, shape=None) + V = View(container=C_, storageID='S0', coord=old_center, + shape=d.shape, psize=new_psize) + S_out.reformat() + # should use the view here, but there is a bug (#74) + S_out.data[0] = d + + return S_out + + def prepare_3d_probe(self, S_2d, auto_center=False, system='cartesian', layer=0): + """ + + Prepare a three-dimensional probe from a two-dimensional incident wavefront. + + Parameters + ---------- + S_2d : Two-dimensional storage holding a (typically complex) + wavefront. The absolute positions are not important, only the + pixel sizes. The first index is interpreted as the quasi- + vertical axis zi (coincident with z at theta=0), and the second + index as yi (coincident always with y and r2). + + center : If true, the input wavefront is centered at the + intensity center of mass. + + system : `cartesian` or `natural`, the coordinate system of the + returned object. + + layer : which layer of the 2d probe to use + + """ + + if auto_center: + raise NotImplementedError + + # storage in natural space to fill with the probe + C = Container(data_type=S_2d.dtype, data_dims=3) + if system == 'natural': + View(C, storageID='S0000', psize=self.resolution, shape=self.shape) + S_3d = C.storages['S0000'] + elif system == 'cartesian': + View(C, storageID='S0000', psize=self.resolution * np.array([self.costheta, 1, 1]), shape=self.shape) + S_3d = C.storages['S0000'] + + # center both storages (meaning that the central pixel is the + # physical origin) + S_3d.center = np.array(S_3d.shape[1:]) / 2 + S_2d.center = np.array(S_2d.shape[1:]) / 2 + + # find the physical coordinates (zi, yi) of each point in the 3d probe + if system == 'natural': + r3, r1, r2 = S_3d.grids() + r3, r1, r2 = r3[0], r1[0], r2[0] # layer 0 + x, z, y = self.transformed_grid((r3, r1, r2), input_space='real', input_system='natural') + zi = x * self.sintheta + z * (1/self.costheta - self.sintheta * self.tantheta) + yi = y + elif system == 'cartesian': + x, z, y = S_3d.grids() + x, z, y = x[0], z[0], y[0] + zi = x * self.sintheta + z * (1/self.costheta - self.sintheta * self.tantheta) + yi = y + + # find the corresponding indices into S.data[layer] + zi[:] = zi / S_2d.psize[1] + S_2d.center[1] + yi[:] = yi / S_2d.psize[1] + S_2d.center[1] + + # interpolate + if np.iscomplexobj(S_2d.data): + S_3d.data[0][:] = map_coordinates(np.abs(S_2d.data[layer]), (zi, yi)) + S_3d.data[0][:] *= np.exp(1j * map_coordinates(np.angle(S_2d.data[layer]), (zi, yi))) + else: + S_3d.data[0][:] = map_coordinates(S_2d.data[layer], (zi, yi)) + + #import ipdb; ipdb.set_trace() + + return S_3d + + def probe_extent_vs_fov(self): + """ + Calculates the extent of the field of view as seen from the + incoming beam. This is the size of the smallest probe (along its + vertical direction zi and horizontal direction yi) which + completely covers the field of view. + + Returns: zi_extent, yi_extent + """ + g = self + b, a, c = g.shape * g.resolution + ap = a + b * g.sintheta + bp = b * g.costheta + y = np.sqrt(ap**2 + bp**2) + gamma = np.arcsin(ap / y) + phi = (np.pi / 2 - gamma - np.deg2rad(g.theta_bragg)) + zi_extent = np.cos(phi) * y + yi_extent = c + + return zi_extent, yi_extent + + +class BasicBragg3dPropagator(object): + """ + Just a wrapper for the n-dimensional FFT, no other Bragg-specific + magic applied here (at the moment). + """ + + def __init__(self, geo=None, ffttype='numpy'): + self.geo = geo + if ffttype == 'numpy': + self.fft = np.fft.fftn + self.ifft = np.fft.ifftn + elif ffttype == 'fftw': + import pyfftw + self.fft = pyfftw.interfaces.numpy_fft.fftn + self.ifft = pyfftw.interfaces.numpy_fft.ifftn + + def update(self): + """ + Update any internal buffers. + """ + return + + def fw(self, a): + return np.fft.fftshift(self.fft(a)) + + def bw(self, a): + return self.ifft(np.fft.ifftshift(a)) diff --git a/ptypy/core/illumination.py b/ptypy/core/illumination.py index 3445c4700..8c2b563fb 100644 --- a/ptypy/core/illumination.py +++ b/ptypy/core/illumination.py @@ -16,133 +16,205 @@ from ..core import geometry from ..utils.verbose import logger from .. import resources +from ..utils.descriptor import EvalDescriptor TEMPLATES = dict() -DEFAULT_aperture = u.Param( - # Aperture form (str). One of: - # - None: no aperture, this may be useful for nearfield - # - 'rect': rectangular aperture - # - 'circ': circular aperture - form='circ', - # Static Noise in the transparent part of the aperture (float). - # Can act like a diffuser but has no wavelength dependency - # Can be either: - # - None : no noise - # - 2-tuple: noise in phase (amplitude (rms), minimum feature size) - # - 4-tuple: noise in phase & modulus (rms, mfs, rms_mod, mfs_mod) - diffuser=None, - # Aperture width or diameter (float). - # May also be a tuple (vertical, horizontal) for size - # in case of an asymmetric aperture - size=None, - # Edge width of aperture in pixel to suppress aliasing (int). - edge=2, - # Size of central stop as a fraction of aperture.size (float). - # If not None: places a central beam stop in aperture. - # The value given here is the fraction of the stop compared to size - central_stop=None, - # Offset between center of aperture and optical axes (float). - # May also be a tuple (vertical, horizontal) for size - # in case of an asymmetric offset - offset=0., - # Rotate aperture by this value (float). - rotate=0., -) -""" Default illumination.aperture parameters. - See :py:data:`.scan.illumination.aperture` and a short listing below """ - -DEFAULT_propagation = u.Param( - # Parameters for propagation after aperture plane - # Propagation to focus takes precedence to parallel propagation - # if focused is not None - # Parallel propagation distance (float). - # If None or 0 : No parallel propagation - parallel=None, - # Propagation distance from aperture to focus (float). - # If None or 0 : No focus propagation - focussed=None, - # Focal spot diameter (float). - spot_size=None, - # Antialiasing factor [not implemented] (float). - # Antialiasing factor used when generating the probe. - # (numbers larger than 2 or 3 are memory hungry) - antialiasing=None, -) - -DEFAULT_diversity = u.Param( - # Noise added on top add the end of initialisation (float). - # Can be either: - # - None : no noise - # - 2-tuple: noise in phase (amplitude (rms), minimum feature size) - # - 4-tuple: noise in phase & modulus (rms, mfs, rms_mod, mfs_mod) - noise=None, - shift=None, - power=1.0, -) - -DEFAULT = u.Param( - override=None, - # User-defined probe (if type is None) (str). - # `None`, path to a *.ptyr file or any python evaluable statement - # yielding a 3d numpy array, If `None` illumination is modeled. - model=None, - # Number of photons in the incident illumination (int, float, None). - photons=None, - recon=u.Param( - rfile='*.ptyr', - ID=None, - layer=None, - ), - stxm=u.Param( - # Label of the scan of whose diffraction data to initialize stxm. - # If None, use own scan_label - label=None, - ), - # Diversity parameters, can be None = no diversity - diversity=DEFAULT_diversity, - # Aperture parameters, can be None = no aperture - aperture=DEFAULT_aperture, - # Propagation parameters, can be None = no propagation - propagation=DEFAULT_propagation, -) -""" Default illumination parameters. See :py:data:`.scan.illumination` - and a short listing below """ - -__all__ = ['DEFAULT', 'init_storage', 'aperture'] - - -def rectangle(grids, dims=None, ew=2): - if dims is None: - dims = (grids.shape[-2] / 2., grids.shape[-1] / 2.) - v, h = dims - V, H = grids - return (u.smooth_step(-np.abs(V) + v/2, ew) - * u.smooth_step(-np.abs(H) + h/2, ew)) - - -def ellipsis(grids, dims=None, ew=2): - if dims is None: - dims = (grids.shape[-2] / 2., grids.shape[-1] / 2.) - v, h = dims - V, H = grids - return u.smooth_step( - 0.5 - np.sqrt(V**2/v**2 + H**2/h**2), ew/np.sqrt(v * h)) +# Local, module-level defaults. These can be appended to the defaults of +# other classes. +illumination_desc = EvalDescriptor('illumination') +illumination_desc.from_string(r""" + [aperture] + type = Param + default = + help = Beam aperture parameters + + [aperture.rotate] + type = float + default = 0. + help = Rotate aperture by this value + doc = + + [aperture.central_stop] + help = size of central stop as a fraction of aperture.size + default = None + doc = If not None: places a central beam stop in aperture. The value given here is the fraction of the beam stop compared to `size` + lowlim = 0. + uplim = 1. + userlevel = 1 + type = float + + [aperture.diffuser] + help = Noise in the transparen part of the aperture + default = None + doc = Can be either: + - ``None`` : no noise + - ``2-tuple`` : noise in phase (amplitude (rms), minimum feature size) + - ``4-tuple`` : noise in phase & modulus (rms, mfs, rms_mod, mfs_mod) + userlevel = 2 + type = tuple + + [aperture.edge] + help = Edge width of aperture (in pixels!) + type = float + default = 2.0 + userlevel = 2 + + [aperture.form] + default = circ + type = None, str + help = One of None, 'rect' or 'circ' + doc = One of: + - ``None`` : no aperture, this may be useful for nearfield + - ``'rect'`` : rectangular aperture + - ``'circ'`` : circular aperture + choices = None,'rect','circ' + userlevel = 2 + + [aperture.offset] + default = 0. + type = float, tuple + help = Offset between center of aperture and optical axes + doc = May also be a tuple (vertical,horizontal) for size in case of an asymmetric offset + userlevel = 2 + + [aperture.size] + default = None + type = float, tuple + help = Aperture width or diameter + doc = May also be a tuple *(vertical,horizontal)* in case of an asymmetric aperture + lowlim = 0. + userlevel = 0 + + [diversity] + default = None + type = Param, None + help = Probe mode(s) diversity parameters + doc = Can be ``None`` i.e. no diversity + userlevel = 1 + + [diversity.noise] + default = (0.5,1.0) + type = tuple + help = Noise in each non-primary mode of the illumination. + doc = Can be either: + - ``None`` : no noise + - ``2-tuple`` : noise in phase (amplitude (rms), minimum feature size) + - ``4-tuple`` : noise in phase & modulus (rms, mfs, rms_mod, mfs_mod) + userlevel = 1 + + [diversity.power] + default = 0.1 + type = tuple, float + help = Power of modes relative to main mode (zero-layer) + uplim = 1.0 + lowlim = 0.0 + userlevel = 1 + + [diversity.shift] + default = None + type = float + help = Lateral shift of modes relative to main mode + doc = **[not implemented]** + userlevel = 2 + + [model] + default = None + type = str, ndarray + help = Type of illumination model + doc = One of: + - ``None`` : model initialitziation defaults to flat array filled with the specified number of photons + - ``'recon'`` : load model from previous reconstruction, see `recon` Parameters + - ``'stxm'`` : Estimate model from autocorrelation of mean diffraction data + - ** : one of ptypys internal image resource strings + - *