How to use the sigpy.linop function in sigpy

To help you get started, we’ve selected a few sigpy examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github mikgroup / sigpy / tests / test_pytorch.py View on Github external
def test_to_pytorch_function_complex(self):
            A = linop.FFT([3])
            x = np.array([1 + 1j, 2 + 2j, 3 + 3j], np.complex)
            y = np.ones([3], np.complex)

            with self.subTest('forward'):
                f = pytorch.to_pytorch_function(
                    A,
                    input_iscomplex=True,
                    output_iscomplex=True).apply
                x_torch = pytorch.to_pytorch(x)
                npt.assert_allclose(f(x_torch).detach().numpy().ravel(),
                                    A(x).view(np.float))

            with self.subTest('adjoint'):
                y_torch = pytorch.to_pytorch(y)
                loss = (f(x_torch) - y_torch).pow(2).sum() / 2
                loss.backward()
github mikgroup / sigpy / tests / test_linop.py View on Github external
self.check_linop_pickleable(A)

                        data_shape = [2, 3, 4]
                        filt = util.randn([4, 2, 2, 3], dtype=dtype)
                        A = linop.ConvolveData(
                            data_shape, filt,
                            mode=mode,
                            multi_channel=True)
                        self.check_linop_linear(A, dtype=dtype, device=device)
                        self.check_linop_adjoint(A, dtype=dtype, device=device)
                        self.check_linop_pickleable(A)

                        data_shape = [2, 3, 4]
                        filt = util.randn([4, 2, 2, 3], dtype=dtype)
                        strides = [2, 2]
                        A = linop.ConvolveData(
                            data_shape, filt,
                            mode=mode, strides=strides,
                            multi_channel=True)
                        self.check_linop_linear(A, dtype=dtype, device=device)
                        self.check_linop_adjoint(A, dtype=dtype, device=device)
                        self.check_linop_pickleable(A)
github mikgroup / sigpy / sigpy / mri / rf / ptx.py View on Github external
xp = device.xp
    with device:
        pulses = xp.zeros((Nc, Nt), xp.complex)

        # set up the system matrix
        if explicit:
            A = rf.linop.PtxSpatialExplicit(sens, coord, dt,
                                            target.shape, b0)
        else:
            A = sp.mri.linop.Sense(sens, coord, weights=None, tseg=tseg,
                                   ishape=target.shape).H

        # handle the Ns * Ns error weighting ROI matrix
        W = sp.linop.Multiply(A.oshape, xp.ones(target.shape))
        if roi is not None:
            W = sp.linop.Multiply(A.oshape, roi)

        # apply ROI
        A = W * A

        # Unconstrained, use conjugate gradient
        if st is None:
            I = sp.linop.Identity((Nc, coord.shape[0]))
            b = A.H * W * target

            alg_method = sp.alg.ConjugateGradient(A.H * A + alpha * I,
                                                  b, pulses, P=None,
                                                  max_iter=max_iter, tol=tol)

        # Constrained case, use SDMM
        else:
            # vectorize target for SDMM
github mikgroup / sigpy / sigpy / mri / rf / linop.py View on Github external
coord[:, 1])))

        # add sensitivities to system matrix
        AFullExplicit = xp.empty(AExplicit.shape)
        for ii in range(nc):
            if three_d:
                tmp = xp.squeeze(sens[ii, :, :, :]).flatten()
            else:
                tmp = sens[ii, :, :].flatten()
            D = xp.transpose(xp.tile(tmp, [coord.shape[0], 1]))
            AFullExplicit = xp.concatenate((AFullExplicit, D * AExplicit),
                                           axis=1)

        # remove 1st empty AExplicit entries
        AFullExplicit = AFullExplicit[:, coord.shape[0]:]
        A = sp.linop.MatMul((coord.shape[0] * nc, 1), AFullExplicit)

        # Finally, adjustment of input/output dimensions to be consistent with
        # the existing Sense linop operator. [nc x nt] in, [dim x dim] out
        Ro = sp.linop.Reshape(ishape=A.oshape, oshape=sens.shape[1:])
        Ri = sp.linop.Reshape(ishape=(nc, coord.shape[0]),
                              oshape=(coord.shape[0] * nc, 1))
        A = Ro * A * Ri

        A.repr_str = 'pTx spatial explicit'

        # output a sigpy linop or a numpy array
        if ret_array:
            return A.linops[1].mat
        else:
            return A
github mikgroup / sigpy / sigpy / app.py View on Github external
return gradh_x

            if self.R is None:
                gamma_primal = self.lamda + self.mu
            else:
                gamma_primal = self.mu

        else:
            gradh = None
            gamma_primal = 0

        if self.G is None:
            proxfc = prox.L2Reg(y.shape, 1, y=y)
            gamma_dual = 1
        else:
            A = linop.Vstack([A, self.G])
            proxf1c = prox.L2Reg(self.y.shape, 1, y=y)
            proxf2c = prox.Conj(self.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:
                self.sigma = 1

            S = linop.Multiply(A.oshape, self.sigma)
            AHA = A.H * S * A
            max_eig = MaxEig(
                AHA,
                dtype=self.x.dtype,
                device=self.x_device,
github mikgroup / sigpy / sigpy / alg.py View on Github external
def Amult(self, x, A, mu, lam):
        M = sp.linop.Multiply(x.shape, np.ones(x.shape) * np.sqrt(1 / mu))
        L = sp.linop.Multiply(x.shape, np.ones(x.shape) * np.sqrt(lam))
        Y = sp.linop.Vstack((A, M, L))
        Ry = sp.linop.Reshape((Y.oshape[0], 1), Y.oshape)
        Y = Ry * Y
        return Y
github mikgroup / sigpy / sigpy / learn / app.py View on Github external
def __init__(self, y, L, lamda=0.005,
                 mode='full', multi_channel=False, device=sp.cpu_device, **kwargs):
        self.y = sp.to_device(y, device)
        self.L = sp.to_device(L, device)
        self.lamda = lamda
        self.mode = mode
        self.multi_channel = multi_channel
        self.device = device

        self._get_params()
        self.A_R = sp.linop.ConvolveInput(self.R_shape, self.L,
                                          mode=self.mode,
                                          input_multi_channel=True,
                                          output_multi_channel=self.multi_channel)
        
        proxg_R = sp.prox.L1Reg(self.R_shape, lamda)
        super().__init__(self.A_R, self.y, proxg=proxg_R, **kwargs)
github mikgroup / sigpy / sigpy / mri / dcf.py View on Github external
"""
    device = sp.Device(device)
    xp = device.xp

    with device:
        w = xp.ones(coord.shape[:-1], dtype=coord.dtype)
        img_shape = sp.estimate_shape(coord)

        # Get kernel
        x = xp.arange(n, dtype=coord.dtype) / n
        kernel = xp.i0(beta * (1 - x**2)**0.5).astype(coord.dtype)
        kernel /= kernel.max()

        G = sp.linop.Gridding(img_shape, coord, width, kernel)
        with tqdm(total=max_iter, desc="PipeMenonDCF",
                  disable=not show_pbar) as pbar:
            for it in range(max_iter):
                GHGw = G.H * G * w
                w /= xp.abs(GHGw)
                resid = xp.abs(GHGw - 1).max().item()

                pbar.set_postfix(resid='{0:.2E}'.format(resid))
                pbar.update()

    return w