memoryerror while creating cartesian product in numpy

#shared dimensions:
   sh = a.shape[1: ]
aba = (a[: , None, None] + b[None,: , None] - a[None, None,: ]).reshape(-1, * sh)
aba
#array([
   [0., 0.03],
   #[-1., 0.16],
   #[1., -0.1],
   #[0., 0.03]
])

    [
       [2],
       [
          [7],
          [
             [2],
             7 + [3], equiv to[7], +[3],
             [4]
          ][7]
       ][4]
    ]

              [
                 [4],
                 [
                    [1, 2, 3],
                    [
                       [4, 4, 4],
                       [
                          [1, 2, 3]
                       ] + [5], equiv to[1, 2, 3], +[5, 5, 5],
                       [6]
                    ][1, 2, 3]
                 ][6, 6, 6]
              ]

idx = cartesian_product(
   np.arange(len(a)),
   np.arange(len(b)) + len(a),
   np.arange(len(a))
)

x = np.concatenate((a, b))
result = x[idx.ravel(),: ].reshape( * idx.shape, -1)

size_in_GB = A.shape[0] ** 2 * A.shape[1] * B.shape[0] * (size_of_datatype) / 1e9

import numpy as np
import tables #register blosc
import h5py as h5
import h5py_cache as h5c

a = np.arange(500 * 50).reshape(500, 50)
b = np.arange(40 * 50).reshape(40, 50)

# isn 't well documented, have a look at https://github.com/Blosc/hdf5-blosc
compression_opts = (0, 0, 0, 0, 5, 1, 1)
compression_opts[4] = 9 #compression level 0.. .9
compression_opts[5] = 1 #shuffle
compression_opts[6] = 1 #compressor(I guess that 's lz4)

      File_Name_HDF5 = 'Test.h5'
      f = h5.File(File_Name_HDF5, 'w', chunk_cache_mem_size = 1024 ** 2 * 300) dset = f.create_dataset('Data', shape = (a.shape[0] ** 2 * b.shape[0], a.shape[1]), dtype = 'd', chunks = (a.shape[0] * b.shape[0], 1), compression = 32001, compression_opts = (0, 0, 0, 0, 9, 1, 1), shuffle = False)

      #Write the data
      for i in range(a.shape[0]):
      sh = a.shape[1: ] aba = (a[i] + b[: , None] - a).reshape(-1, * sh) dset[i * a.shape[0] * b.shape[0]: (i + 1) * a.shape[0] * b.shape[0]] = aba

      f.close()

File_Name_HDF5 = 'Test.h5'
f = h5c.File(File_Name_HDF5, 'r', chunk_cache_mem_size = 1024 ** 2 * 300)
dset = f['Data']
chunks_size = 500
for i in range(0, dset.shape[0], chunks_size):
   #Iterate over the first column
data = dset[i: i + chunks_size,: ] #avoid excessive calls to the hdf5 library
#Do something with the data

f.close()

f = h5c.File(File_Name_HDF5, 'r', chunk_cache_mem_size = 1024 ** 2 * 300)
dset = f['Data']
for i in range(dset.shape[1]):
   # Iterate over the second dimension
# fancy indexing e.g.[: , i] will be much slower
# use np.expand_dims or in this
case np.squeeze after the read operation from the dset
#
if you wan 't to have the same result than [:,i] (1 dim array)
data = dset[: , i: i + 1]
#Do something with the data

f.close()

def cartesian_2d(arrays, out = None): arrays = [np.asarray(x) for x in arrays] dtype = arrays[0].dtype n = np.prod([x.shape[0]
   for x in arrays
]) if out is None: out = np.empty([n, len(arrays), arrays[0].shape[1]], dtype = dtype) m = n // arrays[0].shape[0]     out[:, 0] = np.repeat(arrays[0], m, axis=0)     if arrays[1:]:         cartesian_2d(arrays[1:], out=out[0:m, 1:, :])         for j in range(1, arrays[0].shape[0]):             out[j * m:(j + 1) * m, 1:] = out[0:m, 1:]     return out a = [[ 0, -0.02], [1, -0.15]] b = [[0, 0.03]] result = cartesian_2d([a,b,a]) // array([[[ 0.  , -0.02],     [ 0.  ,  0.03],     [ 0.  , -0.02]],    [[ 0.  , -0.02],     [ 0.  ,  0.03],     [ 1.  , -0.15]],    [[ 1.  , -0.15],     [ 0.  ,  0.03],     [ 0.  , -0.02]],    [[ 1.  , -0.15],     [ 0.  ,  0.03],       [ 1.  , -0.15]]]) 

def cartesian_2d(arrays, out = None): arrays = [np.asarray(x) for x in arrays] dtype = arrays[0].dtypen = np.prod([x.shape[0]
   for x in arrays
]) if out is None: out = np.empty([n, len(arrays), arrays[0].shape[1]], dtype = dtype) m = n // arrays[0].shape[0]out[:, 0] = np.repeat(arrays[0], m, axis=0)if arrays[1:]:  cartesian_2d(arrays[1:], out=out[0:m, 1:, :])  for j in range(1, arrays[0].shape[0]): out[j * m:(j + 1) * m, 1:] = out[0:m, 1:]return out  a = [[ 0, -0.02], [1, -0.15]] b = [[0, 0.03]]  result = cartesian_2d([a,b,a])  // array([[[ 0.  , -0.02],[ 0.  ,  0.03],[ 0.  , -0.02]],[[ 0.  , -0.02],[ 0.  ,  0.03],[ 1.  , -0.15]],[[ 1.  , -0.15],[ 0.  ,  0.03],[ 0.  , -0.02]],[[ 1.  , -0.15],[ 0.  ,  0.03],[ 1.  , -0.15]]]) 

result[: , 0,: ] + result[: , 1,: ] - result[: , 2,: ] //array([[ 0.  ,  0.03], [-1.  ,  0.16], [ 1.  , -0.1 ], [ 0.  ,  0.03]]) 

#shared dimensions: sh = a.shape[1: ] aba = (a[: , None, None] + b[None,: , None] - a[None, None,: ]).reshape(-1, * sh) aba #array([
   [0., 0.03], #[-1., 0.16], #[1., -0.1], #[0., 0.03]
])

    [
       [2],
       [
          [7],
          [
             [2], 7 + [3], equiv to[7], +[3],
             [4]
          ][7]
       ][4]
    ]

  [
     [4],
     [
        [1, 2, 3],
        [
           [4, 4, 4],
           [
              [1, 2, 3]
           ] + [5], equiv to[1, 2, 3], +[5, 5, 5],
           [6]
        ][1, 2, 3]
     ][6, 6, 6]
  ]

The product method in itertools provides an implementation of the Cartesian product that when run on with many arguments quickly gives out of memory errors.The current implementation creates a lot of unnecessary lists in this situation.A more appropriate implementation uses dynamic programming to avoid these out of memory issues.

The trouble is that itertools.product accepts iterators, and there is no guaranteed way of "restarting"
an arbitrary iterator in Python.Consider:

   >>>
   a = iter([1, 2, 3]) >>>
   b = iter([4, 5, 6]) >>>
   next(a)
1
   >>>
   next(b)
4
   >>>
   from itertools
import product
   >>>
   list(product(a, b))[(2, 5), (2, 6), (3, 5), (3, 6)]

Since there 's no way to get back to items you'
ve already consumed, the current approach is to consume all of the iterators to begin with and store their items in arrays, then lazily produce tuples of the items at the right indices of those arrays.

Perhaps one could consume lazily from the iterators, say, only filling up the pools as they 're needed and not storing the contents of the first iterator, but this would mean sometimes producing a product iterator that was doomed to cause a memory error eventually. If you really need this behavior you could do this in Python:

def lazy_product( * iterables):
   if not iterables:
   yield()
return
it0 = iterables[0]
for x in it0:
   print(f "{x=}")
for rest in lazy_product( * iterables[1: ]):
   print(f "{rest=}")
yield(x, ) + rest

The above could surely be optimized as maybe you 're suggesting, but this would be a backward-incompatible change for itertools.product.

Possibly related:

   https: //bugs.python.org/issue10109

   Henry, I 'm not clear at all about what you'
re saying.Please give at least one specific, concrete example of the behavior you 're objecting to, and specify the behavior you want to see instead.

Hey, Tim, I just wanted a note of clarification.I was working on an approximation algorithm
for a network science tool that is being released soon.Initially, I relied on itertools.product(), but when I moved to testing on non - trivial graphs, I immediately had Out of Memory Errors.This was even on the nodes with large memories on the computing cluster that I was using.The issue is a lot of extra lists are added as the number of iterables passed product grows although the actual number of elements in the iterables is relatively small.
Here is the workaround that I used to handle many arguments:

   def product( * args) - > Iterable:
   ''
'
Takes in several iterables and returns a generator
that yields the cartesian product of the elements in the iterables: param args: sets which you would like the cartesian product of
   : return: generator ''
'
if len(args) == 1 and type(args) == list:
   return args
numArgs = len(args)
indices = np.zeros(numArgs, dtype = np.uint)
checkLen = len(args[0])
while indices[0] < checkLen:
   curr_term = [args[0][indices[0]]]
for i in range(1, numArgs):
   curr_term.append(args[i][indices[i]])
indices[-1] += np.uint(1)
for i in range(1, numArgs):
   updated = False
if indices[-i] >= len(args[-i]):
   indices[-i] = np.uint(0)
indices[-i - 1] += np.uint(1)
updated = True
if not updated:
   break
yield curr_term

Henry, I have to ask again: please give at least one specific, concrete example of the behavior you 're objecting to.  English isn'
t helping, and I still have no idea what your complaint is.

What I 'm looking for:

for i in itertools.product( ?? ? ):
   pass

where you replace the ?? ? with executable code(preferably not using numpy or any other extension package) that provokes the MemoryError you 're talking about.

For example, here I 'm constructing a case with a million arguments.  There'
s no problem at all:

   >>>
   import itertools >>>
   args = [range(100) for i in range(1_000_000)] >>>
   for i in itertools.product( * args):
   ...print(len(i))[and it prints 1000000 over & over until I get bored and kill it]

Note
if it _did_ provoke a problem, we probably wouldn 't care - there are no plausible realistic use cases for passing a million arguments to product().

You haven 't given us a way to reproduce your complaint, or even a clue as to the number of arguments you'
re talking about.The above code was my best guess as to what you _might_ be talking about.But since there 's no MemoryError in sight, apparently not.

I 'm suspecting this may be an XY problem[1], and especially because you posted "a solution" instead of an actual problem ;-)

[1] https: //meta.stackexchange.com/questions/66377/what-is-the-xy-problem

> when I moved to testing on non - trivial graphs, I immediately had >
   Out of Memory Errors.

I 'm going to hazard a guess that the input to product() was a graph traversal iterator that got trapped in an undetected cycle.  Feeding an infinite iterator into product() would eat all available memory, just as it would with list().

>>> np.zeros(10) //1
array([-0., -0., -0., -0., -0., -0., -0., -0., -0., -0.])

>>> np.zeros(10) //1
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])

npy_intp
const fixed_dims[] = {
   1,
   2,
   3
};
// no longer complains that the const-qualifier is discarded
npy_intp size = PyArray_MultiplyList(fixed_dims, 3);

NPY_BLAS_ORDER = openblas python setup.py build

NPY_LAPACK_ORDER = openblas python setup.py build

def foo():
   pass