tsadakane commented 3 years ago

I think I found a bug at checkDevices codes. For example, https://github.com/CERN/TIGRE/blob/3977143a6dd8390366b1f33fa82ae8adc325233b/Common/CUDA/voxel_backprojection.cu#L625-L627

dev is not the ID of GPU but an iterator.
The device ID of the GPU is gpuids[dev]

so this code will behave unexpectedly if the users GPUs are

0: Type A, 1: Type B, 2: Type B

and gpuids = [1,2].

Unfortunately my GPUs are

0: Type A, 1: Type A, 2: Type B

I, therefore, cannot test the exact case, but it seems clear.

Additionally, https://github.com/CERN/TIGRE/blob/3977143a6dd8390366b1f33fa82ae8adc325233b/Common/CUDA/voxel_backprojection.cu#L626 is not necessary.

I think this check can be done by GpuIds class.

The bugfix is done at 9fe4ed7c04e9202d2f79813ff21d7119c46f8973. Removing redundant line and refactoring are done at 8e03d5e5a44e32a3eae903bf560007df6749f172.

Tests are WIP

[x] Ax ( [cone, parallel] x [Siddon, interpolated] ) (parallel does not check devices.)
[x] Atb ([cone, parallel] x [matched, FDK] )
[x] POCS_TV
[x] POCS_TV2
[x] tvdenoising
[x] voxel_backprojection2

AnderBiguri commented 3 years ago

All looks good here, and its indeed a real bug.

Let me know if I can help with anything.

tsadakane commented 3 years ago

Thank you for your proposal. The problem is that I have not yet run the unchecked functions.

tsadakane commented 3 years ago

TEST CODE

#%%Initialize
import numpy as np
import itertools
from tigre.utilities.gpu import GpuIds
import tigre
from tigre.utilities import sample_loader
from tigre.utilities import CTnoise
import tigre.algorithms as algs
from tigre.utilities import gpu
from tigre.utilities.im3Dnorm import im3DNORM

gpuids = gpu.getGpuIds("GeForce RTX 2080 Ti")
# gpuids = gpu.getGpuIds()# All GPUs are selected

if False: # True
    list_proj_type     = ["interpolated", "Siddon"]
    list_backproj_type = ["matched"     , "FDK" ]
    list_geomode       = ["cone"        , "parallel"]

    #%%
    for proj_type, backproj_type, geomode in itertools.product(list_proj_type, list_backproj_type, list_geomode):
        print("{}, {}, {}".format(proj_type, backproj_type, geomode))
        #%% Geometry
        geo = tigre.geometry_default(high_resolution=False)
        geo.mode = geomode

        #%% Load data and generate projections
        # define angles
        angles = np.linspace(0, 2 * np.pi, 100, endpoint=False)
        # Load thorax phatom data
        head = sample_loader.load_head_phantom(geo.nVoxel)
        # generate projections
        projections = tigre.Ax(head, geo, angles, projection_type=proj_type, gpuids=gpuids)
        # add noise
        noise_projections = CTnoise.add(projections, Poisson=1e5, Gaussian=np.array([0, 10]))

        #%% Reconstruct image using OS-SART and FDK

        # BP
        imgs = tigre.Atb(noise_projections, geo, angles, krylov=backproj_type, gpuids=gpuids)

    niter = 3
    imgOSSART = algs.ossart(noise_projections, geo, angles, niter, gpuids=gpuids)

    #%% Show the results
    # tigre.plotimg(np.concatenate([imgs, imgOSSART], axis=1), dim="z")
    tigre.plotimg(np.concatenate([imgs], axis=1), dim="z")

else:
    proj_type     = "interpolated"
    geo = tigre.geometry_default(high_resolution=False)
    geo.mode = "cone"

    #%% Load data and generate projections
    # define angles
    angles = np.linspace(0, 2 * np.pi, 16, endpoint=False)
    # Load thorax phatom data
    head = sample_loader.load_head_phantom(geo.nVoxel)
    # generate projections
    projections = tigre.Ax(head, geo, angles, projection_type=proj_type, gpuids=gpuids)
    # add noise
    noise_projections = CTnoise.add(projections, Poisson=1e5, Gaussian=np.array([0, 10]))    

    ## TEST tvdenoising via tvdenoise via im3ddenoise via fista
    imgFISTA_hightv = algs.fista(
        noise_projections, geo, angles, 3, hyper=2.0e6, tviter=3, tvlambda=20, gpuids=gpuids
    )

    ## TEST pocs_tv via minTV via minimizeTV via ASD_POCS
    epsilon = (
        im3DNORM(tigre.Ax(algs.fdk(noise_projections, geo, angles, gpuids=gpuids), geo, angles, gpuids=gpuids) - noise_projections, 2)
        * 0.15
    )
    alpha = 0.002
    ng = 3
    lmbda = 1
    lambdared = 0.9999  # you generally want 1
    alpha_red = 0.95
    ratio = 0.94
    verb = True

    imgASDPOCS = algs.asd_pocs(
        noise_projections,
        geo,
        angles,
        3,  # these are very important
        tviter=ng,
        maxl2err=epsilon,
        alpha=alpha,  # less important.
        lmbda=lmbda,
        lmbda_red=lambdared,
        rmax=ratio,
        verbose=verb,
        gpuids=gpuids,
    )

    ## TEST aw_pocs_tv via AwminTV via minimizeAwTV via AwASD_POCS

    imgAWASDPOCS = algs.awasd_pocs(
        noise_projections,
        geo,
        angles,
        3,  # these are very important
        tviter=ng,
        maxl2err=epsilon,
        alpha=alpha,  # less important.
        lmbda=lmbda,
        lmbda_red=lambdared,
        rmax=ratio,
        verbose=verb,  # AwASD_POCS params
        delta=np.array([-0.005]),
        gpuids=gpuids,
    )

There are two gpuids at the beginning of the code. If the first one is used, no warning displayed and if the second one is used, warnings are displayed. In both cases, the program ended gracefuly.

tsadakane commented 3 years ago

This test checks that

this modification does not induce sudden death of the program,
if the given gpuids points multiple types of GPUs, the warnings are shown
if the given gpuids points single type of GPU, the warnings are not shown

Unfortunately I could not test the case that I mentioned at the first comment of this issue.

The test is done only on Python, but the CUDA code is common.

Modifiled file	Modified CUDA function	Tested function
POCS_TV.cu	pocs_tv	algs.asd_pocs
POCS_TV2.cu	aw_pocs_tv	algs.awasd_pocs
Siddon_projection.cu	siddon_ray_projection	Ax("cone", "Siddon")
ray_interpolated_projection.cu	interpolation_projection	Ax("cone", "interpolated")
tvdenoising.cu	tvdenoising	algs.fista
voxel_backprojection.cu	voxel_backprojection	Atb("cone", "FDK")
voxel_backprojection2.cu	voxel_backprojection2	Atb("cone", "matched")

AnderBiguri commented 3 years ago

Looks good on my side, unfortunately I can not test your first test either, however I see no failure case for that, so I suggest we merge this and deal with the issue in the future if it comes.

A couple of comments:

1 - Can you add a message to the Changelog? 2 - there are few mexWarnMsgIdAndTxt, do you mind changing the last sentences?

Currently they read as:

mexWarnMsgIdAndTxt("minimizeTV:POCS_TV:GPUselect","Detected one (or more) different GPUs.\n This code is not smart enough to separate the memory GPU wise if they have different computational times or memory limits.\n First GPU parameters used. If the code errors you might need to change the way GPU selection is performed. \n Siddon_projection.cu line 275.");

We should just

mexWarnMsgIdAndTxt("minimizeTV:POCS_TV:GPUselect","Detected one (or more) different GPUs.\n This code is not smart enough to separate the memory GPU wise if they have different computational times or memory limits.\n First GPU parameters used. If the code errors you might need to change the way GPU selection is performed");

as now the GPU selection is done via MATLAB/python. (I can change this last thing myself if you prefer)

CERN / TIGRE

Fix and refactor checkDevices #304

TEST CODE