getkeops / keops

KErnel OPerationS, on CPUs and GPUs, with autodiff and without memory overflows
https://www.kernel-operations.io
MIT License
1.04k stars 64 forks source link

Broadcasting problem with complex ? #187

Open chloesrcb opened 3 years ago

chloesrcb commented 3 years ago

It seems that, in both pyKeOps and RKeOps, the use of Real2Complex doesn't support broadcasting when using binary operations like addition. For instance, we can't add a single complex value (e.g. 2 + 3i) to a vector or a matrix, complex or not.

Is that made on purpose ?

With RKeOps:

# works
formula <- "Sum_Reduction(Add(v,s),1)"
args <- c("v=Pm(6)", "s=Pm(1)")
op1 <- keops_kernel(formula, args)

# works
formula <- "Sum_Reduction(Add(v,Real2Complex(s)),1)"
args <- c("v=Pm(6)", "s=Pm(3)")
op3 <- keops_kernel(formula, args)
# doesn't work
formula <- "Sum_Reduction(Add(v,Real2Complex(s)),1)"
args <- c("v=Pm(6)", "s=Pm(1)")
op2 <- keops_kernel(formula, args)
# Error in compile_formula(formula, var_aliases$var_aliases, dllname) : 
#     Error during cmake call. 

Detailed error:

static_assert( ((F::DIM==DIM)||(F::DIM==1)) &&
((G::DIM==DIM)||(G::DIM==1)) , "Incompatible dimensions for vectorized scalar binary operation.");

static_assert(DIM == FB::DIM, "Dimensions must be the same for Add");

With LazyTensors in pyKeOps:

import numpy as np
from pykeops.numpy import Vi, Vj, Pm

Cplx = 2 + 3j
v = [2 + 3j, 7 + 3j, 4 + 1j]
x = np.random.rand(10, 2)

Pm_Scal = Pm(3.14)
Pm_Cplx = Pm(Cplx)
Pm_v = Pm(v)
x_i = Vi(x)

print(Pm_v + Pm_Scal)  # ValueError: incompatible shapes for addition.
print(Pm_v + Pm_Cplx)  # ValueError: incompatible shapes for addition.
print(x_i + Pm_Cplx)  # ValueError: incompatible shapes for addition.

Amélie Vernay and Chloé Serre-Combe

joanglaunes commented 2 years ago

Hello @AmelieVernay , Thanks for noticing about this. We already discussed about this via email, so I am just repeating here what we said : this is indeed a bug / missing feature, and it has been fixed in the branch python_engine, which is planned to be merged into master soon. We will do the checks as soon as the branch will be merged and the R bindings updated accordingly. So I am keeping this issue open for now.

AmelieVernay commented 2 years ago

Ok, thankyou!

Chloé and Amélie

twmitchel commented 2 years ago

@joanglaunes

I'm encountering a more general problem in pyKeOps - ComplexLazyTensor broadcasting is failing entirely. For example

a_i = ComplexLazyTensor(torch.rand(1, 3, 1, 1, 1, 11, 4, dtype=torch.cfloat))
b_i = ComplexLazyTensor(torch.rand(2, 1, 4, 4, 10, 1, 1, dtype=torch.cfloat))

print(a_i * b_i, flush=True)

fails for me.

Also, if I were to print the sizes of a_i and b_i, the sizes of the last dimensions would be 8 and 2, respectively, which correspond to the dimensions referenced in the error message.

Unfortunately, I cannot incorporate KeOps into my project without complex number support (even though it would be perfect).

Has this been addressed already in the python_engine branch?

If so, is there a way I could install that branch directly?

Otherwise, is there any time frame on when we can expect the master branch to be updated with the changes in the python_engine branch?

joanglaunes commented 2 years ago

Hello @twmitchel , In fact it looks like the bug you are mentioning had been addressed previously in the master branch, because it appears to work in master, as well as in python_engine, whereas I can reproduce the bug in the current release. You can install the python_engine branch as follows: pip install git+https://github.com/getkeops/keops.git@python_engine (also same without @python_engine to install master) The python_engine branch is supposed to be ready now, we just need to do more testing (so in fact getting your feedback about this branch will be great !)