Hope you are well. I've been converting an adjacency matrix to a connectivity index. I am having problems trying to get it take the second neighbour, but have been having issues. I will try to explain my problem and goal as best as I can. It will be better to explain by an example.
This is an example of an adjacency matrix.
A0 A1 A2 A3 A4 A5
A0 0.0 1.0 0.0 0.0 0.0 0.0
A1 1.0 0.0 1.0 0.0 0.0 0.0
A2 0.0 1.0 0.0 1.0 0.0 0.0
A3 0.0 0.0 1.0 0.0 1.0 0.0
A4 0.0 0.0 0.0 1.0 0.0 1.0
A5 0.0 0.0 0.0 0.0 1.0 0.0
The titles (A0, A1 etc.) represents an atom. The values represents a connection. So a value of 1 in the row A0 and column A1 means there is a connection between the 2 atoms and they are neighbours.
I am following an equation to find the first order (first neighbour) connectivity index (CF) by the following equation:
CF = Σ(Si * Sj)^(-0.5)
Where Si corresponds to the connectivity of the atom in row i (as in the example above it would be A0) and Sj for j (A1 in above example). Connectivity (S) is defined as the number of bonds, so the sum of the row or column of that atom.
I used the code below to apply this.
def randic_index(A):
n = A.shape[0]
R = 0
for i in range(n):
for j in range(n):
if A[i,j] != 0:
deg_i = np.sum(A[i,:])
deg_j = np.sum(A[:,j])
R += 1 / np.sqrt(deg_i * deg_j)
return 0.5*R
This is all good for the first order connectivity index, but my issues are in retrieving the second order (second neighbour). I have tried different methods, but none seem to be working for me.
The formula for the second order connectivity index (CS) is:
CS = Σ(Si * Sj * Sk)^(-0.5) where k corresponds to the neighbour of the j atom.
I need to code a method of selecting the neighbours of the first neighbour and reiterating for all the potential second neighbours,leaving no combination that hasn't been calculated.
Is it possible anyone has advice on how to tackle this, or provide a modification to solve this? And, perhaps, provide something that will allow higher orders than second to be calculated, though I expect I would be able to adjust that myself.
Many Thanks
If the sum for the second order index is taken over all triples (i, j, k) of connected nodes, including the ones where i=k, then the following should work:
import numpy as np
def randic_index2_path(A):
degs = 1 / A.sum(axis=0)**0.5
B = (A * degs**0.5 * degs[..., None])
return np.tril((B.T@B)).sum()
If, on the other hand, the sum is taken over triples of connected nodes (i, j, k) where i and k are different, then you can try the following:
def randic_index2_trail(A):
degs = 1 / A.sum(axis=0)**0.5
B = (A * degs**0.5 * degs[..., None])
return np.tril((B.T@B), k=-1).sum()
Also, the function computing the first order index included in the question can be rewritten as follows:
def randic_index(A):
degs = 1 / A.sum(axis=0)**0.5
return np.tril(A * degs * degs[..., None]).sum()