Repository: NCAS-CMS/cfdm

Issue:
Incorrect data assignment with some multiple list indices
When multiple list indices are provided to a data assignment (i.e. `Data.__setitem__`), the assignment works if all of the lists are strictly monotonically increasing, but only for some cases cases where they are not:

```python
>>> import cfdm
>>> d = cfdm.example_field(0).data
>>> d[[0, 3, 4], [1, 7, 6]].array 
array([[0.034, 0.029, 0.024],
       [0.059, 0.017, 0.009],
       [0.036, 0.013, 0.034]])
>>> d[[0, 3, 4], [1, 7, 6]] = -1       # all lists strictly monotonically increasing => OK
>>> d[[0, 3, 4], [1, 7, 6]].array 
array([[-1., -1., -1.],
       [-1., -1., -1.],
       [-1., -1., -1.]])
```

```python
>>> import cfdm
>>> d = cfdm.example_field(0).data
>>> d[[0, 3, 4], [1, 6, 7]].array 
array([[0.034, 0.029, 0.024],
       [0.059, 0.017, 0.009],
       [0.036, 0.013, 0.034]])
>>> d[[0, 3, 4], [1, 7, 6]] = -1 
>>> d[[0, 3, 4], [1, 7, 6]].array        # one list not strictly monotonically increasing => OK (edited from "not OK"!)
array([[-1., -1., -1.],
       [-1., -1., -1.],
       [-1., -1., -1.]])
```

```python
>>> d = cfdm.example_field(0).data
>>> d[[4, 3, 0], [7, 6, 1]] = -1
>>> d[[4, 3, 0], [7, 6, 1]].array       # both lists strictly monotonically decreasing => NOT OK
array([[ 0.013,  0.034,  0.036],
       [ 0.017,  0.009,  0.059],
       [ 0.029,  0.024, -1.   ]])
```
```python
>>> d = cfdm.example_field(0).data
>>> d[[4, 3, 0], [7, 6, 1]] = -1
>>> d[[4, 3, 0], [7, 6, 1]].array     # NOT OK
array([[ 0.013,  0.034,  0.036],
       [ 0.017,  0.009,  0.059],
       [ 0.029,  0.024, -1.   ]])
```
```python
>>> d = cfdm.example_field(0).data
>>> d[[0, 4, 3], [1, 7, 6]] = -1
>>> d[[0, 4, 3], [1, 7, 6]].array   # OK, though!
array([[-1., -1., -1.],
       [-1., -1., -1.],
       [-1., -1., -1.]])
```

PR to follow (after cfdm has been transferred to `main`) (edit: https://github.com/davidhassell/cfdm/tree/data-setitem-multiple-list-indices-2, lest I forget)

```python
>>> cfdm.environment(paths=False)
Platform: Linux-5.14.0-1052-oem-x86_64-with-glibc2.31
HDF5 library: 1.12.1
netcdf library: 4.8.1
Python: 3.9.12
netCDF4: 1.6.0
numpy: 1.22.3
cfdm.core: 1.10.0.0
cftime: 1.6.1
netcdf_flattener: 1.2.0
cfdm: 1.10.0.0
```

Discussion:

Patch:
diff --git a/Changelog.rst b/Changelog.rst
index 4c24fc3aa..232c07d6e 100644
--- a/Changelog.rst
+++ b/Changelog.rst
@@ -10,6 +10,8 @@ Version 1.10.0.1
   (https://github.com/NCAS-CMS/cfdm/issues/215)
 * New method: `cfdm.Field.get_original_filenames`
 * New method: `cfdm.Data.get_original_filenames`
+* Fixed bug that caused incorrect data assignment with some multiple
+  list indices (https://github.com/NCAS-CMS/cfdm/issues/217)
 * Fixed bug that caused a failure when printing date-time data with
   the first element masked
   (https://github.com/NCAS-CMS/cfdm/issues/211)
diff --git a/cfdm/cfdmimplementation.py b/cfdm/cfdmimplementation.py
index e14be1d41..a119c88f7 100644
--- a/cfdm/cfdmimplementation.py
+++ b/cfdm/cfdmimplementation.py
@@ -1625,22 +1625,6 @@ class CFDMImplementation(Implementation):
         """
         return cell_measure.get_measure(default=None)
 
-    def get_original_filenames(self, parent):
-        """Get the original names of the files containing the construct.
-
-        :Parameters:
-
-            parent: Construct
-                The construct to query of original filenames.
-
-        :Returns:
-
-            `set`
-                The original filenames
-
-        """
-        return parent.get_original_filenames()
-
     def nc_get_dimension(self, parent, default=None):
         """Return the netCDF variable name.
 
diff --git a/cfdm/data/data.py b/cfdm/data/data.py
index 676ed6945..9699dec6e 100644
--- a/cfdm/data/data.py
+++ b/cfdm/data/data.py
@@ -858,8 +858,83 @@ class Data(Container, NetCDFHDF5, Files, core.Data):
         self._set_Array(array, copy=copy)
 
     @classmethod
-    def _set_subspace(cls, array, indices, value):
-        """Set a subspace of the data array defined by indices."""
+    def _set_subspace(cls, array, indices, value, orthogonal_indexing=True):
+        """Assign to a subspace of an array.
+
+        :Parameters:
+
+            array: array_like
+                The array to be assigned to. Must support
+                `numpy`-style indexing. The array is changed in-place.
+
+            indices: sequence
+                The indices to be applied.
+
+            value: array_like
+                The value being assigned. Must support fancy indexing.
+
+            orthogonal_indexing: `bool`, optional
+                If True then apply 'orthogonal indexing', for which
+                indices that are 1-d arrays or lists subspace along
+                each dimension independently. This behaviour is
+                similar to Fortran but different to, for instance,
+                `numpy` or `dask`.
+
+        :Returns:
+
+            `None`
+
+        **Examples**
+
+        Note that ``a`` is redefined for each example, as it is
+        changed in-place.
+
+        >>> a = np.arange(40).reshape(5, 8)
+        >>> {{package}}.Data._set_subspace(a, [[1, 4 ,3], [7, 6, 1]],
+        ...                    np.array([[-1, -2, -3]]))
+        >>> print(a)
+        [[ 0  1  2  3  4  5  6  7]
+         [ 8 -3 10 11 12 13 -2 -1]
+         [16 17 18 19 20 21 22 23]
+         [24 -3 26 27 28 29 -2 -1]
+         [32 -3 34 35 36 37 -2 -1]]
+
+        >>> a = np.arange(40).reshape(5, 8)
+        >>> {{package}}.Data._set_subspace(a, [[1, 4 ,3], [7, 6, 1]],
+        ...                    np.array([[-1, -2, -3]]),
+        ...                    orthogonal_indexing=False)
+        >>> print(a)
+        [[ 0  1  2  3  4  5  6  7]
+         [ 8  9 10 11 12 13 14 -1]
+         [16 17 18 19 20 21 22 23]
+         [24 -3 26 27 28 29 30 31]
+         [32 33 34 35 36 37 -2 39]]
+
+        >>> a = np.arange(40).reshape(5, 8)
+        >>> value = np.linspace(-1, -9, 9).reshape(3, 3)
+        >>> print(value)
+        [[-1. -2. -3.]
+         [-4. -5. -6.]
+         [-7. -8. -9.]]
+        >>> {{package}}.Data._set_subspace(a, [[4, 4 ,1], [7, 6, 1]], value)
+        >>> print(a)
+        [[ 0  1  2  3  4  5  6  7]
+         [ 8 -9 10 11 12 13 -8 -7]
+         [16 17 18 19 20 21 22 23]
+         [24 25 26 27 28 29 30 31]
+         [32 -6 34 35 36 37 -5 -4]]
+
+        """
+        if not orthogonal_indexing:
+            # --------------------------------------------------------
+            # Apply non-orthogonal indexing
+            # --------------------------------------------------------
+            array[tuple(indices)] = value
+            return
+
+        # ------------------------------------------------------------
+        # Still here? Then apply orthogonal indexing
+        # ------------------------------------------------------------
         axes_with_list_indices = [
             i
             for i, x in enumerate(indices)
@@ -867,55 +942,118 @@ class Data(Container, NetCDFHDF5, Files, core.Data):
         ]
 
         if len(axes_with_list_indices) < 2:
-            # --------------------------------------------------------
             # At most one axis has a list-of-integers index so we can
-            # do a normal numpy assignment
-            # --------------------------------------------------------
+            # do a normal assignment
             array[tuple(indices)] = value
         else:
-            # --------------------------------------------------------
             # At least two axes have list-of-integers indices so we
-            # can't do a normal numpy assignment
-            # --------------------------------------------------------
+            # can't do a normal assignment.
+            #
+            # The brute-force approach would be to do a separate
+            # assignment to each set of elements of 'array' that are
+            # defined by every possible combination of the integers
+            # defined by the two index lists.
+            #
+            # For example, if the input 'indices' are ([1, 2, 4, 5],
+            # slice(0:10), [8, 9]) then the brute-force approach would
+            # be to do 4*2=8 separate assignments of 10 elements each.
+            #
+            # This can be reduced by a factor of ~2 per axis that has
+            # list indices if we convert it to a sequence of "size 2"
+            # slices (with a "size 1" slice at the end if there are an
+            # odd number of list elements).
+            #
+            # In the above example, the input list index [1, 2, 4, 5]
+            # can be mapped to two slices: slice(1,3,1), slice(4,6,1);
+            # the input list index [8, 9] is mapped to slice(8,10,1)
+            # and only 2 separate assignments of 40 elements each are
+            # needed.
             indices1 = indices[:]
-            for i, x in enumerate(indices):
+            for i, (x, size) in enumerate(zip(indices, array.shape)):
                 if i in axes_with_list_indices:
-                    # This index is a list of integers
+                    # This index is a list (or similar) of integers
+                    if not isinstance(x, list):
+                        x = np.asanyarray(x).tolist()
+
                     y = []
                     args = [iter(x)] * 2
                     for start, stop in itertools.zip_longest(*args):
-                        if not stop:
+                        if start < 0:
+                            start += size
+
+                        if stop is None:
                             y.append(slice(start, start + 1))
+                            break
+
+                        if stop < 0:
+                            stop += size
+
+                        step = stop - start
+                        if not step:
+                            # (*) There is a repeated index in
+                            #     positions 2N and 2N+1 (N>=0). Store
+                            #     this as a single-element list
+                            #     instead of a "size 2" slice, mainly
+                            #     as an indicator that a special index
+                            #     to 'value' might need to be
+                            #     created. See below, where this
+                            #     comment is referenced.
+                            #
+                            #     For example, the input list index
+                            #     [1, 4, 4, 4, 6, 2, 7] will be mapped
+                            #     to slice(1,5,3), [4], slice(6,1,-4),
+                            #     slice(7,8,1)
+                            y.append([start])
                         else:
-                            step = stop - start
-                            stop += 1
+                            if step > 0:
+                                stop += 1
+                            else:
+                                stop -= 1
+
                             y.append(slice(start, stop, step))
 
                     indices1[i] = y
                 else:
                     indices1[i] = (x,)
 
-            if numpy.size(value) == 1:
+            if value.size == 1:
+                # 'value' is logically scalar => simply assign it to
+                # all index combinations.
                 for i in itertools.product(*indices1):
                     array[i] = value
-
             else:
+                # 'value' has two or more elements => for each index
+                # combination for 'array' assign the corresponding
+                # part of 'value'.
                 indices2 = []
-                ndim_difference = array.ndim - numpy.ndim(value)
-                for i, n in enumerate(numpy.shape(value)):
-                    if n == 1:
+                ndim_difference = array.ndim - value.ndim
+                for i2, size in enumerate(value.shape):
+                    i1 = i2 + ndim_difference
+                    if i1 not in axes_with_list_indices:
+                        # The input 'indices[i1]' is a slice
                         indices2.append((slice(None),))
-                    elif i + ndim_difference in axes_with_list_indices:
+                        continue
+
+                    index1 = indices1[i1]
+                    if size == 1:
+                        indices2.append((slice(None),) * len(index1))
+                    else:
                         y = []
                         start = 0
-                        while start < n:
+                        for index in index1:
                             stop = start + 2
+                            if isinstance(index, list):
+                                # Two consecutive elements of 'value'
+                                # are assigned to the same integer
+                                # index of 'array'.
+                                #
+                                # See the (*) comment above.
+                                start += 1
+
                             y.append(slice(start, stop))
                             start = stop
 
                         indices2.append(y)
-                    else:
-                        indices2.append((slice(None),))
 
                 for i, j in zip(
                     itertools.product(*indices1), itertools.product(*indices2)
diff --git a/cfdm/read_write/netcdf/netcdfwrite.py b/cfdm/read_write/netcdf/netcdfwrite.py
index d902671b7..c6990ee16 100644
--- a/cfdm/read_write/netcdf/netcdfwrite.py
+++ b/cfdm/read_write/netcdf/netcdfwrite.py
@@ -4288,12 +4288,8 @@ class NetCDFWrite(IOWrite):
                 As for the *format* parameter for initialising a
                 `netCDF.Dataset` instance.
 
-            fields: sequence of `Field` or `Domain`
-                The constructs to be written to the netCDF file. Note
-                that these constructs are only used to ascertain if
-                any data to be written is in *filename*. If this is
-                the case and mode is "w" then an exception is raised
-                to prevent *filename* from being deleted.
+            fields: sequence of `Field`
+                The field constructs to be written.
 
         :Returns:
 
@@ -4301,13 +4297,13 @@ class NetCDFWrite(IOWrite):
                 A `netCDF4.Dataset` object for the file.
 
         """
-        if fields and mode == "w":
+        if fields and mode != "r":
             filename = os.path.abspath(filename)
             for f in fields:
-                if filename in self.implementation.get_original_filenames(f):
+                if filename in self.implementation.get_filenames(f):
                     raise ValueError(
-                        "Can't write with mode 'w' to a file that contains "
-                        f"data that needs to be read: {f!r} uses {filename}"
+                        "Can't write to a file that contains data "
+                        f"that needs to be read: {filename}"
                     )
 
         # mode == 'w' is safer than != 'a' in case of a typo (the letters


Test patch:
diff --git a/cfdm/test/test_Data.py b/cfdm/test/test_Data.py
index 291c3b495..940a0d585 100644
--- a/cfdm/test/test_Data.py
+++ b/cfdm/test/test_Data.py
@@ -6,6 +6,7 @@ import os
 import unittest
 
 import numpy
+import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
 
@@ -75,10 +76,9 @@ class DataTest(unittest.TestCase):
         self.assertTrue(str(d) == "[--, 2000-01-21 00:00:00]")
         self.assertTrue(repr(d) == "<Data(2): [--, 2000-01-21 00:00:00]>")
 
-    #    def test_Data__getitem__(self):
     def test_Data__setitem__(self):
         """Test the assignment of data items on Data."""
-        a = numpy.ma.arange(3000).reshape(50, 60)
+        a = np.ma.arange(3000).reshape(50, 60)
 
         d = cfdm.Data(a.filled(), units="m")
 
@@ -95,47 +95,138 @@ class DataTest(unittest.TestCase):
             n = -n - 1
             for dvalue, avalue in (
                 (n, n),
-                (cfdm.masked, numpy.ma.masked),
+                (cfdm.masked, np.ma.masked),
                 (n, n),
             ):
-                message = f"cfdm.Data[{j}, {i}]={dvalue}={avalue} failed"
                 d[j, i] = dvalue
                 a[j, i] = avalue
                 x = d.array
-                self.assertTrue(
-                    (x == a).all() in (True, numpy.ma.masked), message
-                )
-                m = numpy.ma.getmaskarray(x)
-                self.assertTrue(
-                    (m == numpy.ma.getmaskarray(a)).all(),
-                    "d.mask.array="
-                    + repr(m)
-                    + "\nnumpy.ma.getmaskarray(a)="
-                    + repr(numpy.ma.getmaskarray(a)),
-                )
+                self.assertTrue((x == a).all() in (True, np.ma.masked))
+                m = np.ma.getmaskarray(x)
+                self.assertTrue((m == np.ma.getmaskarray(a)).all())
 
-        a = numpy.ma.arange(3000).reshape(50, 60)
+        a = np.ma.arange(3000).reshape(50, 60)
 
         d = cfdm.Data(a.filled(), "m")
 
         (j, i) = (slice(0, 2), slice(0, 3))
-        array = numpy.array([[1, 2, 6], [3, 4, 5]]) * -1
+        array = np.array([[1, 2, 6], [3, 4, 5]]) * -1
 
-        for dvalue in (array, numpy.ma.masked_where(array < -2, array), array):
-            message = "cfdm.Data[%s, %s]=%s failed" % (j, i, dvalue)
+        for dvalue in (array, np.ma.masked_where(array < -2, array), array):
             d[j, i] = dvalue
             a[j, i] = dvalue
             x = d.array
-            self.assertTrue((x == a).all() in (True, numpy.ma.masked), message)
-            m = numpy.ma.getmaskarray(x)
-            self.assertTrue((m == numpy.ma.getmaskarray(a)).all(), message)
+            self.assertTrue((x == a).all() in (True, np.ma.masked))
+            m = np.ma.getmaskarray(x)
+            self.assertTrue((m == np.ma.getmaskarray(a)).all())
 
         # Scalar numeric array
         d = cfdm.Data(9, units="km")
         d[...] = cfdm.masked
         a = d.array
         self.assertEqual(a.shape, ())
-        self.assertIs(a[()], numpy.ma.masked)
+        self.assertIs(a[()], np.ma.masked)
+
+        # Multiple list indices, scalar value
+        d = cfdm.Data(np.arange(40).reshape(5, 8), units="km")
+
+        value = -1
+        for indices in (
+            ([0, 3, 4], [1, 6, 7]),
+            ([0, 3, 4], [1, 7, 6]),
+            ([0, 4, 3], [1, 6, 7]),
+            ([0, 4, 3], [1, 7, 6]),
+            ([4, 3, 0], [7, 6, 1]),
+            ([4, 3, 0], [1, 6, 7]),
+            ([0, 3, 4], [7, 6, 1]),
+            ([0, 3, -1], [7, 6, 1]),
+            ([0, 3, 4], [-1, 6, 1]),
+            ([0, 3, -1], [-1, 6, 1]),
+        ):
+            d[indices] = value
+            self.assertEqual((d.array == value).sum(), 9)
+            value -= 1
+
+        # Repeated list elements
+        for indices in (
+            ([0, 3, 3], [7, 6, 1]),
+            ([3, 3, 0], [7, 6, 1]),
+            ([0, 4, 3], [7, 6, 7]),
+        ):
+            d[indices] = value
+            self.assertEqual((d.array == value).sum(), 6)
+            value -= 1
+
+        for indices in (
+            ([3, 4, 3], [7, 6, 7]),
+            ([3, 3, 4], [7, 7, 6]),
+            ([4, 3, 3], [6, 7, 7]),
+        ):
+            d[indices] = value
+            self.assertEqual((d.array == value).sum(), 4)
+            value -= 1
+
+        # Multiple list indices, array value
+        a = np.arange(40).reshape(1, 5, 8)
+
+        value = np.arange(9).reshape(3, 3) - 9
+
+        for indices in (
+            (slice(None), [0, 3, 4], [1, 6, 7]),
+            (slice(None), [0, 3, 4], [1, 7, 6]),
+            (slice(None), [0, 4, 3], [1, 6, 7]),
+            (slice(None), [0, 4, 3], [1, 7, 6]),
+            (slice(None), [4, 3, 0], [7, 6, 1]),
+            (slice(None), [4, 3, 0], [1, 6, 7]),
+            (slice(None), [0, 3, 4], [7, 6, 1]),
+            (slice(None), [0, 3, -1], [7, 6, 1]),
+            (slice(None), [0, 3, 4], [-1, 6, 1]),
+            (slice(None), [0, 3, -1], [-1, 6, 1]),
+        ):
+            d = cfdm.Data(a.copy())
+            d[indices] = value
+            self.assertEqual((d.array < 0).sum(), 9)
+
+        # Repeated list elements
+        for indices in (
+            (slice(None), [0, 3, 3], [7, 6, 1]),
+            (slice(None), [0, 4, 3], [7, 6, 7]),
+            (slice(None), [3, 3, 4], [1, 6, 7]),
+            (slice(None), [0, 4, 3], [7, 7, 6]),
+        ):
+            d = cfdm.Data(a.copy())
+            d[indices] = value
+            self.assertEqual((d.array < 0).sum(), 6)
+
+        for indices in (
+            (slice(None), [3, 4, 3], [7, 6, 7]),
+            (slice(None), [4, 3, 3], [6, 7, 7]),
+            (slice(None), [3, 3, 4], [6, 7, 7]),
+            (slice(None), [3, 3, 4], [7, 7, 6]),
+            (slice(None), [4, 3, 3], [7, 7, 6]),
+        ):
+            d = cfdm.Data(a.copy())
+            d[indices] = value
+            self.assertEqual((d.array < 0).sum(), 4)
+
+        # Multiple list indices, array value + broadcasting
+        value = np.arange(3).reshape(1, 3) - 9
+
+        for indices in ((slice(None), [0, 3, 4], [1, 6, 7]),):
+            d = cfdm.Data(a.copy())
+            d[indices] = value
+            self.assertEqual((d.array < 0).sum(), 9)
+
+        # Repeated list elements
+        for indices in ((slice(None), [0, 3, 3], [7, 6, 1]),):
+            d = cfdm.Data(a.copy())
+            d[indices] = value
+            self.assertEqual((d.array < 0).sum(), 6)
+
+        for indices in ((slice(None), [4, 3, 3], [7, 7, 6]),):
+            d = cfdm.Data(a.copy())
+            d[indices] = value
+            self.assertEqual((d.array < 0).sum(), 4)
 
     def test_Data_apply_masking(self):
         """Test the `apply_masking` Data method."""