Calculate value of double summation of function over the pairs of vertices of a graph

Question

Given a directed graph G with node pair (p,q) we have

I want to calculate the value of this recursive function where L(p) denotes the set of undirected link neighbors of node p. This function is for the (k+1)th value.

I know how to calculate L(p) and L(q).

Here is my attempt-

from __future__ import division
import copy

import numpy as np
import scipy
import scipy.sparse as sp
import time
import warnings




class Algo():
    def __init__(self, c=0.8, max_iter=5, is_verbose=False, eps=1e-3):
        self.c = c
        self.max_iter = max_iter
        self.is_verbose = is_verbose
        self.eps = eps

    def compute_similarity(self, a):
        '''

        Note: The adjacency matrix is a (0,1)-matrix with zeros on its 
        diagonal.
        '''                                 
        a_transpose=np.transpose(a)
        a_undirected=np.add(a,a_transpose)

        self.num_nodes = np.shape(a)[0]
        current_sim = np.eye(self.num_nodes)
        prev_sim = np.zeros(shape=(self.num_nodes, self.num_nodes))

        #Determine the set of Neighbors for each node
        nbs = []

        for i in range(self.num_nodes):

            nbs.append(np.nonzero(a_undirected[:, i])[0])


        for itr in range(self.max_iter):

            if scipy.allclose(current_sim, prev_sim, atol=self.eps):
                break

            prev_sim = copy.deepcopy(current_sim)

            for i in range(self.num_nodes):
                for j in range(self.num_nodes):
                    if i == j:
                        continue

                    if len(nbs[i]) == 0 or len(nbs[j]) == 0:
                        current_sim[i][j] = 0
                    else:
                        #commonSet=0.0
                        differenceSetofp_ij=0.0
                        differenceSetofq_ij=0.0
                        commonSet=len(np.intersect1d(nbs[i],nbs[j]))/len(np.union1d(nbs[i],nbs[j]))


                        for x in np.setdiff1d(nbs[i],nbs[j]):
                            for y in nbs[j]:

                                differenceSetofp_ij+=prev_sim[x][y]
                                differenceSetofp_ij*=1/(len(np.union1d(nbs[i],nbs[j]))*len(nbs[j]))
                                #print differenceSetofp



                        for x in np.setdiff1d(nbs[j],nbs[i]):
                            for y in nbs[i]:

                                differenceSetofq_ij+=prev_sim[x][y]
                                differenceSetofq_ij*=1/(len(np.union1d(nbs[j],nbs[i]))*len(nbs[i])) 

                        current_sim[i][j] = self.c*(commonSet+differenceSetofp_ij+differenceSetofq_ij)


        if self.is_verbose:

            print('SimRank: converged after {0} iterations'.format(itr))
        return current_sim                  



if __name__ == "__main__":

    a = np.array([ [0, 1, 1, 0, 0, 1, 0, 0, 0],
                   [0, 0, 0, 0, 1, 0, 0, 0, 0],
                   [0, 0, 0, 1, 0, 0, 0, 0, 0],
                   [0, 0, 0, 0, 0, 0, 0, 0, 0],
                   [0, 0, 0, 0, 0, 0, 1, 0, 1],
                   [0, 0, 0, 1, 0, 0, 0, 0, 0],
                   [0, 0, 0, 0, 0, 1, 0, 0, 0],
                   [0, 0, 0, 0, 0, 0, 1, 0, 1],
                   [0, 0, 0, 0, 0, 0, 0, 0, 0]])                    



    obj = Algo(is_verbose=True, eps=1e-2)
    start_time = time.time()
    S = obj.compute_similarity(a)
    print(S)
    print('Runtime: ', time.time() - start_time, '\n')  

The output I am getting is not correct.