search

SyneRBI / SIRF

Main repository for the CCP SynerBI software
http://www.ccpsynerbi.ac.uk
Other
58 stars 29 forks source link

Incosistency between dot and get_inner_product #629

Open AnderBiguri opened 4 years ago

AnderBiguri commented 4 years ago

If you use: 

examples_path = SIRF_PATH + '/data/examples/Registration'
x = sirf.Reg.NiftiImageData3D(rexamples_path + "/test.nii.gz") #ref_aladin 
T = na.get_transformation_matrix_forward()
y = sirf.Reg.NiftiImageData3D(examples_path + "/test2.nii.gz") #flo_aladin

nr = sirf.Reg.NiftyResample()
nr.set_reference_image(x)
nr.set_floating_image(y)
nr.set_interpolation_type_to_linear()
nr.add_transformation(T)
# Add this on test_pReg.py#L950 to test.

# NOTE: you need PR #626 for this to work. Otherwise shaw domain/range
x = nr.domain_geometry().allocate(value = 'random')
y = nr.range_geometry().allocate( value = 'random')
y_hat = nr.direct(x)

res = y_hat.dot(y)

It returns a NiftiImageData

However, if you instead use

res = y_hat.get_inner_product(y)

The output variable is type float. Slightly confusing, expectations would be that it returns some DataContainer class, or at least that dot and get_inner_product return the same class. This last expectation is stronger, particularly because the way SIRF uses dot is basically by defining it as an inner product (e.g. it also does complex dot products). Mathematical definitions aside, one would expect that for NiftyReg data/operators, dot and get_inner_product to be exactly the same function.

I suspect that this may come from having multiple libraries with different functions and trying to harmonize all of them, however while I personally think that with leaving them as it is goes against user expectations (here I am again telling you what to do 😉!), at least there should be more docs about the output (otherwise you expect users of SIRF to know the source code).

KrisThielemans commented 4 years ago

This is related to our discussion of dot vs vdot. @AnderBiguri can you insert the link to that here (ideally very briefly summarising if necessary). (There's numpy.inner as well by the way, which isn't what we want either).

An inner_product has to return a scalar. Obviously, for complex data, it'd have to be a complex number, while for real data a float (or double).

I guess it's less clear what dot should return, but I agree with @AnderBiguri that if in SIRF (v)dot is equivalent to get_inner_product, then it has to return the same thing.

AnderBiguri commented 4 years ago

There was a discussion at https://github.com/SyneRBI/SIRF/pull/617#discussion_r404330231

The conversation is about naming of dot. The current SIRF behavior of dot is:

It treats the input as a vector (unrolls any matrix), and uses complex conjugate when data is complex.

Arguments are

  1. If we want SIRF to be as equal as possible between python and MATLAB, then it makes sense to use dot as a name. However, MATLAB does not unroll the input, it treats different column as separate dot product inputs (vectorization).

  2. If we want SIRF to be consistent with numpy, then we should call it vdot. But this breaks uniformity with the MATLAB version.

So the question is here still, dot or vdot? and what about get_inner_product? I am good with any. Calling it some version of dot is consistent with computer scientist, calling it inner_product is consistent with mathematicians ( and I am neither! :) ). I guess the important issue is that we have now both in SIRF, and they return different types. All a bit confusing.

Please add any point I may be missing.

KrisThielemans commented 4 years ago

thanks @AnderBiguri. I'll rope in some of the CIL guys.

@epapoutsellis, @jakobsj, @paskino I understand that CIL uses numpy.dot. Are all CIL data containers 1D then? I guess they are real-valued.

My personal opinion is that if there is confusion in terminology and functionality between different *dot versions, we should stay away from it and stick to (get_)inner_product.

evgueni-ovtchinnikov commented 4 years ago

As a mathematician, I very much regret the meaning of dot being totally distorted by extending this function, which in mathematics is only defined for vectors, to matrices, for which it now returns just the matrix product in Python and products. of matrices' columns or rows in MATLAB (the latter is at least useful in coding block iterative algorithms).

Luckily, in mathematics we have another name for this operation on vectors: inner product - so, if we cannot stay with our current dot, then I suggest we replace it with inner (get_inner_product is much too ugly for my taste).

KrisThielemans commented 4 years ago

I agree. The confusion comes that in imaging, the images are multiple-dimension arrays, but they are not operators, and not matrices either.

I also agree that get_inner_product is unwieldy, but would avoid inner as it conflicts with numpy.inner again. So I'd vote for inner_product.

In fact, I'd then make a sirf.inner_product function that calls a.inner_product(b) (I believe that in MATLAB, inner_product(a,b) will automatically try a.inner_product(b) anyway)

evgueni-ovtchinnikov commented 4 years ago

there is no inner_product in MATLAB R2019b

KrisThielemans commented 4 years ago

there is no inner_product in MATLAB R2019b

no, exactly. that's why we won't have a name clash!

I meant that MATLAB will go and hunt for a method if the function doesn't exist, see https://uk.mathworks.com/help/matlab/matlab_oop/method-invocation.html