mikgroup / sigpy Goto Github PK

Describe the bug
An earlier PR attempted to resolve an issue where calling sigpy.mri.samp.poisson changed the numpy state (#54). However, numba only supports a select number of numpy operations and no longer supports the RandomState class in functions decorated with nb.jit. This is likely the source of the error.

To Reproduce

from sigpy.mri.samp import poisson

# Not the issue occurs with any arbitrary image size and acceleration.
mask = poisson((320, 320), accel=6, seed=80)

Expected behavior
Expect to have no error in call while still ensuring that global numpy state does not change when seed is specified in sigpy.mri.samp.poisson.

When running linear regression using the following code :

import numpy as np
from sigpy.learn import app

n = 2
k = 5
m = 4
batch_size = n

X = np.random.randn(n, k)
y = np.random.randn(n, m)

alpha = 1 / np.linalg.svd(X, compute_uv=False)[0]**2
mat = app.LinearRegression(X, y, batch_size, alpha, max_iter=300).run()

It gives me this error :

~/.local/lib/python3.8/site-packages/sigpy/learn/app.py in __init__(self, input, output, batch_size, mu, lamda, max_iter, max_inner_iter, device, checkpoint_path)
    338             self.output_t = xp.empty(
    339                 (batch_size, ) + output.shape[1:], dtype=dtype)
--> 340             self.t_idx = sp.ShuffledNumbers(num_batches)
    341 
    342         self._get_A()

AttributeError: module 'sigpy' has no attribute 'ShuffledNumbers'

Something like the sk-image implementation , but with a matrix representation.

Is your feature request related to a problem? Please describe.
Currently espirit_calib expects calibration region as input, and outputs only one set of maps.

Describe the solution you'd like

• Extract calibration region inside the function
• Output arbitrary number of maps

Hi! I have been trying to reconstruct a radially sampled dataset acquired on a Siemens Biograph mMR scanner. I have a matlab script that is able to decode the proprietary raw data format, and then I can load a kspace in python using h5py.

The full dimension of the dataset is [512,2400,47,6], [#samples, #spokes, #slices, #coils ] respectively. When I load it into python the shape changes to (6, 47, 2400, 512), as expected.

I tried to go straigth to 3D, but either I made something wrong (likely) or your code needs an insane amount of RAM to manage that dataset.

Then I simplified the problem and went to 2D, with a kspace shape of (6, 2400, 512), with the hope of being able to follow the tutorial more closely.
I created a coord structure as:
coord = sp.mri.radial([spokes,nx,2], img_shape=[128,128], golden = True)

and a dcf as suggested in the tutorial:
dcf = (coord[..., 0]**2 + coord[..., 1]**2)**0.5

I am able to reconstruct using a standard:
img_grid = sp.nufft_adjoint(ksp * dcf, coord)

but the results are not even close to the abdominal image I was expecting.

Any idea what could have gone wrong?

Describe the bug
After installing Sigpy on Anaconda (win7), when I try to import sigpy I get the following error:

I dont have GPU, so I did not install Cupy.

AttributeError: module 'cupy' has no attribute 'ElementwiseKernel'

To Reproduce
run the following code on Jupyter:

import sigpy as sp

Expected behavior

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-8-7f0de9ccd830> in <module>
----> 1 import sigpy as sp

~\Anaconda3\envs\Sigpy_MRI\lib\site-packages\sigpy\__init__.py in <module>
      9 
     10 """
---> 11 from sigpy import alg, app, config, linop, prox
     12 
     13 from sigpy import (backend, block, conv, interp,

~\Anaconda3\envs\Sigpy_MRI\lib\site-packages\sigpy\alg.py in <module>
      5 """
      6 import numpy as np
----> 7 from sigpy import backend, util
      8 
      9 

~\Anaconda3\envs\Sigpy_MRI\lib\site-packages\sigpy\util.py in <module>
    475     import cupy as cp
    476 
--> 477     _axpy_cuda = cp.ElementwiseKernel(
    478         'S a, T x',
    479         'T y',

AttributeError: module 'cupy' has no attribute 'ElementwiseKernel'

​

Desktop (please complete the following information):

• OS: win7

I'm trying to reconstruct data acquired by Siemens Skyra scanner with Star VIBE (Golden Angle) sequence. I'm using sigpy.mri.radial for generating the trajectories and then sigpy.linop.NUFFT for reconstructing the data. Please help.

ISODELAY time of a simulated RF pulse is needed to make it really useful.

What's I want
sigpy.mri.rf.slr.dzrf() return designed RF pulse and its isodelay.

Describe the bug
Starting around v1.20, numpy has deprecated generic data types (np.float, np.bool, np.complex) in favor of standard python types (float, bool, complex). This is causing a flurry of deprecation warnings in sigpy, where default types may be any of the generic numpy dtypes.

Recommendation
Set explicit dtypes with explicit precision (e.g. np.float64 instead of np.float)

We installed the following packages

CUDA Version: 10.2
cudatoolkit 10.2.89 h8f6ccaa_8 conda-forge
cudnn 8.2.0.53 h2c0ae14_0 conda-forge
cupy 9.2.0 py38h14045e6_0 conda-forge
sigpy 0.1.23 py38_0 frankong

and ran the following example
x = cupy.array([0.1,0.2,0.3])
sigpy.convolve(x, x)

We get the following error. Can you please let us know what is the issue and how to resolve this ?

RuntimeError Traceback (most recent call last)
in
----> 1 sigpy.convolve(x, x)

~/.conda/envs/py38new/lib/python3.8/site-packages/sigpy/conv.py in convolve(data, filt, mode, strides, multi_channel)
40 if config.cudnn_enabled:
41 if np.issubdtype(data.dtype, np.floating):
---> 42 output = _convolve_cuda(data, filt,
43 mode=mode, strides=strides,
44 multi_channel=multi_channel)

~/.conda/envs/py38new/lib/python3.8/site-packages/sigpy/conv.py in _convolve_cuda(data, filt, mode, strides, multi_channel)
339 output = xp.empty((B, c_o) + p, dtype=data.dtype)
340 filt = util.flip(filt, axes=range(-D, 0))
--> 341 cudnn.convolution_forward(data, filt, None, output,
342 pads, s, dilations, groups,
343 auto_tune=auto_tune,

cupy/cudnn.pyx in cupy.cudnn.convolution_forward()

cupy/cudnn.pyx in cupy.cudnn._find_algorithm_fwd()

RuntimeError: No available algorithm found.

As of 85739ee, MPI.COMPLEX is passed into Bcast. However, this is assuming the input array is complex. Ideally, there should be a dictionary that describes the relationship between data-types in numpy/cupy and the MPI types. I did try something akin to this but I got a dictionary error when trying to index complex types.

Hi, I have been looking at the functions for undersampling described in the docs page, but I cannot see any example with them in the notebooks provided and the description for me is a bit confusing. Could anyone please post a demo on how to use these undersampling functions?

I am trying to replicate compressed sensing experiments using radial and spiral trajectories for my bachelor thesis.
Thank you

Hello,

I have been using SigPy to perform L1-wavelet CS and TV reconstruction of undersampled data stored in h5 files. One of the steps in my workflow involves estimating Sensitivy maps of undersampled data with a certain ACS region using the ESPIRiT app.

My script produces the desired maps and reconstruction results while computing on a local system with no GPU. Because of the longer computation time, I deployed the same script in a remote machine with GPU

The only changes made to the script was the addition of device=sp.Device(0) to the ESPIRiT app and convertion of cupy.ndarray to numpy.ndarray as seen in the lines below:

mps_msk = mr.app.EspiritCalib(Masked_ksp, calib_width=num_low_frequency, device=sp.Device(0)).run()
mps_msk = cupy.asnumpy(mps_msk) 

Expected behavior
Computation of Sensitivity maps (in the local system) for Undersampled data with ACS region = 13 produces the sensitivity maps and its coil combined image (slice number 6) as seen below.

 mps_msk = mr.app.EspiritCalib(Masked_ksp, calib_width=num_low_frequency).run()

Issue
The ESPIRiT computation on the same data (with same ACS region) performed in the remote machine with GPU(device=sp.Device(0) enabled) lead to faster computation but produced unexpected cropping of the maps as seen in the below images.

Problem
This unexpected behavior of cropping in estimated maps leads to cropping or poor reconstruction of the Image as seen below.

Desktop (Local machine):

• Machine: Macbook Pro mid 2012
• OS: macOS Catalina
• Version: 10.15.7
• Graphics: Intel HD Graphics 4000 1536 MB
• Memory: 8 GB 1600 MHz DDR3
• Processor: 2,5 GHz Dual-Core Intel Core i5

Desktop (Remote server):

• OS: Manjaro Linux 64bit / Linux 5.10.2-2-MANJARO
• Graphics: GeForce RTX 2080 Ti 10GB

Additional context
After estimating the sensitivity maps, I convert cupy.ndarray to numpy.ndarray to store the data in the h5 file. As far as my understanding goes, this does not cause any changes in the data.

My input k-space shape to the Espirit app is [coils, RO, PE] and calib_width=num_low_frequency = 13 for 8-fold accelerated data.

I kindly request you to provide some insight and solution to this issue, so that I could make my work reproducible.

PS: I would also like to thank all the developers and contributors of SigPy for providing us with this awesome toolkit.

In https://sigpy.readthedocs.io/en/latest/generated/sigpy.prox.L2Reg.html#sigpy.prox.L2Reg

is the proximal operator a function of y or z. It appears that the base class sigpy.prox.Prox uses y as the input to the prox operator. But y is a parameter for the L2Reg class called "bias".

I think "z" is the bias based on what the description is. My tests confirm this.

See below for a test which passes based on the assumption that z is the bias parameter and y is the input to the prox operator.

Also the statement in the docs above says min instead of argmin.

Hi, when using a non-cartesian trajectory, JSENSE estimates coil sensitivity maps with a shape that I didn't think was intuitive. It seems the shape of the coil sensitivity maps is determined by sp.estimate_shape(coord).

For example, when using a 3D centre-out radial trajectory that corresponds to an image size of 128, the min/max coordinates are something like

(kx_min, kx_max) = (-63.99089, 63.984524)
(ky_min, ky_max) = (-63.988026, 63.93878)
(kz_min, kz_max) = (-64.0, 63.914078)

When I give the corresponding array of k-space coordinates to sp.estimate_shape, the output estimated shape is [127, 127, 127], instead of the expected [128, 128, 128].

I'm not sure what the best fix for this is... The simplest thing to do is to change the int() function in estimate_shape to round():

shape = [round(coord[..., i].max().item() - coord[..., i].min().item())
         for i in range(ndim)]

But this will change historic results, and still might not produce the expected shape for certain trajectories.

Alternatively, all functions and classes that rely on estimate_shape can be modified so that a shape can be provided explicitly as an argument, and only when a shape isn't provided explicitly will estimate_shape be used. Functions and classes that rely on estimate_shape are:

nufft_adjoint (this function can already accept a shape explicitly)
JsenseRecon
ConvSense
ConvImage
pipe_menon_dcf

The change is relatively straight forward... but since so many different functions/classes will be changed, I thought I'd ask here first for opinions and discussion.

Hi SigPy group,

First, thank you for creating this powerful software - I've found it easy to use and very helpful for understanding CS in MRI.

I've been working on a problem for a few weeks now: I want to build an app that is capable of solving a multi-regularizer problem. Specifically, I want to solve the following optimization problem:

where lambdaW and lambdaTV are penalties for the wavelet and total variation components, respectively. This optimization problem can also be found in Miki Lustig's 2007 paper.

I came here to see if you have solved this issue in SigPy before, or if you have any advice? So far, I have tried reformulating the problem to use temporary variables for x in each of the regularizer components. However, I haven't been able to get my head around getting to the final solution (i.e. how to recombine the temporary variables). Perhaps it is my inexperience with Python, but I can't quite figure out where to go next.

Thanks for your help, and I look forward to hearing back from you!
Joey

Copied below is the most recent version of an app I tried to create using the current apps as a template: (please be critical)

class TVWaveletRecon(sp.app.LinearLeastSquares):
    r"""L1 Wavelet and total variation regularized reconstruction.
    
    Wavelet is good at preserving edges and low contrast information while TV 
    is efficient at suppressing noise and streaking artifacts.

    Considers the problem

    .. math::
        \min_x \frac{1}{2} \| A x - y \|_2^2 + \lambdaW \| W x \|_1 + \lambdaTV \| G x \|_1

    where A is the sampling operator,
    W is the wavelet operator,
    x is the image, and y is the k-space measurements.

    Args:
        y (array): k-space measurements.
        mps (array): sensitivity maps.
        lamdaW (float): regularization parameter for the wavelet component.
        lamdaTV (float): regularization parameter for the finite difference component.
        weights (float or array): weights for data consistency.
        coord (None or array): coordinates.
        wave_name (str): wavelet name.
        device (Device): device to perform reconstruction.
        coil_batch_size (int): batch size to process coils.
        Only affects memory usage.
        comm (Communicator): communicator for distributed computing.
        **kwargs: Other optional arguments.

    References:
        Lustig, M., Donoho, D., & Pauly, J. M. (2007).
        Sparse MRI: The application of compressed sensing for rapid MR imaging.
        Magnetic Resonance in Medicine, 58(6), 1082-1195.
        
        Zangen, Z., Khan, W., Babyn, P., Cooper, D., Pratt, I., Carter, Y. (2013)
        Improved Compressed Sensing-Based Algorithm for Sparse-View CT Image Reconstruction.
        Computational and Mathematical Methods in Medicine.
        10.1155/2013/185750

    """

    def __init__(self, y, mps, lamdaW, lamdaTV,
                 weights=None, coord=None,
                 wave_name='db4', device=sp.cpu_device,
                 coil_batch_size=None, comm=None, show_pbar=True,
                 transp_nufft=False, **kwargs):
        weights = _estimate_weights(y, weights, coord)
        if weights is not None:
            y = sp.to_device(y * weights**0.5, device=device)
        else:
            y = sp.to_device(y, device=device)

        A = linop.Sense(mps, coord=coord, weights=weights,
                        comm=comm, coil_batch_size=coil_batch_size,
                        transp_nufft=transp_nufft)
        img_shape = mps.shape[1:]
        
        # Wavelet
        W = sp.linop.Wavelet(img_shape, wave_name=wave_name)
        # Finite difference
        G = sp.linop.FiniteDifference(A.ishape)
        
        proxg1 = sp.prox.UnitaryTransform(sp.prox.L1Reg(W.oshape, lamdaW), W)
        proxg2 = sp.prox.L1Reg(G.oshape, lamdaTV)
        
        def g(input):
            device = sp.get_device(input)
            xp = device.xp
            with device:
                return lamdaW * xp.sum(xp.abs(W(input))).item() + lamdaTV * xp.sum(xp.abs(input)).item()
        if comm is not None:
            show_pbar = show_pbar and comm.rank == 0
            
        # Call super().__init(...) to call the __init(...) of the parent class,
        # sp.app.LinearLeastSquares
        super().__init__(A, y, proxg=proxg1, g=g, show_pbar=show_pbar, **kwargs)


        def h(input):
            device = sp.get_device(input)
            xp = device.xp
            with device:
                return lamdaW * xp.sum(xp.abs(W(input))).item() + lamdaTV * xp.sum(xp.abs(input)).item()
        if comm is not None:
            show_pbar = show_pbar and comm.rank == 0
        
        super().__init__(A, y, proxg=proxg2, g=h, G=G, show_pbar=show_pbar, **kwargs)

Code First: https://github.com/ShannonZ/playground

Why is the reconstructed image so obscure
ksp = np.load('ksp.npy')
img_shape = ksp.shape

F = sp.linop.FFT(ksp.shape, axes=(0,1))
mask = abs(ksp)>0
P = sp.linop.Multiply(ksp.shape, mask)
W = sp.linop.Wavelet(img_shape)
wav = W * F.H * ksp

A = P * F * W.H

lamda = 0.005
proxg = sp.prox.L1Reg(wav.shape, lamda)
alpha = 1
wav_thresh = proxg(alpha, wav)

max_iter = 100
alpha = 1

def gradf(x):
return A.H * (A * x - ksp)

wav_hat = np.zeros(wav.shape, np.complex)
alg = sp.alg.GradientMethod(gradf, wav_hat, alpha, proxg=proxg, max_iter=max_iter)

while not alg.done():
alg.update()
print('\rL1WaveletRecon, Iteration={}'.format(alg.iter), end='')

pl.ImagePlot(W.H(wav_hat))
Any demo about recon from normal single coil sampled data
I can only find tutorials about reconstruction from multi-coil data.

Describe the bug
'sigpy.mri.app.TotalVariationReconandsigpy.mri.app.SenseReconbehave differently if parametercoil_batch_size` is smaller than the total number of coils

Expected behavior
The final reconstructed image should be the same, no matter the coil batch size used. The doc itself states that that parameter affects only performance

Screenshots
RECONSTRUCTED IMAGE WITH coil_batch_size < 20 (with 20 = # coils)

RECONSTRUCTED IMAGE WITH coil_batch_size = 20 (with 20 = # coils)

Additional context
Maybe the coil_batch_size parameter operates on a different dimension than axis=0?
These are the shapes of the arrays I am trying to pass to either apps:

Hello,
I am trying to setup a script to reconstruct a dataset acquired with a radialVIBE sequence. The coordinate sampling scheme is basically a radial sampling in the axial plane x-y, and a full cartesian sampling in the z direction.

I built my own function for this purpose, but the only way I am able to perfome the 3D reconstruction, is if I loop each slice at a time, along the z direction. If I pass the all kspace to the function sp.nufft_adjoint, I get all sorts of errors.

I have been thinking that maybe the problem is in the shape of the kspace.
So far I am trying to pass it as: [n_coils, n_slices, n_spokes, n_ro]

Let me know if you need more information.
Thanks

When nufft is applied to a real image, and then adjoint nufft to bring back to image space, the returned image should have imaginary values close to zero. But when I'm trying to do the same with nufft and adjoint nufft of sigpy, the imaginary part isn't really close to zero.
When I tried to do the same thing, including same angles, with PyNufft, it gave me desired results.
I'm wondering why such thing is happening?

One more thing, when oversampling is not desired, I'm passing oversamp as 1, it returns nan values (happens during the calculation of apodization).

Thanks for your help in advance :)

Hi,

Thanks for making Sigpy available. I would like to use and benchmark Sigpy. I'm not yet familiar with it, I might have missed something.

Installation (and tutorial) run smootly using pip for the libraries and using the command python3 setup.py install. I'm using a laptop on Ubuntu 18.04 with

• Intel® Core™ i7-4910MQ CPU @ 2.90GHz × 8
• Nvidia Quadro 4Go
• 32 Go Ram.

I didn't succeed to use/activate multi-cpu usage for performing operation either IFFT and coil map estimation or the L1WaveletRecon. For instance, EspiritCalib never finished while JsenseRecon was actually running on but one cpu as well as the L1WaveletRecon. Do you have any example to follow ?

The code is here in case: https://github.com/valeryozenne/sigpy-3D-ifft-cms/blob/master/main.py

Thanks in advance,
Valéry

ftype (string): type of filter to use: 'ms' (sinc), 'pm' (Parks-McClellan equal-ripple), 'min' (minphase using factored pm), 'max' (maxphase using factored pm), 'ls' (least squares).

Hi Frank,

I recently wanted to try out sigpy and had a version installed, but didn't know which one. Finding that out was surprisingly difficult, as sigpy does not have a sigpy.version member. Instead, I actually started comparing the source code that I found in my python path.

Of course, with pip (and probably conda), there are other ways of finding out the installed version, but not when installing it manually.

Therefore, I think it'd be a good idea to add a version string to sigpy.

I don't know much about python packaging, but it could probably live in a file and just be read into, for example, setup.py. This would also make sure that that version is always accurate.

Best regards,
Christian

I'm trying to use the sigpy.gridding method but having some issues with this, mainly with regard to specifying the kernel. I'm trying to do spline interpolation but when I specify kernel='spline' I get an error which I understand after looking in the source code where the kernel is specified as an array. Can you provide some guidance for how to use this function by specifying the kernel as an array? In particular, I'd be interested to know if it would be possible to do nearest neighbour gridding.

Thanks!

When running the L1WaveletRecon app on my GPU, python is throwing this error, as though it doesn’t understand what the ’S’ in ’S a’ stands for, on line 482 of util.py. My knowledge of cupy is limited but reading the code, it looks like _xpay_cuda defines a CUDA kernel that calculates x + a*y, where a is a scalar (which is what ’S’ signifies) and x and y are tensors (’T’).

To reproduce:
Invoke L1WaveletRecon on the GPU

Error text:
Traceback (most recent call last):
File "./prom_nufft_test_3d.py", line 69, in
show_pbar=False, device=sp.Device(0)).run()
File "/nas/home/grissowa/sigpy-rf/sigpy/app.py", line 85, in run
self.alg.update()
File "/nas/home/grissowa/sigpy-rf/sigpy/alg.py", line 63, in update
self._update()
File "/nas/home/grissowa/sigpy-rf/sigpy/alg.py", line 187, in _update
util.axpy(self.x, -self.alpha, self.gradf(self.x))
File "/nas/home/grissowa/sigpy-rf/sigpy/util.py", line 441, in axpy
_axpy_cuda(a, x, y)
File "cupy/core/_kernel.pyx", line 540, in cupy.core._kernel.ElementwiseKernel.call
File "cupy/core/_kernel.pyx", line 587, in cupy.core._kernel.ElementwiseKernel._decide_params_type
File "cupy/core/_kernel.pyx", line 299, in cupy.core._kernel._decide_params_type_core
KeyError: ’S'

My nvcc —version output:
nvcc: NVIDIA (R) Cuda compiler driver
Copyright (c) 2005-2019 NVIDIA Corporation
Built on Sun_Jul_28_19:07:16_PDT_2019
Cuda compilation tools, release 10.1, V10.1.243

I am running this in a Conda environment with python 3.6, with CUDA 10.1.

Hi again,
sorry for the "improper" use of Github Issues, but I don't see other ways of getting in touch with you.
I began experimenting with the use of the sp.app.LinearLeastSquares class to build iterative reconstruction algorithms.

What I'd like to do is to perform a regularized reconstruction, while using a 'target' image (let's say a fully sampled anatomy) as prior information.
At the moment, after building a forward operator A which seems to work for my problem, I am doing the following:

super().__init__(A, y, max_iter = max_iter, show_pbar=show_pbar, 
                 z = reference, solver = 'ConjugateGradient', **kwargs)

But the results are not as expected. Am I interpreting the meaning of the z variable in the wrong way, or should I look somewhere else to solve the problem?

Thanks

Is your feature request related to a problem? Please describe.
Currently SigPy functions get the device from the first argument, and perform operations in that device. This is inconsistent with CuPy, which simply uses the device in the current context.

Describe the solution you'd like
Change all SigPy functions to use the device in the current context, with CPU as default. If inputs are inconsistent, an error will occur.

Describe alternatives you've considered
Keep the same interface for compatibility reasons. But this requires users to actively think whether they are using a CuPy function or a SigPy function.

Additional context
None

Describe the bug
Setting the seeds for random number generators in relevant libraries (cupy, numpy, random, etc.) should guarantee that the sensitivity map estimates are identical when using the espirit and jsense calibration apps; but this is not the case.

To Reproduce
Colab: https://colab.research.google.com/drive/113Xkjkuadl476thWy8nlKMjb5VHzBD96?usp=sharing

Expected behavior
Setting the seed for the random number generators for all relevant libraries should guarantee that the maps are the same.

Reproduce
This code is reproducible on colab using the most recent SigPy version (0.1.22) and cupy 7.4.0. Note colab does not support higher versions of cupy just yet, but when this was run locally using sigpy v0.1.22 and cupy 9.0.0, the same issues were seen.

Describe the bug
Output masks produced by sigpy.mri.poisson are not within the expected +/-0.1 bounds of the expected acceleration (passed in as an argument). We also see this in cases where crop_corner=False and in the trivial case acceleration=1.0.

To Reproduce

# desired acceleration
desired_accel = 1.0  # also tried with 4.0 and 6.0

mask = sigpy.mri.poisson(
            (512, 80), 
            accel=desired_accel,
            calib=(16, 16),
            dtype=np.bool,
            crop_corner=False,
        )

acc = 512 * 80 / mask.sum()   # between 7.2-7.6 depending on the seed

Expected behavior
The output mask should correspond to an acceleration closer to the desired_accel

Desktop (please complete the following information):

• OS: Ubuntu 20.04
• SigPy Version: 0.1.20

Describe the bug
sigpy.mri.poisson hangs when both the seed and calib arguments are passed in.

To Reproduce
See this colab below
https://colab.research.google.com/drive/1k0w6KjkPHx8U8K8dXmaFVi9hO1Cg1Z1e?usp=sharing

Expected behavior
Passing both these arguments together should not cause the program to hang

When using 'Sigpy.mri.poission' in pytorch's 'Dataloader' to generate random mask(iterative generation), it is easy to occur 'Aborted (core dumped)' in any step and shut down my process.
I want to know is there any way to solve this problem.
Thanks!

Describe the bug
I am running into issues with computing the calibration matrix for a multicoil 3D image (ndim = 4). The array_to_blocks function only accepts inputs with ndim = {1, 2, 3}.

To Reproduce
Steps to reproduce the behavior:

1. Load multicoil 3D ksp data; I'm working with dimensions [ncoils, RO, PE1, PE2] = [12, 256, 218, 170]
2. Run mr.app.EspiritCalib(ksp).run()
3. Error raised: "Traceback (most recent call last):
   File "", line 1, in
   File "/cluster/home/brian.nghiem/.conda/envs/sigpy_env/lib/python3.7/site-packages/sigpy/mri/app.py", line 422, in init
   mat = sp.array_to_blocks(calib, kernel_shape, kernel_strides)
   File "/cluster/home/brian.nghiem/.conda/envs/sigpy_env/lib/python3.7/site-packages/sigpy/block.py", line 103, in array_to_blocks
   raise ValueError('Only support ndim=1, 2, or 3, got {}'.format(ndim))""

Computer
OS: Linux (remote cluster)

Describe the bug
Throw “ValueError: axes don't match array” when doing fft with 2D ksp data.

The kspace data was acquired using a single slice SE sequence (256*192).

To Reproduce
I made some changes in the demo to reproduce the same error,

import numpy as np
import sigpy as sp
import sigpy.mri as mr
import sigpy.plot as pl

ksp = np.load('data/cartesian_ksp.npy')
ksp = ksp[0,:,:]
img_shape = ksp.shape

F = sp.linop.FFT(ksp.shape, axes=(0,1))
pl.ImagePlot(F.H * ksp, z=0, title=r'$F^H y$')
Expected behavior
FFT should support 2-N dimensional data.

Desktop (please complete the following information):

Nothing about OS.

I have a script that allows me to decode vendor-specific raw data from the Siemens Biograph mMR into .cfl files (kspace and non-cartesian trajectory) that can be fed into BART for the reconstruction.
I was wondering if sigpy has an interface that can allow me to use the cfl python library (e.g. function cfl.readcfl() and cfl.writecfl()) as input and use the way more pythonic tools provided by sigpy for the reconstruction.

In setting up SigPy on a new machine, we found that cupy is casting inputs to the element wise kernel as const. Switching the input to outputs is what cupy source does. I have a working update: https://github.com/kmjohnson3/sigpy/tree/cupy_const_updates and will issue a pull request when more validated.

Hello all,

I make use of data from tutorial to try to test CS for radial data:
01-gridding-reconstruction

The sensitivity map is generated from re-gridding K-Space data.

ksp = np.load('data/projection_ksp.npy')
coord = np.load('data/projection_coord.npy')
dcf = (coord[..., 0]**2 + coord[..., 1]**2)**0.5

img_grid_tune = sp.nufft_adjoint(ksp * dcf, coord, oversamp=1, width=2)
ksp_grid = fourier.fft(img_grid_tune, axes=(-2,-1) )
mps = mr.app.EspiritCalib(ksp_grid).run()

Unfortunately, the result looks bad

lamda = 0.005
img_tv = mr.app.TotalVariationRecon(ksp * dcf, mps, lamda, coord=coord).run()
pl.ImagePlot(img_tv, title='Total Variation Regularized Reconstruction')

The whole procedure is attached.
01-gridding-reconstruction - Jupyter Notebook.pdf

Thanks for any kind help~

Best Regards,
Kingaza

Is your feature request related to a problem? Please describe.
A normal operator for the NUFFT linear operator in conjunction with Pull Request #86 should help in increasing reconstruction speed for cases where the NUFFT psf fits in memory.

Describe the solution you'd like
An NUFFT normal operator akin to BART.

Additional context
I talked with @frankong and we decided it'll be best to have this discussion on here to have it track-able.

What I've tried so far

• Let F be the NUFFT linear operator with input dimension "n". We are working in 1D for simplicity.
• From what I understand, the normal operator N = FH * F is linear and time invariant, and hence can be effectively modeled by a convolution.
• Let d be a delta and p = FH * F * d be the impulse response/ point spread function/ convolution kernel. This has support equal to "n".
• In order to model the linear convolution with FFTs, we have to zero-pad the input to "2n". Let R denote the (center) zero-padding from "n" to "2n".
• Let P = FFT(R(p)). Then, in principle, the following should hold: F^H * F = R.H * FFT.H * P * FFT * R.
• However, in practice, I get higher than expected errors.

1D Example

This gives me an error of ~13%.

import sigpy as sp
import sigpy.plot as pl
import sigpy.mri as mr

n      = 256
devnum = 1
zpad   = 2

device = sp.Device(devnum)
xp = device.xp

# Linear operators
coords = mr.spiral(n * 1E-3, n, 2, 2, 32, 1.5, 0.27, 1.8)[:, 0, None]
nF = sp.linop.NUFFT((n,), coords) # Non-uniform
uF = sp.linop.FFT([int(k * zpad) for k in nF.ishape]) # Uniform
R  = sp.linop.Resize(uF.ishape, nF.ishape)

with device:

    # Calculating PSF
    d = xp.zeros((nF.ishape), dtype=xp.complex64)
    d[d.shape[0]//2] = 1
    h = nF.H * nF * d
    
    # Creating Toeplitz operator
    psf = xp.zeros([int(k * zpad) for k in d.shape], xp.complex64)
    psf[(psf.shape[0] - d.shape[0])//2 : (psf.shape[0] - d.shape[0])//2 + d.shape[0]] = h
    PSF = uF(psf)
    T = R.H * uF.H * sp.linop.Multiply(uF.oshape, PSF) * uF * R 
    
    # Testing
    nrmse = lambda x, y: 100 * xp.linalg.norm(x/xp.linalg.norm(x) - y/xp.linalg.norm(y))
    vec = xp.random.standard_normal(T.ishape) + 1j * xp.random.standard_normal(T.ishape)
    print("Random test: %0.2f%%" % nrmse(nF.H * nF * vec, T * vec))
    print("Delta test:  %0.2f%%" % nrmse(h, T(d)))

2D Example

This gives me an error of ~70%.

import sigpy as sp
import sigpy.plot as pl
import sigpy.mri as mr

n = 256
devnum = 1
zpad = 2

device = sp.Device(devnum)
xp = device.xp

# Linear operators
coords = mr.spiral(n * 1E-3, n, 2, 2, 32, 1.5, 0.27, 1.8)
nF = sp.linop.NUFFT((n, n), coords) # Non-uniform
uF = sp.linop.FFT([int(k * zpad) for k in nF.ishape], axes=(0, 1)) # Uniform
R  = sp.linop.Resize(uF.ishape, nF.ishape)

with device:

    # Calculating PSF
    d = xp.zeros((nF.ishape), dtype=xp.complex64)
    d[d.shape[0]//2, d.shape[1]//2] = 1
    h = nF.H * nF * d
    
    # Creating Toeplitz operator
    psf = xp.zeros([int(k * zpad) for k in d.shape], xp.complex64)
    psf[(psf.shape[0] - d.shape[0])//2 : (psf.shape[0] - d.shape[0])//2 + d.shape[0],
        (psf.shape[1] - d.shape[1])//2 : (psf.shape[1] - d.shape[1])//2 + d.shape[1]] = h
    PSF = uF(psf)
    T = R.H * uF.H * sp.linop.Multiply(uF.oshape, PSF) * uF * R 
    
    # Testing
    nrmse = lambda x, y: 100 * xp.linalg.norm(x/xp.linalg.norm(x) - y/xp.linalg.norm(y))
    vec = xp.random.standard_normal(T.ishape) + 1j * xp.random.standard_normal(T.ishape)
    print("Random test: %0.2f%%" % nrmse(nF.H * nF * vec, T * vec))
    print("Delta test:  %0.2f%%" % nrmse(h, T(d)))

Hello,

I have Anaconda python 3.8.
For installation, I gave the command
$conda install -c frankong sigpy

But it gives the following error:

UnsatisfiableError: The following specifications were found
to be incompatible with the existing python installation in your environment:

Specifications:

• sigpy -> python[version='>=3.6,<3.7.0a0|>=3.7,<3.8.0a0']

Your python: python=3.8

If python is on the left-most side of the chain, that's the version you've asked for.
When python appears to the right, that indicates that the thing on the left is somehow
not available for the python version you are constrained to. Note that conda will not
change your python version to a different minor version unless you explicitly specify
that.

I can try downgrading my python to 3.7 but is it possible to upgrade sigpy to work with python 3.8.

Thanks,
Hemant

Describe the bug
I try to interpolate the radial data from a cartesian k-space, and the coord data is generated by sigpy.mri.radial()
But the reconstruction seems quite strange, and the sampled k-space value is not as expected

To Reproduce

import numpy as np
import matplotlib.pyplot as plt

import sigpy as sp
from sigpy.mri import samp
import sigpy.plot as pl


def fftc_2d(_data):
    data = _data.copy()
    data = np.fft.ifftshift(data)
    data = np.fft.fft2(data, norm="ortho")
    return np.fft.fftshift(data)    

def ifftc_2d(_data):
    data = _data.copy()
    data = np.fft.ifftshift(data)
    data = np.fft.ifft2(data, norm="ortho")
    return np.fft.fftshift(data)    

def show_image(image):
    plt.figure()
    plt.imshow(np.abs(image), cmap='gray')
    plt.show()

def show_kspace(kspace):
    plt.figure()
    plt.imshow(np.log(np.abs(kspace)), cmap='gray')
    plt.show()

# create shepp-logan image and k-space
image_shape = (160, 160)
shepp_image = sp.shepp_logan(image_shape)
show_image(shepp_image)
shepp_kspace = fftc_2d(shepp_image)
show_kspace(shepp_kspace)

# sample as radial style
spoke_num, sample_num = 60, 160
coord_shepp = samp.radial([spoke_num, sample_num, 2], image_shape)
pl.ScatterPlot(coord_shepp, title='Trajectory')

# if no offset to the coord, no image can be reconstructed
# BUT the k-space is very strange (the data in middle column shall have the largest magitude)
coord_offset = 80
kspace_nu = sp.interpolate(shepp_kspace, coord_shepp + coord_offset)
show_kspace(kspace_nu)

dcf_nu = (coord_shepp[..., 0]**2 + coord_shepp[..., 1]**2)**0.5
image_nu = sp.nufft_adjoint(kspace_nu * dcf_nu, coord_shepp)
show_image(image_nu)

Expected behavior
A clear and concise description of what you expected to happen.

Screenshots
If applicable, add screenshots to help explain your problem.

Desktop (please complete the following information):

• OS: Win10
• sigpy 0.1.23

"bloch_forward" is not callable using sp.mri
Generates module 'sigpy.mri' has no attribute 'bloch_forward'

When I call mr.poisson((640, 320), accel = 4, calib = (20, 20), crop_corner = False), I get the following error:

free(): invalid next size (normal)
Aborted (core dumped)

But if I change calib to (19,19) or (21,21) everything works fine. I know this is similar to the other open issues about mri.poisson, but perhaps it may provide some more context for the issue?

In sigpy/mri/util.py the function 'tseg_off_res_b_ct' will always fail in line 69: t = np.linspace(0, T, T/dt)
T/dt is always a float and not an integer, which will return an error.

Is your feature request related to a problem? Please describe.
There is no wavelet transform in GPU. Currently, SigPy moves array to CPU and uses pywavelet to perform the wavelet transform, which is the main bottleneck for compressed sensing MRI reconstructions.

Describe the solution you'd like
A GPU wavelet transform can leverage the GPU multi-channel convolution operations already wrapped in SigPy. Low pass and high pass filtering can be done in one operation using the output channels. Strides can be used to incorporate subsampling.

Describe alternatives you've considered
Implement the GPU kernels for wavelet transforms. But this would be less optimized than cuDNN convolutions.

Describe the bug
Images in your documentation do not render.

To Reproduce
Go the readthedocs and look at the documentation pages. The logo and other images do not render.

Expected behavior
I should see images instead of placeholders.

Screenshots
If applicable, add screenshots to help explain your problem.

Desktop (please complete the following information):

• OS: Red Hat Linux
• Browser: Chrome
• Version 76.0.3809.100 (Official Build) (64-bit)

Is your feature request related to a problem? Please describe.
CuPy is a useful library for GPU-based ndarray operations; however there can often be issues with installation. SigPy has been useful in making the CuPy dependency optional; however, in the current state, the dependency is not really optional.

In the case where the cupy library is installed (pip/conda install) but not properly configured (e.g. wrong nvcc version is active), the importing cupy fails as shown below. Currently, sigpy attempts to set cupy_enabled by checking util.find_spec("cupy") is not None. However, this is not sufficient, as the cupy package can still be installed but can fail on import.

>>> import cupy
Traceback (most recent call last):
  File "/home/usr/miniconda3/envs/dl_ss_env/lib/python3.7/site-packages/cupy/__init__.py", line 16, in <module>
    from cupy import _core  # NOQA
  File "/home/usr/miniconda3/envs/dl_ss_env/lib/python3.7/site-packages/cupy/_core/__init__.py", line 1, in <module>
    from cupy._core import core  # NOQA
ImportError: libcublas.so.10: cannot open shared object file: No such file or directory

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/usr/miniconda3/envs/dl_ss_env/lib/python3.7/site-packages/cupy/__init__.py", line 37, in <module>
    raise ImportError(_msg) from e
ImportError: CuPy is not correctly installed.

If you are using wheel distribution (cupy-cudaXX), make sure that the version of CuPy you installed matches with the version of CUDA on your host.
Also, confirm that only one CuPy package is installed:
  $ pip freeze

If you are building CuPy from source, please check your environment, uninstall CuPy and reinstall it with:
  $ pip install cupy --no-cache-dir -vvvv

Check the Installation Guide for details:
  https://docs.cupy.dev/en/latest/install.html

original error: libcublas.so.10: cannot open shared object file: No such file or directory

Describe the solution you'd like
Even if import cupy fails when cupy is installed, sigpy should not fail. Currently, sigpy passes along the ImportError raised by cupy, which causes sigpy to error out.

Because cupy is an optional dependency, sigpy should catch the exception, print a warning message to the user, and continue on with cupy_enabled=False. I believe the code below should do the trick:

# in sigpy.config
from importlib import util
import warnings

cupy_enabled = util.find_spec("cupy") is not None
if cupy_enabled:
    try:
        import cupy
    except ImportError as e:
        warnings.warn(e)
        cupy_enabled=False
if cupy_enabled:  # pragma: no cover
    cudnn_enabled = util.find_spec("cupy.cuda.cudnn") is not None
    nccl_enabled = util.find_spec("cupy.cuda.nccl") is not None
else:
    cudnn_enabled = False
    nccl_enabled = False

I don't know if this is of interest, or if you are aware of this behavior, but I thought I would mentioned it here, in case anyone was struggling with this "bug" ...

Describe the bug
When creating an instance of an 'app' class, we pass to the __init__ function all the input variables needed to create an instance of an 'alg'. If using a GPU, this sends a lot of data (i.e. coord, weights, ...) to the GPU's RAM memory.
Running the run() method, we can then perform the optimization and get back the result. However, the memory allocated on the GPU is not freed automatically. Furthermore, even if we delete to instance of the app class, after completing the optimization, part of that memory remain allocated, untill we restart the python kernel. At least, this is what happens using python notebooks.

A few notes:

• the only algorithm no affected by this, which is able to correctly free GPU's memory is Conjugate Gradient. I haven't been able to figure out why.
• To check the memory used by your app you may do as follows:
  -- create a subclass of sp.app.LinearLeastSquares
  -- in the __init__() method add the following lines:

import cupy
self.mempool = cupy.get_default_memory_pool()
self.pinned_mempool = cupy.get_default_pinned_memory_pool()

-- then define a _post_update(self) method as:

def _post_update(self):
        print(self.mempool.used_bytes())              
        print(self.mempool.total_bytes())             
        print(self.pinned_mempool.n_free_blocks())
        return

You will se that after each alg.update(), the memory occupancy keeps growing.
After the run() is finished you can delete the instance of the app, but such memory will remain flagged.

The only way I found so far to overcome this issue is to add a del self.alg to the __output(self) method of the app:

def _output(self):
        x = sp.to_device(self.x, sp.cpu_device)
        del self.x
        del self.alg
        return x

this seems like the only way to me to deallocate all the structures created by the algorithm on the GPU.

Describe the bug
Calling the sigpy.mri.poisson function modifies the default state for numpy.random. This can be an issue if the user has specified a different default random seed earlier in their code.

To Reproduce
Follow call trace for sigpy.mri.poisson -> sigpy.mri._poisson
https://github.com/mikgroup/sigpy/blob/master/sigpy/mri/samp.py#L139

Expected behavior
If seed is specified, sigpy.mri.poisson should use a temporary random state based on that seed. This way the default state will not be modified

I updated sigpy to 0.1.21 (from 0.1.18) and now I can't run the JsenseRecon app on the GPU. I've tried downgrading sigpy again, and installing in a fresh environment and still get the same error, so I'm stuck, even though it was working fine before I upgraded (via conda). Maybe someone can help me figure out what's going on.

Code used to reproduce:

import numpy as np
import sigpy as sp
import sigpy.mri as mri
import sigpy.plot as plt

dataInDir = "/home/ltorres/data/forMatt/"
ksp = np.load("{}ksp.npy".format(dataInDir))
coord = np.load("{}coord.npy".format(dataInDir))
dcf = np.load("{}dcf.npy".format(dataInDir))

device = 0
ksp = sp.to_device(ksp, device)
coord = sp.to_device(coord, device)
dcf = sp.to_device(dcf, device)
mps = mri.app.JsenseRecon(
    ksp,
    coord=coord,
    weights=dcf,
    device=device,
).run()

Resulting Error:

JsenseRecon:   0%|                                                         | 0/10 [00:00<?, ?it/s]Traceback (most recent call last):
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/linop.py", line 95, in apply
    output = self._apply(input)
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/linop.py", line 1657, in _apply
    multi_channel=self.multi_channel)
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/conv.py", line 143, in convolve_filter_adjoint
    multi_channel=multi_channel)
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/conv.py", line 311, in _complex
    outputr = func(data1r, data2r, *kargs, **kwargs)
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/conv.py", line 412, in _convolve_filter_adjoint_cuda
    tensor_core=tensor_core)
  File "cupy/cudnn.pyx", line 1817, in cupy.cudnn.convolution_backward_filter
  File "cupy/cudnn.pyx", line 1469, in cupy.cudnn._find_algorithm_bwd_filter
RuntimeError: No available algorithm found.

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/linop.py", line 95, in apply
    output = self._apply(input)
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/linop.py", line 362, in _apply
    output = linop(output)
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/linop.py", line 122, in __call__
    return self.__mul__(input)
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/linop.py", line 131, in __mul__
    return self.apply(input)
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/linop.py", line 98, in apply
    raise RuntimeError('Exceptions from {}.'.format(self)) from e
RuntimeError: Exceptions from <[8, 1, 16, 16, 16]x[8, 23, 23, 23]> ConvolveFilterAdjoint Linop>.

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/ltorres/projects/extreme_mri/bug.py", line 22, in <module>
    device=device,
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/app.py", line 85, in run
    self.alg.update()
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/alg.py", line 63, in update
    self._update()
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/alg.py", line 414, in _update
    self.min1()
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/mri/app.py", line 336, in min_mps_ker
    show_pbar=False).run()
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/app.py", line 213, in __init__
    self._get_alg()
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/app.py", line 247, in _get_alg
    self._get_ConjugateGradient()
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/app.py", line 264, in _get_ConjugateGradient
    AHy = self.A.H(self.y)
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/linop.py", line 122, in __call__
    return self.__mul__(input)
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/linop.py", line 131, in __mul__
    return self.apply(input)
  File "/home/ltorres/anaconda3/envs/sigpyUpdated/lib/python3.7/site-packages/sigpy/linop.py", line 98, in apply
    raise RuntimeError('Exceptions from {}.'.format(self)) from e
RuntimeError: Exceptions from <[8, 16, 16, 16]x[8, 14827820]> Reshape * ConvolveFilterAdjoint * FFT * NUFFTAdjoint * Reshape * Sum * Multiply Linop>.
JsenseRecon:   0%|                                                         | 0/10 [00:01<?, ?it/s]```

In case it's relevant, data is non-cartesian, downsampled to fit on GPU and organized as:

ksp: [coils, spokes, readouts]
coord:[spokes, readouts, dimension]
dcf: [spokes, readouts]

Package versions:

cudatoolkit - 10.2.89 
cudnn - 7.6.5 
cupy - 8.0.0
numba - 0.51.2
numpy - 1.19.1
sigpy - 0.1.18/0.1.21

Similar to #77 - Newer versions of cupy don't fail on import cupy but rather when cudnn and nccl libraries are fetched.

Because cupy is an optional dependency, and there are use cases where users share virtualenv's on a cluster, but each machine on the cluster has different hardware, it may make sense to make failed package importing a warning.

Solution

Add try/except to protect the import of cudnn and nccl components

# in sigpy.config
from importlib import util
import warnings

cupy_enabled = util.find_spec("cupy") is not None
if cupy_enabled:
    try:
        import cupy
    except ImportError as e:
        warnings.warn(e)
        cupy_enabled=False
if cupy_enabled:  # pragma: no cover
    try:
        cudnn_enabled = util.find_spec("cupy.cuda.cudnn") is not None
    except ImportError as e:
        warnings.warn(e)
        cudnn_enabled = False
    try:
        nccl_enabled = util.find_spec("cupy.cuda.nccl") is not None
    except ImportError as e:
        warnings.warn(e)
        nccl_enabled = False
else:
    cudnn_enabled = False
    nccl_enabled = False

Hey folks,

Thanks for creating this awesome tool!

I was wondering if you'd be interested in including a general method that may facilitate validation of adjoint linear operations.
Something as simple as checking if <d, Ap> == <A^H d, p> holds for a pair of primal/dual random vectors {p, d}.
I thought such feature could ease common developments, e.g., checking if the adjoint of a multi-dim spatial difference operator is correctly implemented.
(My apologize if it already exists and I've missed it.)

I feel like the following patch can do the trick, but am not so sure if it's compatible with things like different backends, etc.

From e88c3fdf0e8dd3db5e38b722316453b0d3ea8d74 Mon Sep 17 00:00:00 2001
From: <[email protected]>
Date: Thu, 19 May 2022 19:56:24 -0400
Subject: [PATCH] A tool for adjoint operator validation.

---
 sigpy/linop.py | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/sigpy/linop.py b/sigpy/linop.py
index 0efd1af..cdfd533 100644
--- a/sigpy/linop.py
+++ b/sigpy/linop.py
@@ -109,6 +109,23 @@ class Linop():
     def _normal_linop(self):
         return self.H * self

+    def _valid_adjoint(self, rtol=1e-05, atol=1e-08, equal_nan=False) -> bool:
+        """Validate if adjoint is correctly implemented
+
+        For random primal/dual vectors, {p, d}, expect <d, Ap> == <Aᴴd, p>.
+        """
+        p1 = np.random.randn(self.ishape)+1j*np.random.randn(self.ishape)
+        d1 = self.apply(p1)
+
+        d2 = np.random.randn(self.oshape)+1j*np.random.randn(self.oshape)
+        p2 = self.H.apply(d2)
+
+        res_fwd = (d2.ravel()).conj().T @ d1.ravel()  # <  d, Ap>
+        res_adj = (p2.ravel()).conj().T @ p1.ravel()  # <Aᴴd,  p>
+
+        return np.isclose(res_fwd, res_adj,
+                          rtol=rtol, atol=atol, equal_nan=equal_nan)
+
     @property
     def H(self):
         r"""Return adjoint linear operator.
--
2.25.1

Thanks!