code: def _cartesian(arraySizes, out=None):
      
      arrays = []
      for i in range(len(arraySizes)):
          arrays.append(nu.arange(0, arraySizes[i]))
  
      arrays = [nu.asarray(x) for x in arrays]
      dtype = arrays[0].dtype
     
      n = nu.prod([x.size for x in arrays])
      if out is None:
          out = nu.zeros([n, len(arrays)], dtype=dtype)
      
      m = n // arrays[0].size
      out[:,0] = nu.repeat(arrays[0], m)
      if arrays[1:]:
          _cartesian(arraySizes[1:], out=out[0:m,1:])
          for j in range(1, arrays[0].size):
              out[j*m:(j+1)*m,1:] = out[0:m,1:]
      return out
docs: NAME:
          cartesian
      PURPOSE:
          Generate a cartesian product of input arrays.
      INPUT: 
          arraySizes - list of size of arrays
          out - Array to place the cartesian product in.
      OUTPUT: 
          2-D array of shape (product(arraySizes), len(arraySizes)) containing cartesian products
          formed of input arrays.
      HISTORY:
          2016-06-02 - Obtained from
          http://stackoverflow.com/questions/1208118/using-numpy-to-build-an-array-of-all-combinations-of-two-arrays