Restriction with repeated indices

[mrava87]

Problem

The adjoint of the Restriction operator is broken when iava has repeated values and is applied to a CuPy array, eg:

from pylops import Restriction

R = Restriction(10, np.array([0,2, 2]))
y = R * np.ones(10)
R.H * y
>> ValueError: operands could not be broadcast together with shape (3,) and requested shape (2,)

This is due to the fact that the adjoint is implemented differently for Numpy/Cupy with the latter using self.iavamask. The problem is in _compute_iavamask at

iavamask[iava] = 1
...
iavamask = np.where(iavamask.ravel() == 1)[0]

since repeated indices are found only once.

However, whilst this works for NumPy the operator does not pass the dot-test anymore as the adjoint instead of summing the repeated values it just overwrites the last one.

Proposal

Add values of repeated indices in adjoint. For both Numpy/Cupy, it is possible to simply use

ncp.add.at(
     y,
     tuple(
           self.iava if ax == self.axis else slice(None)
            for ax in range(y.ndim)
     ),
    x,
)