python - Broadcasting with reduction or extension in Numpy -

in following code calculate magnitudes of vectors between pairs of given points. speed operation in numpy can use broadcasting

import numpy np points = np.random.rand(10,3)  pair_vectors = points[:,np.newaxis,:] - points[np.newaxis,:,:] pair_dists = np.linalg.norm(pair_vectors,axis=2).shape 

or outer product iteration

it = np.nditer([points,points,none], flags=['external_loop'], op_axes=[[0,-1,1],[-1,0,1],none]) a,b,c in it:     c[...] = b - pair_vectors = it.operands[2] pair_dists = np.linalg.norm(pair_vectors,axis=2) 

my question how 1 use broadcasting or outer product iteration create array form 10x10x6 last axis contains coordinates of both points in pair (extension). , in related way, possible calculate pair distances using broadcasting or outer product iteration directly, i.e. produce matrix of form 10x10 without first calculating difference vectors (reduction).

to clarify, following code creates desired matrices using slow looping.

pair_coords = np.zeros(10,10,6) pair_dists = np.zeros(10,10) in range(10):     j in range(10):         pair_coords[i,j,0:3] = points[i,:]         pair_coords[i,j,3:6] = points[j,:]         pair_dists[i,j] = np.linalg.norm(points[i,:]-points[j,:]) 

this failed attempt calculate distanced (or apply other function takes 6 coordinates of both points in pair , produce scalar) using outer product iteration.

res = np.zeros((10,10)) = np.nditer([points,points,res], flags=['reduce_ok','external_loop'], op_axes=[[0,-1,1],[-1,0,1],none]) a,b,c in it: c[...] = np.linalg.norm(b-a) pair_dists = it.operands[2] 

here's approach produce arrays in vectorized ways -

from itertools import product scipy.spatial.distance import pdist, squareform  n = points.shape[0]  # indices selecting rows off points array , stacking them  idx = np.array(list(product(range(n),repeat=2))) p_coords = np.column_stack((points[idx[:,0]],points[idx[:,1]])).reshape(n,n,6)  # distances upper triangular elements.  # create symmetric 1 final dists array. p_dists = squareform(pdist(points)) 

few other vectorized approaches discussed in this post, have there too!