Merge branch 'master' of https://github.com/mikgroup/sigpy

conda.recipe/meta.yaml

@@ -1,5 +1,5 @@
 {% set name = "sigpy" %}
-{% set version = "0.1.15" %}
+{% set version = "0.1.16" %}
 
 package:
   name: '{{ name|lower }}'

docs/conf.py

@@ -26,7 +26,7 @@
 # The short X.Y version
 version = ''
 # The full version, including alpha/beta/rc tags
-release = '0.1.15'
+release = '0.1.16'
 
 
 # -- General configuration ---------------------------------------------------
@@ -83,7 +83,7 @@
 #
 html_theme = 'sphinx_rtd_theme'
 
-html_logo = 'figures/sigpy_logo_white.pdf'
+html_logo = 'figures/sigpy_logo_white.png'
 
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the

docs/guide_basic.rst

@@ -35,7 +35,7 @@ For example to create an array on GPU 1, we can do:
 Note that this can also be done with ``cupy.cuda.Device``, and you can choose to use it as well.
 The main difference is that :class:`sigpy.Device` maps -1 to CPU, and makes it easier to develop CPU/GPU generic code.
 
-.. image:: figures/device.pdf
+.. image:: figures/device.png
    :align: center
 
 To transfer an array between device, we can use :class:`sigpy.to_device`. For example, to transfer a numpy array to GPU 1, we can do:

docs/guide_iter.rst

@@ -3,7 +3,7 @@ Building iterative methods
 
 SigPy provides four abstraction classes (Linop, Prox, Alg, and App) for optimization based iterative methods. Such abstraction is inspired by similar structure in BART.
 
-.. image:: figures/architecture.pdf
+.. image:: figures/architecture.png
    :align: center
 
 The Linop class abstracts a linear operator, and supports adjoint, addition, composing, and stacking. Prepackaged Linops include FFT, NUFFT, and wavelet, and common array manipulation functions. In particular, given a Linop ``A``, the following operations can be performed:

docs/guide_multi_devices.rst

@@ -12,12 +12,12 @@ Another benefit is that an MPI parallelized code can run on both shared memory a
 For example, if we consider the following shared memory configuration (one multi-core CPU and two GPUs),
 and want to run the blue and red tasks concurrently:
 
-.. image:: figures/multiprocess_desired.pdf
+.. image:: figures/multiprocess_desired.png
    :align: center
 
 Then, using MPI, we can split the tasks to two MPI nodes as follows:
 
-.. image:: figures/multiprocess_mpi.pdf
+.. image:: figures/multiprocess_mpi.png
    :align: center
 
 Note that tasks on each MPI node can run on any CPU/GPU device, and in our example, the blue task uses CPU and GPU 0, and

docs/index.rst

@@ -25,6 +25,7 @@
    mri
    mri_linop
    mri_app
+   mri_rf
 
 .. toctree::
    :hidden:

docs/make.bat

@@ -1,36 +1,36 @@
-@ECHO OFF
-
-pushd %~dp0
-
-REM Command file for Sphinx documentation
-
-if "%SPHINXBUILD%" == "" (
-	set SPHINXBUILD=sphinx-build
-)
-set SOURCEDIR=.
-set BUILDDIR=_build
-set SPHINXPROJ=sigpy
-
-if "%1" == "" goto help
-
-%SPHINXBUILD% >NUL 2>NUL
-if errorlevel 9009 (
-	echo.
-	echo.The 'sphinx-build' command was not found. Make sure you have Sphinx
-	echo.installed, then set the SPHINXBUILD environment variable to point
-	echo.to the full path of the 'sphinx-build' executable. Alternatively you
-	echo.may add the Sphinx directory to PATH.
-	echo.
-	echo.If you don't have Sphinx installed, grab it from
-	echo.http://sphinx-doc.org/
-	exit /b 1
-)
-
-%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS%
-goto end
-
-:help
-%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS%
-
-:end
-popd
+@ECHO OFF
+
+pushd %~dp0
+
+REM Command file for Sphinx documentation
+
+if "%SPHINXBUILD%" == "" (
+	set SPHINXBUILD=sphinx-build
+)
+set SOURCEDIR=.
+set BUILDDIR=_build
+set SPHINXPROJ=sigpy
+
+if "%1" == "" goto help
+
+%SPHINXBUILD% >NUL 2>NUL
+if errorlevel 9009 (
+	echo.
+	echo.The 'sphinx-build' command was not found. Make sure you have Sphinx
+	echo.installed, then set the SPHINXBUILD environment variable to point
+	echo.to the full path of the 'sphinx-build' executable. Alternatively you
+	echo.may add the Sphinx directory to PATH.
+	echo.
+	echo.If you don't have Sphinx installed, grab it from
+	echo.http://sphinx-doc.org/
+	exit /b 1
+)
+
+%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS%
+goto end
+
+:help
+%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS%
+
+:end
+popd

docs/mri.rst

@@ -46,6 +46,7 @@ Sampling Functions
 
    sigpy.mri.poisson
    sigpy.mri.radial
+   sigpy.mri.spiral
 
 Simulation Functions
 --------------------

docs/mri_rf.rst

@@ -0,0 +1,36 @@
+MRI RF Design (`sigpy.mri.rf`)
+==============================
+
+.. automodule::
+   sigpy.mri.rf
+
+SLR Pulse Design
+--------------------------
+.. autosummary::
+    :toctree: generated
+    :nosignatures:
+
+    sigpy.mri.rf.slr.dzrf
+    sigpy.mri.rf.slr.root_flip
+    sigpy.mri.rf.slr.dz_gslider_rf
+    sigpy.mri.rf.slr.dz_gslider_b
+    sigpy.mri.rf.slr.dz_hadamard_b
+
+RF Pulse Simulation
+--------------------------
+.. autosummary::
+    :toctree: generated
+    :nosignatures:
+
+    sigpy.mri.rf.sim.abrm
+    sigpy.mri.rf.sim.abrm_nd
+    sigpy.mri.rf.sim.abrm_hp
+
+RF Utility
+--------------------------
+.. autosummary::
+    :toctree: generated
+    :nosignatures:
+
+    sigpy.mri.rf.util.dinf
+

setup.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.1.15
+current_version = 0.1.16
 commit = True
 tag = True
 

setup.py

@@ -10,7 +10,7 @@
     long_description = f.read()
 
 setup(name='sigpy',
-      version='0.1.15',
+      version='0.1.16',
       description='Python package for signal reconstruction.',
       long_description=long_description,
       long_description_content_type="text/x-rst",

sigpy/alg.py

@@ -220,6 +220,7 @@ class ConjugateGradient(Alg):
 
     def __init__(self, A, b, x, P=None, max_iter=100, tol=0):
         self.A = A
+        self.b = b
         self.P = P
         self.x = x
         self.tol = tol
@@ -348,14 +349,13 @@ def __init__(self, proxfc, proxg, A, AH, x, u,
         super().__init__(max_iter)
 
     def _update(self):
-        x_old = self.x.copy()
-
         # Update dual.
         util.axpy(self.u, self.sigma, self.A(self.x_ext))
         backend.copyto(self.u, self.proxfc(self.sigma, self.u))
 
         # Update primal.
         with self.x_device:
+            x_old = self.x.copy()
             util.axpy(self.x, -self.tau, self.AH(self.u))
             backend.copyto(self.x, self.proxg(self.tau, self.x))
 
@@ -553,7 +553,6 @@ def _update(self):
             gradf_x = self.gradf(self.x)
             p = -self.inv_hessf(self.x)(gradf_x)
             self.lamda2 = -xp.real(xp.vdot(p, gradf_x)).item()
-
             if self.lamda2 < 0:
                 raise ValueError(
                     'Direction is not descending. Got lamda2={}. '

sigpy/app.py

@@ -324,16 +324,17 @@ def _get_PrimalDualHybridGradient(self):
             else:
                 proxg = self.proxg
 
-        if self.G is None:
-            proxfc = prox.L2Reg(self.y.shape, 1, y=-self.y)
-            gamma_dual = 1
-        else:
-            A = linop.Vstack([A, self.G])
-            proxf1c = prox.L2Reg(self.y.shape, 1, y=-self.y)
-            proxf2c = prox.Conj(proxg)
-            proxfc = prox.Stack([proxf1c, proxf2c])
-            proxg = prox.NoOp(self.x.shape)
-            gamma_dual = 0
+        with self.y_device:
+            if self.G is None:
+                proxfc = prox.L2Reg(self.y.shape, 1, y=-self.y)
+                gamma_dual = 1
+            else:
+                A = linop.Vstack([A, self.G])
+                proxf1c = prox.L2Reg(self.y.shape, 1, y=-self.y)
+                proxf2c = prox.Conj(proxg)
+                proxfc = prox.Stack([proxf1c, proxf2c])
+                proxg = prox.NoOp(self.x.shape)
+                gamma_dual = 0
 
         if self.tau is None:
             if self.sigma is None:

sigpy/linop.py

@@ -555,20 +555,23 @@ class Diag(Linop):
     Args:
         linops (list of Linops): list of linops with the same input and
             output shape.
-        axis (int or None): If None, inputs/outputs are vectorized
+        iaxis (int or None): If None, inputs are vectorized
+            and concatenated.
+        oaxis (int or None): If None, outputs are vectorized
             and concatenated.
 
     """
 
-    def __init__(self, linops, axis=None):
+    def __init__(self, linops, oaxis=None, iaxis=None):
         self.nops = len(linops)
 
         self.linops = linops
-        self.axis = axis
+        self.oaxis = oaxis
+        self.iaxis = iaxis
         ishape, self.iindices = _hstack_params(
-            [linop.ishape for linop in self.linops], axis)
+            [linop.ishape for linop in self.linops], iaxis)
         oshape, self.oindices = _vstack_params(
-            [linop.oshape for linop in self.linops], axis)
+            [linop.oshape for linop in self.linops], oaxis)
 
         super().__init__(oshape, ishape)
 
@@ -592,26 +595,32 @@ def _apply(self, input):
                     iend = self.iindices[n]
                     oend = self.oindices[n]
 
-                if self.axis is None:
-                    output[ostart:oend] = linop(
+                if self.iaxis is None:
+                    output_n = linop(
                         input[istart:iend].reshape(linop.ishape)).ravel()
                 else:
-                    ndim = len(linop.oshape)
-                    axis = self.axis % ndim
-                    oslc = tuple([slice(None)] * axis + [slice(ostart, oend)] +
-                                 [slice(None)] * (ndim - axis - 1))
-
                     ndim = len(linop.ishape)
-                    axis = self.axis % ndim
+                    axis = self.iaxis % ndim
                     islc = tuple([slice(None)] * axis + [slice(istart, iend)] +
                                  [slice(None)] * (ndim - axis - 1))
 
-                    output[oslc] = linop(input[islc])
+                    output_n = linop(input[islc])
+
+                if self.oaxis is None:
+                    output[ostart:oend] = output_n
+                else:
+                    ndim = len(linop.oshape)
+                    axis = self.oaxis % ndim
+                    oslc = tuple([slice(None)] * axis + [slice(ostart, oend)] +
+                                 [slice(None)] * (ndim - axis - 1))
+
+                    output[oslc] = output_n
 
         return output
 
     def _adjoint_linop(self):
-        return Diag([op.H for op in self.linops], axis=self.axis)
+        return Diag([op.H for op in self.linops],
+                    oaxis=self.iaxis, iaxis=self.oaxis)
 
 
 class Reshape(Linop):
@@ -1155,7 +1164,8 @@ def __init__(self, ishape, axes=None, wave_name='db4', level=None):
         self.wave_name = wave_name
         self.axes = axes
         self.level = level
-        oshape, _ = wavelet.get_wavelet_shape(ishape, wave_name, axes, level)
+        oshape = wavelet.get_wavelet_shape(ishape, wave_name=wave_name,
+                                           axes=axes, level=level)
 
         super().__init__(oshape, ishape)
 
@@ -1188,14 +1198,14 @@ def __init__(self, oshape, axes=None, wave_name='db4', level=None):
         self.wave_name = wave_name
         self.axes = axes
         self.level = level
-        ishape, self.coeff_slices = wavelet.get_wavelet_shape(
-            oshape, wave_name, axes, level)
+        ishape = wavelet.get_wavelet_shape(oshape, wave_name=wave_name,
+                                           axes=axes, level=level)
         super().__init__(oshape, ishape)
 
     def _apply(self, input):
         return wavelet.iwt(
-            input, self.oshape, self.coeff_slices,
-            wave_name=self.wave_name, axes=self.axes, level=self.level)
+            input, self.oshape, wave_name=self.wave_name,
+            axes=self.axes, level=self.level)
 
     def _adjoint_linop(self):
         return Wavelet(self.oshape, axes=self.axes,

sigpy/mri/rf/__init__.py

@@ -0,0 +1,23 @@
+"""This MRI submodule contains functions and classes for MRI pulse design.
+
+It currently provides tools for SLR pulse design and RF pulse simulation.
+Tools for designing pTx, adiabatic, and other types of rf pulses,
+as well as gradient and trajectory designers will be added soon.
+
+Explore RF design tutorials at `sigpy-rf-tutorials`_.
+
+See in-progress features at `sigpy-rf`_.
+
+.. _sigpy-rf-tutorials: https://github.com/jonbmartin/sigpy-rf-tutorials
+.. _sigpy-rf: https://github.com/jonbmartin/sigpy-rf
+
+"""
+from sigpy.mri.rf import slr, sim, util
+from sigpy.mri.rf.sim import *  # noqa
+from sigpy.mri.rf.slr import *  # noqa
+from sigpy.mri.rf.util import *  # noqa
+
+__all__ = []
+__all__.extend(slr.__all__)
+__all__.extend(sim.__all__)
+__all__.extend(util.__all__)