BlockDataContainers with (1,1) shape

[samdporter]

I'm currently doing some reconstructions where I end up with BlockDataContainers with shape=(1,1).

A 1xn BlockOperator needs to take a nx1 BlockDataContainer and output a 1x1 DataContainer so that this can be used in a CIL function (KullbackLeibler in this case).

Unfortunately, this results in a 1x1 BlockDataContainer, which causes an error due to the lack of an as_array() method.

Would it be possible to add a method such as

 def __new__(cls, *args, **kwargs):
    '''
    Create a new BlockDataContainer object 
    or return a DataContainer if the shape is (1,1)
    '''
    shape = kwargs.get('shape', None)
    if shape is None:
        shape = (len(args),1)
    if shape == (1,1):
        return DataContainer(args[0].as_array(), *args, **kwargs)
    else:
        return super(BlockDataContainer, cls).__new__(cls)`

So that a 1x1 `BlockDataContainer` always results in a `DataContainer

Of course, I might just be using the Block Framework in the wrong way so any help would be much appreciated. I'll comment some example code below

[paskino]

Thanks @samdporter, it looks like the direct method of BlockOperator is the culprit as it always returns a BlockDataContainer.

https://github.com/TomographicImaging/CIL/blob/45e0b235ac79efc03dc5cb46b8ddda9d16f6331b/Wrappers/Python/cil/optimisation/operators/BlockOperator.py#L181

We could fix it as done in adjoint where in the case of a 1x1 we just return the data container: https://github.com/TomographicImaging/CIL/blob/45e0b235ac79efc03dc5cb46b8ddda9d16f6331b/Wrappers/Python/cil/optimisation/operators/BlockOperator.py#L227-L231

Your suggestion is more elegant, but I'm not sure about the consequences it may have.

 def __new__(cls, *args, **kwargs):
    '''
    Create a new BlockDataContainer object 
    or return a DataContainer if the shape is (1,1)
    '''
    shape = kwargs.get('shape', None)
    if shape is None:
        shape = (len(args),1)
    if shape == (1,1):
        return args[0]
    else:
        return super(BlockDataContainer, cls).__new__(cls)

[samdporter]

Hey @paskino,

Thanks for looking at this

This did have some unintended consequences, unfortunately, so I stopped using it in my code and instead added a as_array() Function for a (1,1) BlockDataContainer which is working so far

[MargaretDuff]

Discussed with @samdporter! We should make this a priority before the SIRF challenge and the SIRF/CIL training at PSMR

[MargaretDuff]

@paskino suggested looking at his PR https://github.com/TomographicImaging/CIL/pull/1766

[epapoutsellis]

I think that BlockOperator needs a BlockFunction , so KullbackLeibler should be inside a BlockFunction. Could you share some code? If you end up with (1,1) BlockDataContainer then you may not need a BlockOperator.

[samdporter]

I'm doing some synergistic deconvolution. I'm working with SIRF which means I need a BlockDataContainerfor multiple channels. The problem comes when I define a operator like this:

k = op.BlockOperator(operator0, operator0, shape = (1,2))

I want this to take a BDC with two images and return one image - operator0(image0) + operator1(image1) but instead it returns a (1,1) BDC which SIRF doesn't recognise.

PR 1766nearly fixes this but only returns a DataContainer for the adjoint method if necessary

This is what I'm using it for

[epapoutsellis]

@samdporter can you do

spdhg_f1(k1.direct(image_dict['OSEM']))

if shape=(1,2), then you are in a product space, $X{1}\times X{2}$, so you need (2,1) vector, but the domain of spdhg_f1 is only one, e.g., $X_{1}$. What is the shape of initial=initial_estimate? If it is a DataContainer then the problem is with domain/range of operators in the BlockOperator and the corresponding domains of functions.

[samdporter]

initial_estimate = BlockDataContainer(image0, image1) but spdhg_f1 = KullbackLeibler(b=data0) needs a DataContainer. This is required so that image0, image1 can have separate data fidelity functions but share the same prior.

[MargaretDuff]

I'm doing some synergistic deconvolution. I'm working with SIRF which means I need a BlockDataContainerfor multiple channels. The problem comes when I define a operator like this:

k = op.BlockOperator(operator0, operator0, shape = (1,2))

I want this to take a BDC with two images and return one image - operator0(image0) + operator1(image1) but instead it returns a (1,1) BDC which SIRF doesn't recognise.

PR 1766nearly fixes this but only returns a DataContainer for the adjoint method if necessary

This is what I'm using it for

Hi @samdporter - this issue is popping it's head up again - could you point out where in the code it was causing an issue?