numpy -- Transform non-contiguous data to contiguous data in place

Question

Consider the following code:

import numpy as np
a = np.zeros(50)
a[10:20:2] = 1
b = c = a[10:40:4]
print b.flags  # You'll see that b and c are not C_CONTIGUOUS or F_CONTIGUOUS

My question:

Is there a way (with only a reference to b) to make both b and c contiguous? It is completely fine if np.may_share_memory(b,a) returns False after this operation.

Things which are close, but don't quite work out are: np.ascontiguousarray/np.asfortranarray as they will return a new array.

---

My use case is that I have very large 3D fields stored in a subclass of a numpy.ndarray. In order to save memory, I would like to chop those fields down to the portion of the domain that I am actually interested in processing:

a = a[ix1:ix2,iy1:iy2,iz1:iz2]

Slicing for the subclass is somewhat more restricted than slicing of ndarray objects, but this should work and it will "do the right thing" -- the various custom meta-data attached on the subclass will be transformed/preserved as expected. Unfortunately, since this returns a view, numpy won't free the big array afterward so I don't actually save any memory here.

To be completely clear, I'm looking to accomplish 2 things:

• preserve the metadata on my class instance. slicing will work, but I'm not sure about other forms of copying.
• make it so that the original array is free to be garbage collected

Answer 1

According to Alex Martelli:

"The only really reliable way to ensure that a large but temporary use of memory DOES return all resources to the system when it's done, is to have that use happen in a subprocess, which does the memory-hungry work then terminates."

However, the following appears to free at least some of the memory: Warning: my way of measuring free memory is Linux-specific:

import time
import numpy as np

def free_memory():
    """
    Return free memory available, including buffer and cached memory
    """
    total = 0
    with open('/proc/meminfo', 'r') as f:
        for line in f:
            line = line.strip()
            if any(line.startswith(field) for field in ('MemFree', 'Buffers', 'Cached')):
                field, amount, unit = line.split()
                amount = int(amount)
                if unit != 'kB':
                    raise ValueError(
                        'Unknown unit {u!r} in /proc/meminfo'.format(u=unit))
                total += amount
    return total

def gen_change_in_memory():
    """
    https://stackoverflow.com/a/14446011/190597 (unutbu)
    """
    f = free_memory()
    diff = 0
    while True:
        yield diff
        f2 = free_memory()
        diff = f - f2
        f = f2
change_in_memory = gen_change_in_memory().next

Before allocating the large array:

print(change_in_memory())
# 0

a = np.zeros(500000)
a[10:20:2] = 1
b = c = a[10:40:4]

After allocating the large array:

print(change_in_memory())
# 3844 # KiB

a[:len(b)] = b
b = a[:len(b)]
a.resize(len(b), refcheck=0)
time.sleep(1)

Free memory increases after resizing:

print(change_in_memory())
# -3708 # KiB

Answer 2

You can do this in cython:

In [1]:
%load_ext cythonmagic

In [2]:
%%cython
cimport numpy as np

np.import_array()

def to_c_contiguous(np.ndarray a):
    cdef np.ndarray new
    cdef int dim, i
    new = a.copy()
    dim = np.PyArray_NDIM(new)
    for i in range(dim):
        np.PyArray_STRIDES(a)[i] = np.PyArray_STRIDES(new)[i]
    a.data = new.data
    np.PyArray_UpdateFlags(a, np.NPY_C_CONTIGUOUS)
    np.set_array_base(a, new)

In [8]:
import sys
import numpy as np
a = np.random.rand(10, 10, 10)
b = c = a[::2, 1::3, 2::4]
d = a[::2, 1::3, 2::4]
print sys.getrefcount(a)
to_c_contiguous(b)
print sys.getrefcount(a)
print np.all(b==d)

The output is:

4
3
True

to_c_contiguous(a) will create a c_contiguous copy of a, and make it as the base of a.

After the call of to_c_contiguous(b), the refcount of a is decreased, and when the refcount of a become 0, it will be freed.

Answer 3

I would claim the correct way to accomplish the 2 things you listed is by np.copying the slices you create.

Of course, in order for that to work correctly, you would need to define an appropriate __array_finalize__. You weren't very clear about why you decided to avoid it in the first place, but my feeling is that you should define it. (how did you solve the bx**2 problem without using __array_finalize__?)