From 572fb38ce4fddc5c776481b24b72c012311234ed Mon Sep 17 00:00:00 2001 From: Robert Blackwell Date: Tue, 11 Apr 2023 12:07:26 -0400 Subject: [PATCH] Merge embedding function files #8 --- ManifoldEM/DMembeddingII.py | 103 +++++++++++++++++++++++++++++++++- ManifoldEM/sembeddingonFly.py | 37 ------------ ManifoldEM/slaplacianonFly.py | 80 -------------------------- 3 files changed, 100 insertions(+), 120 deletions(-) delete mode 100644 ManifoldEM/sembeddingonFly.py delete mode 100644 ManifoldEM/slaplacianonFly.py diff --git a/ManifoldEM/DMembeddingII.py b/ManifoldEM/DMembeddingII.py index 2eb6d39..458e11c 100644 --- a/ManifoldEM/DMembeddingII.py +++ b/ManifoldEM/DMembeddingII.py @@ -9,10 +9,107 @@ Copyright (c) UWM, Ali Dashti 2016 (original matlab version) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Copyright (c) Columbia University Hstau Liao 2019 (python version) -Copyright (c) Columbia University Evan Seitz 2019 (python version) - +Copyright (c) Columbia University Evan Seitz 2019 (python version) + ''' +def sembedding(yVal, yCol, yRow, nS, options1): + """ + Laplacian eigenfunction embedding using sparse arrays + """ + # Copyright (c) UWM, Ali Dashti 2016 (original matlab version) + # Copyright (c) Columbia University Hstau Liao 2018 (python version) + + options = namedtuple('options', 'sigma alpha visual nEigs autotune') + options.sigma = options1.sigma + options.alpha = options1.alpha + options.nEigs = options1.nEigs + options.autotune = 0 + + l, sigmaTune = slaplacian(yVal, yCol, yRow, nS, options) + try: + vals, vecs = eigsh(l, k=options.nEigs + 1, maxiter=300) + except ArpackNoConvergence as e: + + vals = e.eigenvalues + vecs = e.eigenvectors + print("eigsh not converging in 300 iterations...") + ix = np.argsort(vals)[::-1] + vals = np.sort(vals)[::-1] + vecs = vecs[:, ix] + + return (vals, vecs) + + +def slaplacian(*arg): + """ + Given a set of nS data points, and the dinstances to nN nearest neighbors + for each data point, slaplacian computes a sparse, nY by nY symmetric + graph Laplacian matrix l. + + The input data are supplied in the column vectors yVal and yInd of length + nY * nN such that + + yVal( ( i - 1 ) * nN + ( 1 : nN ) ) contains the distances to the + nN nearest neighbors of data point i sorted in ascending order, and + + yInd( ( i - 1 ) * nN + ( 1 : nN ) ) contains the indices of the nearest + neighbors. + + yVal and yInd can be computed by calling nndist + + slaplacian admits a number of options passed as name-value pairs + + alpha : normalization, according to Coifman & Lafon + + nAutotune : number of nearest neighbors for autotuning. Set to zero if no + autotuning is to be performed + + sigma: width of the Gaussian kernel + + Copyright (c) UWM, Ali Dashti 2016 (original matlab version) + Copyright (c) Columbia University Hstau Liao 2019 (python version) + Copyright (c) Columbia University Evan Seitz 2019 (python version) + """ + yVal = arg[0] + yCol = arg[1] + yRow = arg[2] + nS = arg[3] #dataset size + options = arg[4] #options.sigma: Gaussian width + + nNZ = len(yVal) #number of nonzero elements + + # if required, compute autotuning distances: + if options.autotune > 0: + print('Autotuning is not implemented in this version of slaplacian' + '\n') + else: + sigmaTune = options.sigma + + + yVal = yVal / sigmaTune**2 + + # compute the unnormalized weight matrix: + yVal = np.exp(-yVal) #apply exponential weights (yVal is distance**2) + l = csc_matrix((yVal, (yRow, yCol)), shape=(nS, nS)) + d = np.array(l.sum(axis=0)).T + + if options.alpha != 1: #apply non-isotropic normalization + d = d**options.alpha + + yVal = yVal / (d[yRow].flatten('C') * d[yCol].flatten('C')) + l = csc_matrix((yVal, (yRow, yCol)), shape=(nS, nS)) + + # normalize by the degree matrix to form normalized graph Laplacian: + d = np.array(l.sum(axis=0)) + d = np.sqrt(d).T + + yVal = yVal / (d[yRow].flatten('C') * d[yCol].flatten('C')) + l = csc_matrix((yVal, (yRow, yCol)), shape=(nS, nS)) + l = np.abs(l + l.T) / 2.0 #iron out numerical wrinkles + temp = l - l.T + + return (l, sigmaTune) + def get_yColVal(params): @@ -159,7 +256,7 @@ def op(D, k, tune, prefsigma): #*arg options.visual = visual options.nEigs = nEigs - lamb, v = sembeddingonFly.op(yVal, yCol, yRow, nS, options) + lamb, v = sembedding(yVal, yCol, yRow, nS, options) #psi = v[:, 1 : nEigs+1]/np.tile(v[:, 0 ].reshape((-1,1)), (1, nEigs)) true_shape = v.shape[1] - 1 diff --git a/ManifoldEM/sembeddingonFly.py b/ManifoldEM/sembeddingonFly.py deleted file mode 100644 index acf2746..0000000 --- a/ManifoldEM/sembeddingonFly.py +++ /dev/null @@ -1,37 +0,0 @@ -import numpy as np - -from collections import namedtuple -from scipy.sparse.linalg import eigsh, ArpackNoConvergence - -from ManifoldEM import slaplacianonFly -""" -%SEMBEDDING Laplacian eigenfunction embedding using sparse arrays - -""" -''' -Copyright (c) UWM, Ali Dashti 2016 (original matlab version) -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Copyright (c) Columbia University Hstau Liao 2018 (python version) -''' - - -def op(yVal, yCol, yRow, nS, options1): - options = namedtuple('options', 'sigma alpha visual nEigs autotune') - options.sigma = options1.sigma - options.alpha = options1.alpha - options.nEigs = options1.nEigs - options.autotune = 0 - - l, sigmaTune = slaplacianonFly.op(yVal, yCol, yRow, nS, options) - try: - vals, vecs = eigsh(l, k=options.nEigs + 1, maxiter=300) - except ArpackNoConvergence as e: - - vals = e.eigenvalues - vecs = e.eigenvectors - print("eigsh not converging in 300 iterations...") - ix = np.argsort(vals)[::-1] - vals = np.sort(vals)[::-1] - vecs = vecs[:, ix] - - return (vals, vecs) diff --git a/ManifoldEM/slaplacianonFly.py b/ManifoldEM/slaplacianonFly.py deleted file mode 100644 index aea9583..0000000 --- a/ManifoldEM/slaplacianonFly.py +++ /dev/null @@ -1,80 +0,0 @@ -import numpy as np -import time -from scipy.sparse import csc_matrix -from collections import namedtuple -""" -function [ l, sigmaTune ] = slaplacianonFly( yVal,yCol,yRow, nS, varargin ) -%SLAPLACIAN Sparse Laplacian matrix -% -% Given a set of nS data points, and the dinstances to nN nearest neighbors -% for each data point, slaplacian computes a sparse, nY by nY symmetric -% graph Laplacian matrix l. -% -% The input data are supplied in the column vectors yVal and yInd of length -% nY * nN such that -% -% yVal( ( i - 1 ) * nN + ( 1 : nN ) ) contains the distances to the -% nN nearest neighbors of data point i sorted in ascending order, and -% -% yInd( ( i - 1 ) * nN + ( 1 : nN ) ) contains the indices of the nearest -% neighbors. -% -% yVal and yInd can be computed by calling nndist -% -% slaplacian admits a number of options passed as name-value pairs -% -% alpha : normalization, according to Coifman & Lafon -% -% nAutotune : number of nearest neighbors for autotuning. Set to zero if no -% autotuning is to be performed -% -% sigma: width of the Gaussian kernel -% -% -Copyright (c) UWM, Ali Dashti 2016 (original matlab version) -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Copyright (c) Columbia University Hstau Liao 2019 (python version) -Copyright (c) Columbia University Evan Seitz 2019 (python version) - -""" - - -def op(*arg): - yVal = arg[0] - yCol = arg[1] - yRow = arg[2] - nS = arg[3] #dataset size - options = arg[4] #options.sigma: Gaussian width - - nNZ = len(yVal) #number of nonzero elements - - # if required, compute autotuning distances: - if options.autotune > 0: - print('Autotuning is not implemented in this version of slaplacian' + '\n') - else: - sigmaTune = options.sigma - - - yVal = yVal / sigmaTune**2 - - # compute the unnormalized weight matrix: - yVal = np.exp(-yVal) #apply exponential weights (yVal is distance**2) - l = csc_matrix((yVal, (yRow, yCol)), shape=(nS, nS)) - d = np.array(l.sum(axis=0)).T - - if options.alpha != 1: #apply non-isotropic normalization - d = d**options.alpha - - yVal = yVal / (d[yRow].flatten('C') * d[yCol].flatten('C')) - l = csc_matrix((yVal, (yRow, yCol)), shape=(nS, nS)) - - # normalize by the degree matrix to form normalized graph Laplacian: - d = np.array(l.sum(axis=0)) - d = np.sqrt(d).T - - yVal = yVal / (d[yRow].flatten('C') * d[yCol].flatten('C')) - l = csc_matrix((yVal, (yRow, yCol)), shape=(nS, nS)) - l = np.abs(l + l.T) / 2.0 #iron out numerical wrinkles - temp = l - l.T - - return (l, sigmaTune)