Define setter for Image.array to allow inplace operations

[sidneymau]

(Duplicate of #1273 due to rebasing/etc.)
This pull request addresses #1272 by defining a setter for the Image.array attribute. This protects against any case where a user does an inplace operation on the image array attribute in an interactive environment: without the setter, this would raise an error while still performing the operation, leading to some confusion. This change also obviates the need to operate on a view of the underling array; e.g., instead of

image.array[:] += 1

we can now do

image.array += 1

[sidneymau]

Copying comment from Mike:

I'd like a few more tests that guarantee that the setter doesn't allocate new memory or change fundamental things about the array. I think the implementation is ok, but there are other ways to have written the setter which would pass this test, but would break the image object.
e.g.

image1.array = 1  # should keep array shape and dtype
image1.array = 1+2j  # should only work if array is already complex
image1.array = np.zeros((3,3))  # should only work if current shape is (3,3)

Also, all of those should fail if the image is a cont view.
In other words, try to come up with as many stress tests as possible and make sure it behaves as we would want it to behave.

[sidneymau]

An alternative implementation of the setter is the following:

    @array.setter
    def array(self, other):
        """Set the numpy array in-place.
        """
        self_array = self.array
        other_array = np.asarray(other)
        if self_array.shape != other_array.shape:
            raise GalSimError("Other array shape is not equal to current array shape")
        if self_array.dtype != other_array.dtype:
            raise GalSimError("Other array dtype is not equal to current array dtype")
        self._array = other_array

This has an advantage in explicitly guarding against mutations to shape and/or dtype. I believe this is more robust and closer to what we want

[rmjarvis]

This looks worse to me. I think numpy already guarded against these things. And you lost the feature that writing to a const view would raise an exception.

[rmjarvis]

Also, the dtype check is wrong. For instance, if dtype is float, you should be allowed to assign integers to it and let it remain a float array. (Which I believe the previous implementation did correctly.)

[rmjarvis]

Oh, another problem with this implementation, which might merit a test is that it loses the careful 128 bit alignment that we do for image arrays. The array should always satisfy:

assert image.array.ctypes.data % 16 == 0

Looks like we don't currently have tests for this, but it makes the SSE stuff in C++ work faster.

[sidneymau]

I suppose there's a third version: self._array = np.asarray(other) without the type checking:

>>> image = galsim.Image(np.zeros((2, 2), dtype=int))
>>> image.array = np.array([[1.2, 3.8], [5.6, 6.7]])
>>> image.array
array([[1.2, 3.8],
       [5.6, 6.7]])

This is probably preferred to each of the above? Unless we want to enforce that the dtype of an Image array doesn't mutate.

We would still want to test against shapes, I think, otherwise this leads to

>>> image.array = 1
>>> image.array
array(1)

[sidneymau]

Oh, another problem with this implementation, which might merit a test is that it loses the careful 128 bit alignment that we do for image arrays. The array should always satisfy:

assert image.array.ctypes.data % 16 == 0

Looks like we don't currently have tests for this, but it makes the SSE stuff in C++ work faster.

There's currently a commented assert checking this elsewhere in Image—is this worth raising an error/warning for in the class, or just writing a test?

followup: I have not succeeded in making a numpy array that does not satisfy array.ctypes.data % 16 == 0, so will have to figure out how to do that to write this test

[rmjarvis]

I have not succeeded in making a numpy array that does not satisfy array.ctypes.data % 16 == 0, so will have to figure out how to do that to write this test

Maybe numpy added that as a feature. They didn't used to guarantee that.

[rmjarvis]

>>> image.array = 1
>>> image.array
array(1)

This should set all values to 1, which your original implementation did correctly.

[rmjarvis]

>>> image = galsim.Image(np.zeros((2, 2), dtype=int))
>>> image.array = np.array([[1.2, 3.8], [5.6, 6.7]])
>>> image.array
array([[1.2, 3.8],
       [5.6, 6.7]])

This one is definitely wrong. Assignment should absolutely not change the dtype. To be consistent with how things like += work currently, I think assigning e.g. float to int should just convert them. The implementation would be

self.array[:] = self._safe_cast(other)

cf. methods like __iadd__ etc. in the Image class.

[sidneymau]

Should this be a warning or an error? We could also try realigning the array, but I'm not sure if that's within scope here...

[rmjarvis]

I don't think you're getting my point. I don't want the array object to change upon setting. Use the one that is already in the image object. The original implementation you had was basically right. Here is another test that I think might make my point clearer.

im = galsim.ImageI(3,3)
ar = np.ones((3,3))
im.array = ar   # This test currently works if use im.array[:] = ar instead
assert np.all(im.array == 1)
ar *= 2
assert np.all(im.array == 1)

In other words, IMO the setter should work exactly like if you used [:].

The only exception is that floats should round rather than truncate to match how +=, etc. work.

Another probably useful test is that if the image starts out all 0's, then im.array = ar should be equivalent to im.array += ar. For instance, the above test should also work if you use replace the assignment with

im.array += ar

This is close to the original reason for this PR -- you had wanted to do that += step, and it gave an error, but still worked.

[sidneymau]

These will break if numpy is not guaranteeing that arrays are memory-aligned by default, which I gather was not previously the case. We can get around this by writing a function to explicitly make memory-aligned arrays, but doing so may not be necessary

[rmjarvis]

The code that used to be necessary, but apparently might not be anymore is here:
https://github.com/GalSim-developers/GalSim/blob/releases/2.5/galsim/image.py#L593

[sidneymau]

I'm not sure what the preference for broadcasting Image arrays is—this tests that an int array gets successfully cast to a float array, which may not be the desired behavior

[rmjarvis]

It is not desired behavior.

[rmjarvis]

>>> image.array = 1
>>> image.array
array(1)

This should set all values to 1, which your original implementation did correctly.

[rmjarvis]

>>> image = galsim.Image(np.zeros((2, 2), dtype=int))
>>> image.array = np.array([[1.2, 3.8], [5.6, 6.7]])
>>> image.array
array([[1.2, 3.8],
       [5.6, 6.7]])

This one is definitely wrong. Assignment should absolutely not change the dtype. To be consistent with how things like += work currently, I think assigning e.g. float to int should just convert them. The implementation would be

self.array[:] = self._safe_cast(other)

cf. methods like __iadd__ etc. in the Image class.

[rmjarvis]

I don't think you're getting my point. I don't want the array object to change upon setting. Use the one that is already in the image object. The original implementation you had was basically right. Here is another test that I think might make my point clearer.

im = galsim.ImageI(3,3)
ar = np.ones((3,3))
im.array = ar   # This test currently works if use im.array[:] = ar instead
assert np.all(im.array == 1)
ar *= 2
assert np.all(im.array == 1)

In other words, IMO the setter should work exactly like if you used [:].

The only exception is that floats should round rather than truncate to match how +=, etc. work.

Another probably useful test is that if the image starts out all 0's, then im.array = ar should be equivalent to im.array += ar. For instance, the above test should also work if you use replace the assignment with

im.array += ar

This is close to the original reason for this PR -- you had wanted to do that += step, and it gave an error, but still worked.

[rmjarvis]

The code that used to be necessary, but apparently might not be anymore is here:
https://github.com/GalSim-developers/GalSim/blob/releases/2.5/galsim/image.py#L593

[rmjarvis]

It is not desired behavior.

[rmjarvis]

@sidneymau This has been dormant for a couple months now. Did you want to try to finish it, or should we abandon this effort?

[sidneymau]

Thanks for bumping this—I'm happy to consider this abandoned for now and have the PR closed and the branch pruned. I may pick it up again when I have more time, but I don't expect that to be any time soon.