set contiguous array before conversion from array to tensor

README.rst

@@ -1,25 +1,6 @@
 SigPy
 =====
 
-.. image:: https://img.shields.io/badge/License-BSD%203--Clause-blue.svg
-	:target: https://opensource.org/licenses/BSD-3-Clause
-
-.. image:: https://travis-ci.com/mikgroup/sigpy.svg?branch=master
-	:target: https://travis-ci.com/mikgroup/sigpy
-
-.. image:: https://readthedocs.org/projects/sigpy/badge/?version=latest
-	:target: https://sigpy.readthedocs.io/en/latest/?badge=latest
-	:alt: Documentation Status
-
-.. image:: https://codecov.io/gh/mikgroup/sigpy/branch/master/graph/badge.svg
-	:target: https://codecov.io/gh/mikgroup/sigpy
-
-.. image:: https://zenodo.org/badge/139635485.svg
-   :target: https://zenodo.org/badge/latestdoi/139635485
-
-
-`Source Code <https://github.com/mikgroup/sigpy>`_ | `Documentation <https://sigpy.readthedocs.io>`_ | `MRI Recon Tutorial <https://github.com/mikgroup/sigpy-mri-tutorial>`_ | `MRI Pulse Design Tutorial <https://github.com/jonbmartin/open-source-pulse-design>`_
-
 SigPy is a package for signal processing, with emphasis on iterative methods. It is built to operate directly on NumPy arrays on CPU and CuPy arrays on GPU. SigPy also provides several domain-specific submodules: ``sigpy.plot`` for multi-dimensional array plotting, ``sigpy.mri`` for MRI reconstruction, and ``sigpy.mri.rf`` for MRI pulse design.
 
 Installation
@@ -48,15 +29,15 @@ SigPy can also be installed through ``pip``::
 	# (optional for plot support) pip install matplotlib
 	# (optional for CUDA support) pip install cupy
         # (optional for MPI support) pip install mpi4py
-	
+
 Installation for Developers
 ***************************
 
 If you want to contribute to the SigPy source code, we recommend you install it with ``pip`` in editable mode::
 
 	cd /path/to/sigpy
 	pip install -e .
-	
+
 To run tests and contribute, we recommend installing the following packages::
 
 	pip install coverage ruff sphinx sphinx_rtd_theme black isort
@@ -101,4 +82,3 @@ Want to do machine learning without giving up signal processing? SigPy has conve
 
 	  x_torch = sigpy.to_pytorch(x)
 	  A_torch = sigpy.to_pytorch_function(A)
-

conda.recipe/install.sh

@@ -0,0 +1,27 @@
+# to install sigpy from scratch:
+
+git clone https://github.com/ZhengguoTan/sigpy.git sigpy_github
+
+cd sigpy_github
+
+conda create -n sigpy_github python=3.10
+conda activate sigpy_github
+
+conda install -c anaconda pip
+
+python -m pip install torch torchvision torchaudio
+
+python -m pip install tqdm
+python -m pip install numba
+python -m pip install scipy
+python -m pip install pywavelets
+python -m pip install h5py
+python -m pip install matplotlib
+
+# please -
+# (1) check your cuda version here
+# (2) log into gpu machine when you are in hpc
+# https://docs.cupy.dev/en/stable/install.html
+pip install cupy-cuda11x
+
+pip install -e .

sigpy/__init__.py

@@ -8,23 +8,13 @@
 that can be used in conjuction with Alg to form an App.
 
 """
-from sigpy import (
-    alg,
-    app,
-    backend,
-    block,
-    config,
-    conv,
-    fourier,
-    interp,
-    linop,
-    prox,
-    pytorch,
-    sim,
-    thresh,
-    util,
-    wavelet,
-)
+
+from .version import __version__ # noqa
+from sigpy import alg, app, config, linop, prox, nlop
+
+from sigpy import (backend, block, conv, interp,
+                   fourier, pytorch, sim, thresh,
+                   util, wavelet)
 from sigpy.backend import *  # noqa
 from sigpy.block import *  # noqa
 from sigpy.conv import *  # noqa
@@ -36,9 +26,7 @@
 from sigpy.util import *  # noqa
 from sigpy.wavelet import *  # noqa
 
-from .version import __version__  # noqa
-
-__all__ = ["alg", "app", "config", "linop", "prox"]
+__all__ = ['alg', 'app', 'config', 'linop', 'prox', 'nlop']
 __all__.extend(backend.__all__)
 __all__.extend(block.__all__)
 __all__.extend(conv.__all__)

sigpy/alg.py

@@ -228,12 +228,13 @@ class ConjugateGradient(Alg):
 
     """
 
-    def __init__(self, A, b, x, P=None, max_iter=100, tol=0):
+    def __init__(self, A, b, x, P=None, max_iter=100, tol=0, verbose=False):
         self.A = A
         self.b = b
         self.P = P
         self.x = x
         self.tol = tol
+        self.verbose = verbose
         self.device = backend.get_device(x)
         with self.device:
             xp = self.device.xp
@@ -278,7 +279,12 @@ def _update(self):
                 util.xpay(self.p, beta, z)
                 self.rzold = rznew
 
-            self.resid = self.rzold.item() ** 0.5
+            self.resid = self.rzold.item()**0.5
+
+            if self.verbose:
+                print("  cg iter: " + "%2d" % (self.iter)
+                      + "; resid: " + "%13.6f" % (self.resid)
+                      + "; norm: " + "%13.6f" % (xp.linalg.norm(self.x.flatten())))
 
     def _done(self):
         return (
@@ -840,6 +846,106 @@ def _done(self):
         return self.iter >= self.max_iter or self.residual <= self.tol
 
 
+class IRGNM(Alg):
+    r"""Iteratively Regularized Gauss-Newton Method (IRGNM)
+
+    Args:
+        A (nlop): Non-linear forward model.
+        y (array): Observation.
+        x (array): Variable.
+        x0 (array): bias for L2 regularization.
+        max_iter (int): maximum number of iterations.
+        alpha (int): regularization strength.
+        redu (float): reduction factor along iterations.
+        inner_iter (int): maximum number of inner iterations.
+        inner_tol (float): tolerance for the inner iterations.
+
+    References:
+        Bauer F., Kannengiesser S. (2007).
+        An alternative approach to the image reconstruction
+        for parallel data acquisition in MRI.
+        Math. Meth. Appl. Sci. 30, 1437-1451.
+
+        Uecker M., Hohage T., Block K. T., Frahm J. (2008)
+        Image reconstruction by regularized nonlinear inversion
+        -- joint estimation of coil sensitivities and image content.
+        Magn. Reson. Med. 60, 674-682.
+
+        Tan Z., Roeloffs V., Voit D., Joseph A. A., Untenberger M.,
+        Merboldt K. D., Frahm J. (2016).
+        Model-based reconstruction for real-time phase-contrast flow MRI:
+        Improved spatiotemporal accuracy.
+        Magn. Reson. Med. 77, 1082-1093.
+
+    """
+    def __init__(self, A, y, x, x0=None,
+                 max_iter=6, alpha=1, redu=2,
+                 inner_iter=100, inner_tol=0.01, verbose=True):
+        self.A = A
+        self.y = y
+        self.device = backend.get_device(y)
+
+        # outer iteration
+        self.max_iter = max_iter
+        self.alpha = alpha
+        self.redu = redu
+
+        # inner iteration
+        self.inner_iter = inner_iter
+        self.inner_tol = inner_tol
+
+        self.verbose = verbose
+
+        xp = self.device.xp
+        with self.device:
+            self.x = sp.to_device(x, device=self.device)
+
+            self.x0 = xp.zeros(x.shape, dtype=y.dtype)
+            if x0 is not None:
+                self.x0 = sp.to_device(x0, device=self.device)
+
+        super().__init__(max_iter)
+
+    def _update(self):
+        xp = self.device.xp
+
+        with self.device:
+            dx = xp.zeros_like(self.x)
+
+            r = self.y - self.A(self.x)
+
+            resid = xp.linalg.norm(r).item()
+
+            if self.verbose:
+                print("gn iter: " + "%2d" % (self.iter)
+                      + "; alpha: " + "%.6f" % (self.alpha)
+                      + "; resid: " + "%.6f" % (resid))
+
+            p = self.A.adjoint(self.x, r)
+            p += self.alpha * (self.x0 - self.x)
+
+            # update dx
+            def AHA(x):
+                return self.A.adjoint(self.x, self.A.derivative(self.x, x)) \
+                    + self.alpha * x
+
+            inner_tol = self.inner_tol * xp.linalg.norm(p).item()
+
+            inner_alg = ConjugateGradient(AHA, p, dx,
+                                          max_iter=self.inner_iter,
+                                          tol=inner_tol)
+            while not inner_alg.done():
+                inner_alg.update()
+
+            # update x
+            self.x += 1. * dx
+            self.alpha /= self.redu
+
+    def _done(self):
+        over_iter = self.iter >= self.max_iter
+        return over_iter
+
+
 class GerchbergSaxton(Alg):
     """Gerchberg-Saxton method, also called the variable exchange method.
     Iterative method for recovery of a signal from the amplitude of linear