TomographicImaging / CIL

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

VectorData/ImageData mulitplication #1977

Open epapoutsellis opened 14 hours ago

epapoutsellis commented 14 hours ago

@bosschmidt

There is a problem with ImageData multiplied with VectorData

from cil.framework import VectorGeometry, ImageGeometry
vg = VectorGeometry(505)
ig  = ImageGeometry(505,505)

x1 = ig.allocate("random")
x2 = vg.allocate("random")

res1 = x1.array.dot(x2.array)
print(res1.shape)

res2 = x1.dot(x2)

returns

(505,)

ValueError                                Traceback (most recent call last)
Cell In[3], line 11
      8 res1 = x1.array.dot(x2.array)
      9 print(res1.shape)
---> 11 res2 = x1.dot(x2)

File /opt/conda/lib/python3.11/site-packages/cil/framework/data_container.py:718, in DataContainer.dot(self, other, *args, **kwargs)
    716         return sf
    717 else:
--> 718     raise ValueError('Shapes are not aligned: {} != {}'.format(self.shape, other.shape))

ValueError: Shapes are not aligned: (505, 505) != (505,)

So when comparing shapes in https://github.com/TomographicImaging/CIL/blob/master/Wrappers/Python/cil/framework/data_container.py#L705

we should check self.shape[1] == other.shape[0].

jakobsj commented 9 hours ago

Thanks for spotting this. Agree this can be seen as an inconsistency with but not sure if it is a bug though.

Numpy .dot when called with a 2D array (such as x1.array) does matrix-matrix multiplication (here matrix-vector since x2.array is 1D), see https://numpy.org/devdocs/reference/generated/numpy.dot.html which also suggests that .dot is not preferred for matrix-matrix multiplication, instead recommends matmul or @.

CIL datacontainer's .dot considers each element a vector and so determines that they have different length and thus cannot compute their inner product, which I think is our intended behaviour. Forcing the same behaviour of datacontainer's .dot would imply that a datacontainer can represent a matrix, something that is generally represented by LinearOperator.

It seems there is a discrepancy between numpy's convention of .dot and the mathematical differences between a vector and linear operator and it is not clear that we'd want to align with numpy's convention here. Can I ask what you were trying to achieve by carrying out the dot operation?