Coding Isomap (& MDS) function using only numpy and scipy in python

Question

I have coded Isomap function starting with computing the eulidean distance matrix (using scipy.spatial.distance.cdist), next basing on K-nearest neighbors method and Dijkstra algorithm (to determinate the shortest path) I have Computed the full distance matrix over all paths, finally I have did map computations, following by the dimensionality reduction. BUT, I want to use epsilon instead of K-nearest neighbors like in the following :

Y = isomap (X, epsilon, d)

• X is an n × m matrix which corresponds to n points with m attributes.

• epsilon is an anonymous function of the distance matrix used to find the parameters of neighborhood. (The neighborhood graph must be formed by eliminating the edges whose width is greater to epsilon of the complete distance graph).

• d is a parameter which signifies the output dimension.

• Y is an n × d matrix, which signifies the embedding resulting from isomap.

THANKS in advance

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist


def distance_Matrix(X):
    return cdist(X,X,'euclidean')

def Dijkstra(h):
    q = h.copy()
    for i in range(ndata):
        for j in range(ndata):
            k = np.argmin(q[i,:])
            while not(np.isinf(q[i,k])):
                q[i,k] = np.inf
                for l in neighbours[k,:]:
                    possible = h[i,l] + h[l,k]
                    if possible < h[i,k]:
                        h[i,k] = possible
                k = np.argmin(q[i,:])
    return h

def MDS(D,newdim=2):  

    n = D.shape[0]
    # Torgerson formula
    I = np.eye(n)
    J = np.ones(D.shape)
    J = I-(1/n)*J
    B = (-1/2)*np.dot(np.dot(J,D),np.dot(D,J))  # B = -(1/2).JD²J

    # 
    eigenval, eigenvec = np.linalg.eig(B)
    indices = np.argsort(eigenval)[::-1]   
    eigenval = eigenval[indices]
    eigenvec = eigenvec[:, indices]

    # dimension reduction
    K = eigenvec[:, :newdim]
    L = np.diag(eigenval[:newdim])  
    # result
    Y = K @ L **(1/2)    
    return np.real(Y)

def isomap(data,newdim=2,K=12):

    ndata = np.shape(data)[0]
    ndim = np.shape(data)[1]

    d = distance_Matrix(X)

    # replace begin 
    # K-nearest neighbours
    indices = d.argsort()
    #notneighbours = indices[:,K+1:]
    neighbours = indices[:,:K+1]
    # replace end

    h = np.ones((ndata,ndata),dtype=float)*np.inf
    for i in range(ndata):
        h[i,neighbours[i,:]] = d[i,neighbours[i,:]]
    h = Dijkstra(h)
    return MDS(h,newdim)

Answer 1

Try sklearn.neighbors.radius_neighbors_graph for your distance matrix