# Path: megaman/geometry/rmetric.py
# def riemann_metric_lazy( Y, sample,laplacian, n_dim, invert_h=False, mode_inv = 'svd'):
#     """
#     Parameters
#     ----------
#     Y: array-like, shape = (n_samples, mdimY )
#         The embedding coordinates of the points
#     sample : array-like, shape (n_samples)
#         array of indices over which to calculate the riemannian metric.
#     laplacian: array-like, shape = (n_samples, n_samples)
#         The Laplacian of the data. It is recommended to use the "geometric"
#         Laplacian (default) option from geometry.graph_laplacian()
#     n_dim : integer, optional
#         Use only the first n_dim <= mdimY dimensions.All dimensions
#         n_dim:mdimY are ignored.
#     invert_h: boolean, optional
#         if False, only the "dual Riemannian metric" is computed
#         if True, the dual metric matrices are inverted to obtain the
#         Riemannian metric G.
#     mode_inv: string, optional
#        How to compute the inverses of h_dual_metric, if invert_h
#         "inv", use numpy.inv()
#         "svd" (default), use numpy.linalg.svd(), then invert the eigenvalues
#         (possibly a more numerically stable method with H is symmetric and
#         ill conditioned)
# 
#     Returns
#     -------
#     h_dual_metric : array, shape=(n_samples, n_dim, n_dim)
#     Optionally :
#     g_riemann_metric : array, shape=(n_samples, n_dim, n_dim )
#     Hvv : singular vectors of H, transposed, shape = ( n_samples, n_dim, n_dim )
#     Hsvals : singular values of H, shape = ( n_samples, n_dim )
#     Gsvals : singular values of G, shape = ( n_samples, n_dim )
# 
#     Notes
#     -----
#     References
#     ----------
#     "Non-linear dimensionality reduction: Riemannian metric estimation and
#     the problem of geometric discovery",
#     Dominique Perraul-Joncas, Marina Meila, arXiv:1305.7255
#     """
#     n_samples = laplacian.shape[0]
#     laplacian = laplacian[sample,:]
#     n_dim_Y = Y.shape[1]
#     h_dual_metric_full = np.zeros((len(sample), n_dim_Y, n_dim_Y))
#     h_dual_metric = np.zeros((len(sample), n_dim, n_dim))
#     for i in np.arange(n_dim_Y):
#         for j in np.arange(i, n_dim_Y):
#             yij = Y[:,i]*Y[:,j]
#             h_dual_metric_full[ :, i, j] = 0.5*(laplacian.dot(yij)-Y[sample,j]*laplacian.dot(Y[:,i])-Y[sample,i]*laplacian.dot(Y[:,j]))
#     for j in np.arange(n_dim_Y - 1):
#         for i in np.arange(j+1, n_dim_Y):
#             h_dual_metric_full[ :,i,j] = h_dual_metric_full[:,j,i]
# 
#     riemann_metric, h_dual_metric, Hvv, Hsvals, Gsvals = compute_G_from_H(h_dual_metric_full)
#     return h_dual_metric,riemann_metric, Hvv, Hsvals, Gsvals
# 
# Path: megaman/geometry/affinity.py
# def compute_affinity_matrix(adjacency_matrix, method='auto', **kwargs):
#     """Compute the affinity matrix with the given method"""
#     if method == 'auto':
#         method = 'gaussian'
#     return Affinity.init(method, **kwargs).affinity_matrix(adjacency_matrix)
# 
# Path: megaman/geometry/laplacian.py
# def compute_laplacian_matrix(affinity_matrix, method='auto', **kwargs):
#     """Compute the laplacian matrix with the given method"""
#     if method == 'auto':
#         method = 'geometric'
#     return Laplacian.init(method, **kwargs).laplacian_matrix(affinity_matrix)


import numpy as np
import sys, time
import multiprocessing as mp
from megaman.geometry.rmetric import riemann_metric_lazy
from megaman.geometry.affinity import compute_affinity_matrix
from megaman.geometry.laplacian import compute_laplacian_matrix
from multiprocessing.pool import ThreadPool
    for i in range(n):
        p = PS[i]
        nbr = nbrs[i]
        if len(nbr) > 1:
            X_t = project_tangent(X, p, nbr, d)
            X_tangent = X.copy()
            X_tangent = X_tangent[:, :d]
            X_tangent[nbr, :] = X_t # rmetric only depends on nbrs
            H = riemann_metric_lazy(X_tangent, sample, L, d)[0]
            dist += np.linalg.norm(H[i, :, :] - np.eye(d), ord=2)
        else:
            nsum -= 1
    if nsum > 0:
        distortion = dist / nsum
    else:
        distortion = float('Inf')
    return distortion

def evaluate_radius(radius, d, sample, rad_bw_ratio=3.0):
    global DIST
    global X

    t0 = time.time()
    D = DIST.copy()
    (n, dim) = X.shape
    D.data[D.data > radius] = 0.0
    D.eliminate_zeros()
    h = radius / rad_bw_ratio
    A = compute_affinity_matrix(D, 'gaussian', radius=h)
    (PS, nbrs) = compute_nbr_wts(A, sample)

# Path: megaman/geometry/rmetric.py
# def riemann_metric_lazy( Y, sample,laplacian, n_dim, invert_h=False, mode_inv = 'svd'):
#     """
#     Parameters
#     ----------
#     Y: array-like, shape = (n_samples, mdimY )
#         The embedding coordinates of the points
#     sample : array-like, shape (n_samples)
#         array of indices over which to calculate the riemannian metric.
#     laplacian: array-like, shape = (n_samples, n_samples)
#         The Laplacian of the data. It is recommended to use the "geometric"
#         Laplacian (default) option from geometry.graph_laplacian()
#     n_dim : integer, optional
#         Use only the first n_dim <= mdimY dimensions.All dimensions
#         n_dim:mdimY are ignored.
#     invert_h: boolean, optional
#         if False, only the "dual Riemannian metric" is computed
#         if True, the dual metric matrices are inverted to obtain the
#         Riemannian metric G.
#     mode_inv: string, optional
#        How to compute the inverses of h_dual_metric, if invert_h
#         "inv", use numpy.inv()
#         "svd" (default), use numpy.linalg.svd(), then invert the eigenvalues
#         (possibly a more numerically stable method with H is symmetric and
#         ill conditioned)
# 
#     Returns
#     -------
#     h_dual_metric : array, shape=(n_samples, n_dim, n_dim)
#     Optionally :
#     g_riemann_metric : array, shape=(n_samples, n_dim, n_dim )
#     Hvv : singular vectors of H, transposed, shape = ( n_samples, n_dim, n_dim )
#     Hsvals : singular values of H, shape = ( n_samples, n_dim )
#     Gsvals : singular values of G, shape = ( n_samples, n_dim )
# 
#     Notes
#     -----
#     References
#     ----------
#     "Non-linear dimensionality reduction: Riemannian metric estimation and
#     the problem of geometric discovery",
#     Dominique Perraul-Joncas, Marina Meila, arXiv:1305.7255
#     """
#     n_samples = laplacian.shape[0]
#     laplacian = laplacian[sample,:]
#     n_dim_Y = Y.shape[1]
#     h_dual_metric_full = np.zeros((len(sample), n_dim_Y, n_dim_Y))
#     h_dual_metric = np.zeros((len(sample), n_dim, n_dim))
#     for i in np.arange(n_dim_Y):
#         for j in np.arange(i, n_dim_Y):
#             yij = Y[:,i]*Y[:,j]
#             h_dual_metric_full[ :, i, j] = 0.5*(laplacian.dot(yij)-Y[sample,j]*laplacian.dot(Y[:,i])-Y[sample,i]*laplacian.dot(Y[:,j]))
#     for j in np.arange(n_dim_Y - 1):
#         for i in np.arange(j+1, n_dim_Y):
#             h_dual_metric_full[ :,i,j] = h_dual_metric_full[:,j,i]
# 
#     riemann_metric, h_dual_metric, Hvv, Hsvals, Gsvals = compute_G_from_H(h_dual_metric_full)
#     return h_dual_metric,riemann_metric, Hvv, Hsvals, Gsvals
# 
# Path: megaman/geometry/affinity.py
# def compute_affinity_matrix(adjacency_matrix, method='auto', **kwargs):
#     """Compute the affinity matrix with the given method"""
#     if method == 'auto':
#         method = 'gaussian'
#     return Affinity.init(method, **kwargs).affinity_matrix(adjacency_matrix)
# 
# Path: megaman/geometry/laplacian.py
# def compute_laplacian_matrix(affinity_matrix, method='auto', **kwargs):
#     """Compute the laplacian matrix with the given method"""
#     if method == 'auto':
#         method = 'geometric'
#     return Laplacian.init(method, **kwargs).laplacian_matrix(affinity_matrix)


# Path: megaman/utils/estimate_radius.py
import numpy as np
import sys, time
import multiprocessing as mp
from megaman.geometry.rmetric import riemann_metric_lazy
from megaman.geometry.affinity import compute_affinity_matrix
from megaman.geometry.laplacian import compute_laplacian_matrix
from multiprocessing.pool import ThreadPool

    for i in range(n):
        p = PS[i]
        nbr = nbrs[i]
        if len(nbr) > 1:
            X_t = project_tangent(X, p, nbr, d)
            X_tangent = X.copy()
            X_tangent = X_tangent[:, :d]
            X_tangent[nbr, :] = X_t # rmetric only depends on nbrs
            H = riemann_metric_lazy(X_tangent, sample, L, d)[0]
            dist += np.linalg.norm(H[i, :, :] - np.eye(d), ord=2)
        else:
            nsum -= 1
    if nsum > 0:
        distortion = dist / nsum
    else:
        distortion = float('Inf')
    return distortion

def evaluate_radius(radius, d, sample, rad_bw_ratio=3.0):
    global DIST
    global X

    t0 = time.time()
    D = DIST.copy()
    (n, dim) = X.shape
    D.data[D.data > radius] = 0.0
    D.eliminate_zeros()
    h = radius / rad_bw_ratio
    A = compute_affinity_matrix(D, 'gaussian', radius=h)
    (PS, nbrs) = compute_nbr_wts(A, sample)

    L = compute_laplacian_matrix(A, method='geometric', scaling_epps=h)