Predictinglinks.py

import numpy as np
from scipy.linalg import eigvalsh
from scipy.optimize import minimize

def structural_hamiltonian(A, beta):
    """
    Calculate the structural Hamiltonian H(A).
    
    Args:
    A (np.array): Adjacency matrix
    beta (np.array): Parameters β_k
    
    Returns:
    float: H(A) value
    """
    eigenvalues = eigvalsh(A)
    H = 0
    for k in range(3, len(beta) + 3):
        H -= beta[k-3] * np.log(np.sum(eigenvalues**k))
    return H

def conditional_probability(A, x, y, beta):
    """
    Calculate the conditional probability of link (x,y).
    
    Args:
    A (np.array): Adjacency matrix
    x, y (int): Indices of the link
    beta (np.array): Parameters β_k
    
    Returns:
    float: Conditional probability
    """
    A_xy = A.copy()
    A_xy[x, y] = A_xy[y, x] = 1
    delta_H = structural_hamiltonian(A_xy, beta) - structural_hamiltonian(A, beta)
    return 1 / (1 + np.exp(delta_H))

def log_pseudo_likelihood(beta, A):
    """
    Calculate the log pseudo-likelihood.
    
    Args:
    beta (np.array): Parameters β_k
    A (np.array): Adjacency matrix
    
    Returns:
    float: Negative log pseudo-likelihood (for minimization)
    """
    N = A.shape[0]
    ll = 0
    for x in range(N):
        for y in range(x+1, N):
            p = conditional_probability(A, x, y, beta)
            ll += A[x, y] * np.log(p) + (1 - A[x, y]) * np.log(1 - p)
    return -ll

def estimate_parameters(A, initial_beta):
    """
    Estimate parameters β_k using maximum pseudo-likelihood.
    
    Args:
    A (np.array): Adjacency matrix
    initial_beta (np.array): Initial guess for β_k
    
    Returns:
    np.array: Estimated β_k
    """
    result = minimize(lambda beta: log_pseudo_likelihood(beta, A), 
                      initial_beta, method='L-BFGS-B')
    return result.x

def calculate_link_scores(A, beta):
    """
    Calculate scores S_xy for all possible links.
    
    Args:
    A (np.array): Adjacency matrix
    beta (np.array): Estimated parameters β_k
    
    Returns:
    dict: Scores for each possible link
    """
    N = A.shape[0]
    scores = {}
    for x in range(N):
        for y in range(x+1, N):
            if A[x, y] == 0:  # Only for non-observed links
                scores[(x, y)] = conditional_probability(A, x, y, beta)
    return scores

# Example usage
N = 10  # Number of nodes
A = np.random.randint(0, 2, size=(N, N))
A = (A + A.T) // 2  # Make it symmetric
np.fill_diagonal(A, 0)  # No self-loops

initial_beta = np.zeros(3)  # Assuming we're using k=3,4,5
estimated_beta = estimate_parameters(A, initial_beta)
print("Estimated parameters:", estimated_beta)

link_scores = calculate_link_scores(A, estimated_beta)
print("Link scores for non-observed links:", link_scores)