# 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
    nbrs = list()
    for ii in range(len(sample)):
        w = np.array(A[sample[ii],:].todense()).flatten()
        p = w / np.sum(w)
        nbrs.append(np.where(p > 0)[0])
        Ps.append(p[nbrs[ii]])
    return Ps, nbrs

def project_tangent(X, p, nbr, d):
    Z = (X[nbr, :] - np.dot(p, X[nbr, :]))*p[:, np.newaxis]
    sig = np.dot(Z.transpose(), Z)
    e_vals, e_vecs = np.linalg.eigh(sig)
    j = e_vals.argsort()[::-1]  # Returns indices that will sort array from greatest to least.
    e_vec = e_vecs[:, j]
    e_vec = e_vec[:, :d]  # Gets d largest eigenvectors.
    X_tangent = np.dot(X[nbr, :], e_vec)  # Project local points onto tangent plane
    return X_tangent

def distortion(X, L, sample, PS, nbrs, n, d):
    n = len(sample)
    dist = 0.0
    nsum = n
    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

# 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

    nbrs = list()
    for ii in range(len(sample)):
        w = np.array(A[sample[ii],:].todense()).flatten()
        p = w / np.sum(w)
        nbrs.append(np.where(p > 0)[0])
        Ps.append(p[nbrs[ii]])
    return Ps, nbrs

def project_tangent(X, p, nbr, d):
    Z = (X[nbr, :] - np.dot(p, X[nbr, :]))*p[:, np.newaxis]
    sig = np.dot(Z.transpose(), Z)
    e_vals, e_vecs = np.linalg.eigh(sig)
    j = e_vals.argsort()[::-1]  # Returns indices that will sort array from greatest to least.
    e_vec = e_vecs[:, j]
    e_vec = e_vec[:, :d]  # Gets d largest eigenvectors.
    X_tangent = np.dot(X[nbr, :], e_vec)  # Project local points onto tangent plane
    return X_tangent

def distortion(X, L, sample, PS, nbrs, n, d):
    n = len(sample)
    dist = 0.0
    nsum = n
    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]