Release 0.5 : merging in GPU efforts

.clang-format

@@ -0,0 +1,65 @@
+---
+Language:        Cpp
+# BasedOnStyle:  Google
+AccessModifierOffset: -2
+AlignAfterOpenBracket: true
+AlignEscapedNewlinesLeft: false
+AlignOperands:   true
+AlignTrailingComments: true
+AllowAllParametersOfDeclarationOnNextLine: true
+AllowShortBlocksOnASingleLine: false
+AllowShortCaseLabelsOnASingleLine: false
+AllowShortIfStatementsOnASingleLine: false
+AllowShortLoopsOnASingleLine: true
+AllowShortFunctionsOnASingleLine: All
+AlwaysBreakAfterDefinitionReturnType: false
+AlwaysBreakTemplateDeclarations: true
+AlwaysBreakBeforeMultilineStrings: true
+BreakBeforeBinaryOperators: None
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializersBeforeComma: false
+BinPackParameters: false
+BinPackArguments: false
+ColumnLimit:     80
+ConstructorInitializerAllOnOneLineOrOnePerLine: true
+ConstructorInitializerIndentWidth: 4
+DerivePointerAlignment: true
+ExperimentalAutoDetectBinPacking: false
+IndentCaseLabels: true
+IndentWrappedFunctionNames: false
+IndentFunctionDeclarationAfterType: false
+MaxEmptyLinesToKeep: 1
+KeepEmptyLinesAtTheStartOfBlocks: false
+NamespaceIndentation: None
+ObjCBlockIndentWidth: 2
+ObjCSpaceAfterProperty: false
+ObjCSpaceBeforeProtocolList: false
+PenaltyBreakBeforeFirstCallParameter: 1
+PenaltyBreakComment: 300
+PenaltyBreakString: 1000
+PenaltyBreakFirstLessLess: 120
+PenaltyExcessCharacter: 1000000
+PenaltyReturnTypeOnItsOwnLine: 200
+PointerAlignment: Left
+SpacesBeforeTrailingComments: 2
+Cpp11BracedListStyle: true
+Standard:        Auto
+IndentWidth:     2
+TabWidth:        8
+UseTab:          Never
+BreakBeforeBraces: Allman
+SpacesInParentheses: false
+SpacesInSquareBrackets: false
+SpacesInAngles:  false
+SpaceInEmptyParentheses: false
+SpacesInCStyleCastParentheses: false
+SpaceAfterCStyleCast: false
+SpacesInContainerLiterals: true
+SpaceBeforeAssignmentOperators: true
+ContinuationIndentWidth: 4
+CommentPragmas:  '^ IWYU pragma:'
+ForEachMacros:   [ foreach, Q_FOREACH, BOOST_FOREACH ]
+SpaceBeforeParens: ControlStatements
+DisableFormat:   false
+...
+

.github/workflows/test.yml

@@ -21,37 +21,39 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       max-parallel: 5
+      matrix:
+        python-version: [3.7,3.8,3.9]
 
     steps:
     - uses: actions/checkout@v2
-    - name: Set up Python 3.7
+    - name: Set up Python ${{ matrix.python-version }}
       uses: actions/setup-python@v2
       with:
-        python-version: 3.7
+        python-version:  ${{ matrix.python-version }}
     - name: Add conda to system path
       run: |
         # $CONDA is an environment variable pointing to the root of the miniconda directory
         echo $CONDA/bin >> $GITHUB_PATH
     - name: Install dependencies
       run: |
-        conda env update --file ptypy_core_dependencies.yml --name base
+        conda env update --file dependencies_core.yml --name base
     - name: Prepare ptypy
       run: |
         # Dry install to create ptypy/version.py
         python setup.py install -n
-#    - name: Lint with flake8
-#      run: |
-#        conda install flake8
-#        # stop the build if there are Python syntax errors or undefined names
-#        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
-#        # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
-#        flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+    - name: Lint with flake8
+      run: |
+        conda install flake8
+        # stop the build if there are Python syntax errors or undefined names
+        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
+        # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
+        # flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
     - name: Test with pytest
       run: |
         conda install pytest
         conda install pytest-cov
-        # pytest ptypy/test --doctest-modules --junitxml=junit/test-results.xml --cov=ptypy --cov-report=xml --cov-report=html --cov-config=.coveragerc
-        pytest
+        # pytest ptypy/test -v --doctest-modules --junitxml=junit/test-results.xml --cov=ptypy --cov-report=xml --cov-report=html --cov-config=.coveragerc
+        pytest -v
 #    - name: cobertura-report
 #      if: github.event_name == 'pull_request' && (github.event.action == 'opened' || github.event.action == 'reopened' || github.event.action == 'synchronize')
 #      uses: 5monkeys/cobertura-action@v7

.gitignore

@@ -28,3 +28,4 @@ ptypy/version.py
 /env
 *.egg-info
 .DS_Store
+.ipynb_checkpoints

.travis.yml

@@ -3,6 +3,8 @@ sudo: true
 language: python
 python:
   - 3.7
+compiler:
+  - gcc
 before_install:
   - wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh; # grab miniconda
   - bash miniconda.sh -b -p $HOME/miniconda # install miniconda
@@ -12,6 +14,21 @@ before_install:
   - conda install pyyaml
   - conda info -a # and print the info
 
+  - CUDA=10.1.105-1
+  - CUDA_SHORT=10.1
+  - UBUNTU_VERSION=ubuntu1804
+  - INSTALLER=cuda-repo-${UBUNTU_VERSION}_${CUDA}_amd64.deb
+  - wget http://developer.download.nvidia.com/compute/cuda/repos/${UBUNTU_VERSION}/x86_64/${INSTALLER}
+  - sudo dpkg -i ${INSTALLER}
+  - wget https://developer.download.nvidia.com/compute/cuda/repos/${UBUNTU_VERSION}/x86_64/7fa2af80.pub
+  - sudo apt-key add 7fa2af80.pub
+  - sudo apt update -qq
+  - sudo apt install -y cuda-core-${CUDA_SHORT/./-} cuda-cudart-dev-${CUDA_SHORT/./-} cuda-cufft-dev-${CUDA_SHORT/./-} cuda-curand-dev-${CUDA_SHORT/./-}
+  - sudo apt clean
+  - CUDA_HOME=/usr/local/cuda-${CUDA_SHORT}
+  - LD_LIBRARY_PATH=${CUDA_HOME}/lib64:${LD_LIBRARY_PATH}
+  - PATH=${CUDA_HOME}/bin:${PATH}
+
 env:
   - TEST_ENV_NAME=ptypy_core_dependencies
   - TEST_ENV_NAME=ptypy_full_dependencies
@@ -27,7 +44,7 @@ script:
   - 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
+  - py.test test -v --ignore=ptypy/test/accelerate_tests --cov ptypy --cov-report term-missing # now run the tests
 
 after_script:
   - coveralls

README.rst

@@ -33,15 +33,18 @@ To get started quickly, please find the official documentation at the project pa
 Features
 --------
 
-* **Difference Map** [#dm]_ algorithm engine with power bound constraint
+* **Difference Map** [#dm]_ algorithm engine with power bound constraint [#power]_.
 * **Maximum Likelihood** [#ml]_ engine with preconditioners and regularizers.
+* A few more engines (RAAR, sDR, ePIE, ...).
 
 * **Fully parallelized** (CPU only) using the Massage Passing Interface 
   (`MPI <https://en.wikipedia.org/wiki/Message_Passing_Interface>`_). 
   Simply execute your script with::
 
     $ mpiexec -n [nodes] python <your_ptypy_script>.py
 
+* **GPU acceleration** based on custom kernels, pycuda, and reikna.
+
 * A **client-server** approach for visualization and control based on 
   `ZeroMQ <http://www.zeromq.org>`_ .
   The reconstruction may run on a remote hpc cluster while your desktop
@@ -60,11 +63,11 @@ Installation
 
 Installation should be as simple as ::
 
-   $ sudo python setup.py install
+   $ sudo pip install .
 
 or, as a user, ::
 
-   $ python setup.py install --user
+   $ pip install . --user
 
 
 Dependencies
@@ -130,3 +133,5 @@ References
 .. [#dm] P.Thibault, M.Dierolf *et al.*, *Science* **321**, 7 (2009), `doi <http://dx.doi.org/10.1126/science.1158573>`_
 
 .. [#ml] P.Thibault and M.Guizar-Sicairos, *New J. of Phys.* **14**, 6 (2012), `doi <http://dx.doi.org/10.1088/1367-2630/14/6/063004>`_
+
+.. [#power] K.Giewekemeyer *et al.*, **PNAS 108**, 2 (2007), `suppl. material <https://www.pnas.org/doi/10.1073/pnas.0905846107#supplementary-materials>`__, `doi <https://doi.org/10.1073/pnas.0905846107>`__

archive/array_based/__init__.py

@@ -0,0 +1,7 @@
+'''
+A module for gpu acceleration
+
+'''
+import numpy as np
+COMPLEX_TYPE= np.complex64
+FLOAT_TYPE = np.float32

archive/array_based/array_utils.py

@@ -0,0 +1,87 @@
+'''
+useful utilities from ptypy that should be ported to gpu. These don't ahve external dependencies
+'''
+import numpy as np
+from scipy import ndimage as ndi
+
+
+def dot(A, B, acc_dtype=np.float64):
+    assert A.dtype == B.dtype, "Input arrays must of same data type"
+    if np.iscomplexobj(B):
+        out = np.sum(np.multiply(A, B.conj()).real, dtype=acc_dtype)
+    else:
+        out = np.sum(np.multiply(A, B), dtype=acc_dtype)
+    return out
+
+
+def norm2(A):
+    return dot(A, A)
+
+
+def abs2(input):
+    '''
+
+    :param input. An array that we want to take the absolute value of and square. Can be inplace. Can be complex or real. 
+    :return: The real valued abs**2 array
+    '''
+    return np.multiply(input, input.conj()).real
+
+def sum_to_buffer(in1, outshape, in1_addr, out1_addr, dtype):
+    '''
+    :param in1. An array . Can be inplace. Can be complex or real.
+    :param outshape. An array.
+    :param in1_addr. An array . Can be inplace. Can be complex or real.
+    :param out1_addr. An array . Can be inplace. Can be complex or real.
+    :return: The real valued abs**2 array
+    '''
+    out1 = np.zeros(outshape, dtype=dtype)
+    inshape = in1.shape
+    for i1, o1 in zip(in1_addr, out1_addr):
+        out1[o1[0], o1[1]:(o1[1] + inshape[1]), o1[2]:(o1[2] + inshape[2])] += in1[i1[0]]
+    return out1
+
+def norm2(input):
+    '''
+    Input here could be a variety of 1D, 2D, 3D complex or real. all will be single precision at least.
+    return should be real
+    '''
+    return np.sum(abs2(input))
+
+def complex_gaussian_filter(input, mfs):
+    '''
+    takes 2D and 3D arrays. Complex input, complex output. mfs has len 0<x<=2
+    '''
+    if len(mfs)>2:
+        raise NotImplementedError("Only batches of 2D arrays allowed!")
+
+    if input.ndim == 3:
+        mfs = np.insert(mfs, 0, 0)
+
+    return (ndi.gaussian_filter(np.real(input), mfs) +1j *ndi.gaussian_filter(np.imag(input), mfs)).astype(input.dtype)
+
+def mass_center(A):
+    '''
+    Input will always be real, and 2d or 3d, single precision here
+    '''
+    return np.array(ndi.measurements.center_of_mass(A), dtype=A.dtype)
+
+def interpolated_shift(c, shift, do_linear=False):
+    '''
+    complex bicubic interpolated shift.
+    complex output. This shift should be applied to 2D arrays. shift should have len=c.ndims 
+
+    '''
+    if not do_linear:
+        return ndi.interpolation.shift(np.real(c), shift, order=3, prefilter=True) + 1j*ndi.interpolation.shift(np.imag(c), shift, order=3, prefilter=True)
+    else:
+        return ndi.interpolation.shift(np.real(c), shift, order=1, mode='constant', cval=0, prefilter=False) + 1j * ndi.interpolation.shift(np.imag(c), shift, order=1, mode='constant', cval=0, prefilter=False)
+
+
+def clip_complex_magnitudes_to_range(complex_input, clip_min, clip_max):
+    '''
+    This takes a single precision 2D complex input, clips the absolute magnitudes to be within a range, but leaves the phase untouched.
+    '''
+    ampl = np.abs(complex_input)
+    phase = np.exp(1j * np.angle(complex_input))
+    ampl = np.clip(ampl, clip_min, clip_max)
+    complex_input[:] = ampl * phase

archive/array_based/base.py

@@ -0,0 +1,22 @@
+from collections import OrderedDict
+
+class Adict(object):
+
+    def __init__(self):
+        pass
+
+
+class BaseKernel(object):
+
+    def __init__(self, queue_thread=None, verbose=False):
+
+        self.queue = queue_thread
+        self.verbose = False
+        self.npy = Adict()
+        self.ocl = Adict()
+        self.benchmark = OrderedDict()
+
+
+    def log(self, x):
+        if self.verbose:
+            print(x)