[FIXED] Vectorizing calculating the hydrodynamic radius in numpy

Issue

I have a polymer with coordinates stored in Nx3 numpy array, where N is the number of particles in the polymer (degree of polymerization).

I am trying to calculate the hydrodynamic radius of this polymer. The hydrodynamic radius is given by the first expression found in this link. The hydrodynamic radius, Rh is essentially a harmonic average over pair-wise distance.

Given that P is the Nx3 array, this is my current numpy-pythonic implementation:

inv_dist = 0
for i in range(N-1):
    for j in range(i+1, N): 
        inv_dist += 1/(np.linalg.norm (P[i,:]-P[j,:], 2))


Rh = 1/(inv_dist/(N**2) )

np is numpy in this case. I am aware that the wikipedia formula asks for an ensemble average. This would mean that I would loop over EVERY possible configuration of the polymer in my simulation. In any event, the two loops mentioned above will still be computed.

This is a nested for loop with Nx(N-1)/2 iterations. As N gets large, this computation becomes increasingly taxing. How can I vectorize this code to bypass the for loops to an extent?

I would appreciate any advice you have for me.

Solution

You can use scipy.spatial.distance.pdist:

from scipy.spatial.distance import pdist

inv_dist = (1/pdist(P)).sum()

Answered By - user2246849