TomographicImaging / CIL

A versatile python framework for tomographic imaging
https://tomographicimaging.github.io/CIL/
Apache License 2.0
91 stars 41 forks source link

BlockOperator.domain_geometry().allocate() not compatible with in place calls to BlockOperator.direct #1863

Open MargaretDuff opened 2 weeks ago

MargaretDuff commented 2 weeks ago

Description

Due to #1802, if you use out=myblockoperator.domain_geometry().allocate(0) where myblockoperator is a BlockOperator of size (1,1) then out will be a DataContainer and not a BlockDataContainer.

If you then call myblockoperator.direct(something, out=out) you get an error because you try and call getitem on the data container out.

This leads to errors when you try and run PDHG or SPDHG on partitioned data with just one batch e.g.


data = dataexample.SIMPLE_PHANTOM_2D.get(size=(10, 10))

        subsets = 1

        ig = data.geometry
        ig.voxel_size_x = 0.1
        ig.voxel_size_y = 0.1

        detectors = ig.shape[0]
        angles = np.linspace(0, np.pi, 90)
        ag = AcquisitionGeometry.create_Parallel2D().set_angles(
            angles, angle_unit='radian').set_panel(detectors, 0.1)
        # Select device
        dev = 'cpu'

        Aop = ProjectionOperator(ig, ag, dev)

        sin = Aop.direct(data)
        partitioned_data = sin.partition(subsets, 'sequential')
        A = BlockOperator(
            *[IdentityOperator(partitioned_data[i].geometry) for i in range(subsets)])

        # block function
        F = BlockFunction(*[L2NormSquared(b=partitioned_data[i])
                                 for i in range(subsets)])
        alpha = 0.025
        G = alpha * FGP_TV()

        spdhg = SPDHG(f=F, g=G, operator=A,  update_objective_interval=10)

        spdhg.run(7)   
        pdhg = PDHG(f=F, g=G, operator=A, update_objective_interval=10)

        pdhg.run(7)       
MargaretDuff commented 2 weeks ago

@paskino - would appreciate your help on this one, it is a bit mind-boggling