Module pyfocusr.eigsort
Expand source code
import numpy as np
from scipy.stats import wasserstein_distance
from scipy.optimize import linear_sum_assignment
from scipy.spatial import KDTree
from .main import *
class eigsort(object):
"""
Class to sort eigen vectors/ spectral coordinates necessary before aligning/registering spectral coordinates.
"""
def __init__(self,
graph_target,
graph_source,
n_features,
target_as_reference=True, # if `True``, then target order used as reference and source permuted
# to match target. If `False`, then source order used as reference and
# target is permuted.
):
# Input variables
self.graph_target = graph_target # Target mesh (match source eigs to this)
self.graph_source = graph_source # Source mesh (eigs to re-order/flip to match target).
self.n_features = n_features # number of features (eigenvecs/values) to consider during alignment.
self.target_as_reference = target_as_reference # If true, then the target mesh is used as the reference
# Build/create specified data needed for alignment.
# Points
self.rand_target_points = self.graph_target.get_rand_normalized_points() # Get rand target points for alignment
self.rand_source_points = self.graph_source.get_rand_normalized_points() # Get rand source points for alignment
# Eigenvectors
self.rand_target_eig_vecs = self.graph_target.get_rand_eig_vecs() # Get eigvecs of rand target points
self.rand_source_eig_vecs = self.graph_source.get_rand_eig_vecs() # Get eigvecs of rand source points
# Interim values used in class.
self.c_lambda = np.zeros((self.n_features, self.n_features)) # eigenvalue cost matrix
self.c_hist = np.zeros_like(self.c_lambda) # histogram comparison matrix
self.c_hist_f = np.zeros_like(self.c_lambda) # flipped histogram cost matrix.
self.c_spatial = np.zeros_like(self.c_lambda) # spatial cost matrix
self.c_spatial_f = np.zeros_like(self.c_lambda) # flipped spatial cost matrix
# Returns.
self.Q = None # Cost of making each particular "pairing" of eigenvectors/values between two graphs.
def eigen_sort(self):
"""
- Calculate a metric of similarity c for straight matches, and matches from flipped eigenvectors c_f for source
mesh.
- Get the dissimilarity matrix Q by taking the min of c and c_f
- Get matrix S of positions which were flipped to get minimum
- Get list of `row_flipped` and `cols_flipped`
- Get list of matching columns for each row by using hungarian algorithm (linear_sum_assignment) on Q
- From the matches, and the identified flipped rows/cols, identify which of the pairs are flipped
- Flip the eigen vector values for the necessary eigen vectors on the source mesh
"""
# find best matching (negative/positive)
c = self.c_spatial * self.c_lambda * self.c_hist
c_f = self.c_spatial_f * self.c_lambda * self.c_hist_f
self.Q = np.min((c, c_f), axis=0)
# telling us if sign flipped or not. True = same sign, False = flipped
S = c > c_f
(target_flipped, source_flipped) = np.where(S == True)
# Get matches
if self.target_as_reference is True:
target_matches, source_matches = linear_sum_assignment(self.Q)
elif self.target_as_reference is False:
source_matches, target_matches = linear_sum_assignment(self.Q.T)
# The original Matlab code calculates Q as:
# Q = sum(M'*Q')'; which is equivalent to np.sum(np.matmul(M.T, Q.T), axis=0)
# Where M = nXn matrix where n is number of spectral features (or the length of target_matches). and
# has 1s where there is a "match", i.e.
# M = np.zeros((len(target_matches), len(target_matches)))
# M[target_matches, source_matches] = 1
# This results in a permutation cost equal to all of the "costs" for the entire row. This seems wrong,
# it seems like the cost of using a particular "pair" should come down to the cost of just that particular
# match and not all of the matches that it "didnt" use. In fact, bad matches would have high numbers and screw
# things up. We are therefore doing the below.
self.Q = self.Q[target_matches, source_matches]
# This will identify the source & target eigen #s (pairs after matching) that need to be flipped.
# So, if we're on eigen 3 for target, but its pair is eigen 5, then this will identify the pair (3, 5).
# The for loop then iterates through and turns the 5 (from the source) to be negative and leave it where
# it is. Finally, we will then select the final eigenvecs in the appropriate order (after flipping)
# using their indices from target_matches, source_matches.
# If a "flipped pair" is also a "best match", then store the pairs in flipped_pairs so that the corresponding
# source eig_vecs can be flipped to match the target eig_vec orientation.
flipped_pairs = [p2 for p1 in zip(target_flipped, source_flipped) for p2 in zip(target_matches, source_matches) if
p2 == p1]
# Flipped pairs are used to flip the sign of the source eig_vecs.
for mode_0, mode_1 in flipped_pairs:
if self.target_as_reference is True:
self.graph_source.eig_vecs[:, mode_1] = self.graph_source.eig_vecs[:, mode_1] * -1
elif self.target_as_reference is False:
self.graph_target.eig_vecs[:, mode_0] = self.graph_target.eig_vecs[:, mode_0] * -1
# The source eig_vecs are re-ordered to match the order from the target eig_vecs - based on the best matches
# Identified using Hungarian algorithm (linear_sum_assignment) on the dissimilarity matrix Q.
if self.target_as_reference is True:
self.graph_source.eig_vecs[:, target_matches] = self.graph_source.eig_vecs[:, source_matches]
elif self.target_as_reference is False:
self.graph_target.eig_vecs[:, source_matches] = self.graph_target.eig_vecs[:, target_matches]
print_header('Eigenvector Sorting Results')
if self.target_as_reference is True:
print('Using target eigenmaps as the reference')
elif self.target_as_reference is False:
print('Using source eigenmaps as the reference')
print('The matches for eigenvectors were as follows:')
print('Target\t| Source')
for matched_pair in zip(target_matches, source_matches):
source_value = str(matched_pair[1])
target_value = str(matched_pair[0])
if matched_pair in flipped_pairs:
if self.target_as_reference is True:
source_value = '-' + source_value
elif self.target_as_reference is False:
target_value = '-' + target_value
print('{:6}\t| {:6}'.format(target_value, source_value))
print('*Negative source values means those eigenvectors were flipped*\n ')
def calc_c_lambda(self):
"""
Get average gap between successive eigenvalues (used to normalize any errors/differences between
the eigenvalues of the two meshes.
:return:
"""
for graph in [self.graph_source, self.graph_target]:
if graph.eig_val_gap is None:
graph.get_eig_val_gap()
# average the gap of mesh 1 and mesh 2.
eigen_gap = (self.graph_target.eig_val_gap + self.graph_source.eig_val_gap)/2
# Calculate difference in eigenvalues between meshes normalizing to eigen_gap.
for i in range(self.n_features):
for j in range(self.n_features):
self.c_lambda[i, j] = np.exp((self.graph_target.eig_vals[i] - self.graph_source.eig_vals[j]) ** 2 /
(2 * eigen_gap ** 2)
)
def calc_c_hist(self):
"""
Compare the histograms of all of the eigenvectors (eigenfunctions) for each of the
meshes - only upto the number of meshes we are interested in. The paper says to do
this for the number of features we are interested (i.e. 5 or 6), however, this implementation
does that for more than the requested number of features to ensure that we dont get the wrong ones.
So, we'll get more spectral coordinates, order them, and then select the appropriate number from
these re-ordered coordinates.
:return:
"""
# Initially tried using straight values (not log) but eig vec 1 got accentuated too much (in the weighting)
# Then tried just log, but becuase there are negative values it creates erro.
# need to add .5 + a small value to ensure there are no 0 values entered to log.
# wasserstein_distance is the same as earth movers distance, and is the minimum "work" needed to
# transform u (first entry) into v (second entry).
eps = np.finfo(float).eps
for i in range(self.n_features):
for j in range(self.n_features):
self.c_hist[i, j] = wasserstein_distance(np.log(self.rand_target_eig_vecs[:, i]
+ 0.5
+ eps),
np.log(self.rand_source_eig_vecs[:, j]
+ 0.5
+ eps))
self.c_hist_f[i, j] = wasserstein_distance(np.log(self.rand_target_eig_vecs[:, i]
+ 0.5
+ eps),
np.log(-self.rand_source_eig_vecs[:, j]
+ 0.5
+ eps))
def calc_c_spatial(self):
"""
Next, we calculate the spectral distance of points that are nearest to one another on the mesh.
This must assume that the meshes are crudely aligned. This will mean that for knees, we will
always have to flip any left knees to match right knees, we should then use an ICP algoritm to
get a crude starting alignment. However, even without the ICP (as long as left are flipped)
this algorithm scales the sizes of the bounding boxes for the points, so they should be very
crudely aligned.
:return:
"""
source_kd_tree = KDTree(self.rand_source_points)
_, idx_source_for_each_target_pt = source_kd_tree.query(self.rand_target_points)
# Using the same error as the Matlab code., which is the average euclidean error in n-dimensions.
# Paper says that it used the sum of squares.
for i in range(self.n_features):
for j in range(self.n_features):
self.c_spatial[i, j] = np.sqrt(np.sum((self.rand_source_eig_vecs[idx_source_for_each_target_pt, j]
- self.rand_target_eig_vecs[:, i])**2)
)/self.rand_target_eig_vecs.shape[0]
self.c_spatial_f[i, j] = np.sqrt(np.sum((-self.rand_source_eig_vecs[idx_source_for_each_target_pt, j]
- self.rand_target_eig_vecs[:, i])**2)
)/self.rand_target_eig_vecs.shape[0]
def sort_eigenmaps(self):
"""
Run functions necessary to sort eigenvalues.
Seems to do a good job identifying flips.
However, if there are abnormal eigenmaps it doesnt do very well.
:return:
"""
self.calc_c_lambda() # lambda = eigenvalues.
self.calc_c_hist()
self.calc_c_spatial()
self.eigen_sort()
return self.Q # Q was calculated in eigen_sort and is the cost associated with each of the final pairings.
Classes
class eigsort (graph_target, graph_source, n_features, target_as_reference=True)
-
Class to sort eigen vectors/ spectral coordinates necessary before aligning/registering spectral coordinates.
Expand source code
class eigsort(object): """ Class to sort eigen vectors/ spectral coordinates necessary before aligning/registering spectral coordinates. """ def __init__(self, graph_target, graph_source, n_features, target_as_reference=True, # if `True``, then target order used as reference and source permuted # to match target. If `False`, then source order used as reference and # target is permuted. ): # Input variables self.graph_target = graph_target # Target mesh (match source eigs to this) self.graph_source = graph_source # Source mesh (eigs to re-order/flip to match target). self.n_features = n_features # number of features (eigenvecs/values) to consider during alignment. self.target_as_reference = target_as_reference # If true, then the target mesh is used as the reference # Build/create specified data needed for alignment. # Points self.rand_target_points = self.graph_target.get_rand_normalized_points() # Get rand target points for alignment self.rand_source_points = self.graph_source.get_rand_normalized_points() # Get rand source points for alignment # Eigenvectors self.rand_target_eig_vecs = self.graph_target.get_rand_eig_vecs() # Get eigvecs of rand target points self.rand_source_eig_vecs = self.graph_source.get_rand_eig_vecs() # Get eigvecs of rand source points # Interim values used in class. self.c_lambda = np.zeros((self.n_features, self.n_features)) # eigenvalue cost matrix self.c_hist = np.zeros_like(self.c_lambda) # histogram comparison matrix self.c_hist_f = np.zeros_like(self.c_lambda) # flipped histogram cost matrix. self.c_spatial = np.zeros_like(self.c_lambda) # spatial cost matrix self.c_spatial_f = np.zeros_like(self.c_lambda) # flipped spatial cost matrix # Returns. self.Q = None # Cost of making each particular "pairing" of eigenvectors/values between two graphs. def eigen_sort(self): """ - Calculate a metric of similarity c for straight matches, and matches from flipped eigenvectors c_f for source mesh. - Get the dissimilarity matrix Q by taking the min of c and c_f - Get matrix S of positions which were flipped to get minimum - Get list of `row_flipped` and `cols_flipped` - Get list of matching columns for each row by using hungarian algorithm (linear_sum_assignment) on Q - From the matches, and the identified flipped rows/cols, identify which of the pairs are flipped - Flip the eigen vector values for the necessary eigen vectors on the source mesh """ # find best matching (negative/positive) c = self.c_spatial * self.c_lambda * self.c_hist c_f = self.c_spatial_f * self.c_lambda * self.c_hist_f self.Q = np.min((c, c_f), axis=0) # telling us if sign flipped or not. True = same sign, False = flipped S = c > c_f (target_flipped, source_flipped) = np.where(S == True) # Get matches if self.target_as_reference is True: target_matches, source_matches = linear_sum_assignment(self.Q) elif self.target_as_reference is False: source_matches, target_matches = linear_sum_assignment(self.Q.T) # The original Matlab code calculates Q as: # Q = sum(M'*Q')'; which is equivalent to np.sum(np.matmul(M.T, Q.T), axis=0) # Where M = nXn matrix where n is number of spectral features (or the length of target_matches). and # has 1s where there is a "match", i.e. # M = np.zeros((len(target_matches), len(target_matches))) # M[target_matches, source_matches] = 1 # This results in a permutation cost equal to all of the "costs" for the entire row. This seems wrong, # it seems like the cost of using a particular "pair" should come down to the cost of just that particular # match and not all of the matches that it "didnt" use. In fact, bad matches would have high numbers and screw # things up. We are therefore doing the below. self.Q = self.Q[target_matches, source_matches] # This will identify the source & target eigen #s (pairs after matching) that need to be flipped. # So, if we're on eigen 3 for target, but its pair is eigen 5, then this will identify the pair (3, 5). # The for loop then iterates through and turns the 5 (from the source) to be negative and leave it where # it is. Finally, we will then select the final eigenvecs in the appropriate order (after flipping) # using their indices from target_matches, source_matches. # If a "flipped pair" is also a "best match", then store the pairs in flipped_pairs so that the corresponding # source eig_vecs can be flipped to match the target eig_vec orientation. flipped_pairs = [p2 for p1 in zip(target_flipped, source_flipped) for p2 in zip(target_matches, source_matches) if p2 == p1] # Flipped pairs are used to flip the sign of the source eig_vecs. for mode_0, mode_1 in flipped_pairs: if self.target_as_reference is True: self.graph_source.eig_vecs[:, mode_1] = self.graph_source.eig_vecs[:, mode_1] * -1 elif self.target_as_reference is False: self.graph_target.eig_vecs[:, mode_0] = self.graph_target.eig_vecs[:, mode_0] * -1 # The source eig_vecs are re-ordered to match the order from the target eig_vecs - based on the best matches # Identified using Hungarian algorithm (linear_sum_assignment) on the dissimilarity matrix Q. if self.target_as_reference is True: self.graph_source.eig_vecs[:, target_matches] = self.graph_source.eig_vecs[:, source_matches] elif self.target_as_reference is False: self.graph_target.eig_vecs[:, source_matches] = self.graph_target.eig_vecs[:, target_matches] print_header('Eigenvector Sorting Results') if self.target_as_reference is True: print('Using target eigenmaps as the reference') elif self.target_as_reference is False: print('Using source eigenmaps as the reference') print('The matches for eigenvectors were as follows:') print('Target\t| Source') for matched_pair in zip(target_matches, source_matches): source_value = str(matched_pair[1]) target_value = str(matched_pair[0]) if matched_pair in flipped_pairs: if self.target_as_reference is True: source_value = '-' + source_value elif self.target_as_reference is False: target_value = '-' + target_value print('{:6}\t| {:6}'.format(target_value, source_value)) print('*Negative source values means those eigenvectors were flipped*\n ') def calc_c_lambda(self): """ Get average gap between successive eigenvalues (used to normalize any errors/differences between the eigenvalues of the two meshes. :return: """ for graph in [self.graph_source, self.graph_target]: if graph.eig_val_gap is None: graph.get_eig_val_gap() # average the gap of mesh 1 and mesh 2. eigen_gap = (self.graph_target.eig_val_gap + self.graph_source.eig_val_gap)/2 # Calculate difference in eigenvalues between meshes normalizing to eigen_gap. for i in range(self.n_features): for j in range(self.n_features): self.c_lambda[i, j] = np.exp((self.graph_target.eig_vals[i] - self.graph_source.eig_vals[j]) ** 2 / (2 * eigen_gap ** 2) ) def calc_c_hist(self): """ Compare the histograms of all of the eigenvectors (eigenfunctions) for each of the meshes - only upto the number of meshes we are interested in. The paper says to do this for the number of features we are interested (i.e. 5 or 6), however, this implementation does that for more than the requested number of features to ensure that we dont get the wrong ones. So, we'll get more spectral coordinates, order them, and then select the appropriate number from these re-ordered coordinates. :return: """ # Initially tried using straight values (not log) but eig vec 1 got accentuated too much (in the weighting) # Then tried just log, but becuase there are negative values it creates erro. # need to add .5 + a small value to ensure there are no 0 values entered to log. # wasserstein_distance is the same as earth movers distance, and is the minimum "work" needed to # transform u (first entry) into v (second entry). eps = np.finfo(float).eps for i in range(self.n_features): for j in range(self.n_features): self.c_hist[i, j] = wasserstein_distance(np.log(self.rand_target_eig_vecs[:, i] + 0.5 + eps), np.log(self.rand_source_eig_vecs[:, j] + 0.5 + eps)) self.c_hist_f[i, j] = wasserstein_distance(np.log(self.rand_target_eig_vecs[:, i] + 0.5 + eps), np.log(-self.rand_source_eig_vecs[:, j] + 0.5 + eps)) def calc_c_spatial(self): """ Next, we calculate the spectral distance of points that are nearest to one another on the mesh. This must assume that the meshes are crudely aligned. This will mean that for knees, we will always have to flip any left knees to match right knees, we should then use an ICP algoritm to get a crude starting alignment. However, even without the ICP (as long as left are flipped) this algorithm scales the sizes of the bounding boxes for the points, so they should be very crudely aligned. :return: """ source_kd_tree = KDTree(self.rand_source_points) _, idx_source_for_each_target_pt = source_kd_tree.query(self.rand_target_points) # Using the same error as the Matlab code., which is the average euclidean error in n-dimensions. # Paper says that it used the sum of squares. for i in range(self.n_features): for j in range(self.n_features): self.c_spatial[i, j] = np.sqrt(np.sum((self.rand_source_eig_vecs[idx_source_for_each_target_pt, j] - self.rand_target_eig_vecs[:, i])**2) )/self.rand_target_eig_vecs.shape[0] self.c_spatial_f[i, j] = np.sqrt(np.sum((-self.rand_source_eig_vecs[idx_source_for_each_target_pt, j] - self.rand_target_eig_vecs[:, i])**2) )/self.rand_target_eig_vecs.shape[0] def sort_eigenmaps(self): """ Run functions necessary to sort eigenvalues. Seems to do a good job identifying flips. However, if there are abnormal eigenmaps it doesnt do very well. :return: """ self.calc_c_lambda() # lambda = eigenvalues. self.calc_c_hist() self.calc_c_spatial() self.eigen_sort() return self.Q # Q was calculated in eigen_sort and is the cost associated with each of the final pairings.
Methods
def calc_c_hist(self)
-
Compare the histograms of all of the eigenvectors (eigenfunctions) for each of the meshes - only upto the number of meshes we are interested in. The paper says to do this for the number of features we are interested (i.e. 5 or 6), however, this implementation does that for more than the requested number of features to ensure that we dont get the wrong ones. So, we'll get more spectral coordinates, order them, and then select the appropriate number from these re-ordered coordinates.
:return:
Expand source code
def calc_c_hist(self): """ Compare the histograms of all of the eigenvectors (eigenfunctions) for each of the meshes - only upto the number of meshes we are interested in. The paper says to do this for the number of features we are interested (i.e. 5 or 6), however, this implementation does that for more than the requested number of features to ensure that we dont get the wrong ones. So, we'll get more spectral coordinates, order them, and then select the appropriate number from these re-ordered coordinates. :return: """ # Initially tried using straight values (not log) but eig vec 1 got accentuated too much (in the weighting) # Then tried just log, but becuase there are negative values it creates erro. # need to add .5 + a small value to ensure there are no 0 values entered to log. # wasserstein_distance is the same as earth movers distance, and is the minimum "work" needed to # transform u (first entry) into v (second entry). eps = np.finfo(float).eps for i in range(self.n_features): for j in range(self.n_features): self.c_hist[i, j] = wasserstein_distance(np.log(self.rand_target_eig_vecs[:, i] + 0.5 + eps), np.log(self.rand_source_eig_vecs[:, j] + 0.5 + eps)) self.c_hist_f[i, j] = wasserstein_distance(np.log(self.rand_target_eig_vecs[:, i] + 0.5 + eps), np.log(-self.rand_source_eig_vecs[:, j] + 0.5 + eps))
def calc_c_lambda(self)
-
Get average gap between successive eigenvalues (used to normalize any errors/differences between the eigenvalues of the two meshes. :return:
Expand source code
def calc_c_lambda(self): """ Get average gap between successive eigenvalues (used to normalize any errors/differences between the eigenvalues of the two meshes. :return: """ for graph in [self.graph_source, self.graph_target]: if graph.eig_val_gap is None: graph.get_eig_val_gap() # average the gap of mesh 1 and mesh 2. eigen_gap = (self.graph_target.eig_val_gap + self.graph_source.eig_val_gap)/2 # Calculate difference in eigenvalues between meshes normalizing to eigen_gap. for i in range(self.n_features): for j in range(self.n_features): self.c_lambda[i, j] = np.exp((self.graph_target.eig_vals[i] - self.graph_source.eig_vals[j]) ** 2 / (2 * eigen_gap ** 2) )
def calc_c_spatial(self)
-
Next, we calculate the spectral distance of points that are nearest to one another on the mesh. This must assume that the meshes are crudely aligned. This will mean that for knees, we will always have to flip any left knees to match right knees, we should then use an ICP algoritm to get a crude starting alignment. However, even without the ICP (as long as left are flipped) this algorithm scales the sizes of the bounding boxes for the points, so they should be very crudely aligned.
:return:
Expand source code
def calc_c_spatial(self): """ Next, we calculate the spectral distance of points that are nearest to one another on the mesh. This must assume that the meshes are crudely aligned. This will mean that for knees, we will always have to flip any left knees to match right knees, we should then use an ICP algoritm to get a crude starting alignment. However, even without the ICP (as long as left are flipped) this algorithm scales the sizes of the bounding boxes for the points, so they should be very crudely aligned. :return: """ source_kd_tree = KDTree(self.rand_source_points) _, idx_source_for_each_target_pt = source_kd_tree.query(self.rand_target_points) # Using the same error as the Matlab code., which is the average euclidean error in n-dimensions. # Paper says that it used the sum of squares. for i in range(self.n_features): for j in range(self.n_features): self.c_spatial[i, j] = np.sqrt(np.sum((self.rand_source_eig_vecs[idx_source_for_each_target_pt, j] - self.rand_target_eig_vecs[:, i])**2) )/self.rand_target_eig_vecs.shape[0] self.c_spatial_f[i, j] = np.sqrt(np.sum((-self.rand_source_eig_vecs[idx_source_for_each_target_pt, j] - self.rand_target_eig_vecs[:, i])**2) )/self.rand_target_eig_vecs.shape[0]
def eigen_sort(self)
-
- Calculate a metric of similarity c for straight matches, and matches from flipped eigenvectors c_f for source mesh.
- Get the dissimilarity matrix Q by taking the min of c and c_f
- Get matrix S of positions which were flipped to get minimum
- Get list of
row_flipped
andcols_flipped
- Get list of matching columns for each row by using hungarian algorithm (linear_sum_assignment) on Q
- From the matches, and the identified flipped rows/cols, identify which of the pairs are flipped
- Flip the eigen vector values for the necessary eigen vectors on the source mesh
Expand source code
def eigen_sort(self): """ - Calculate a metric of similarity c for straight matches, and matches from flipped eigenvectors c_f for source mesh. - Get the dissimilarity matrix Q by taking the min of c and c_f - Get matrix S of positions which were flipped to get minimum - Get list of `row_flipped` and `cols_flipped` - Get list of matching columns for each row by using hungarian algorithm (linear_sum_assignment) on Q - From the matches, and the identified flipped rows/cols, identify which of the pairs are flipped - Flip the eigen vector values for the necessary eigen vectors on the source mesh """ # find best matching (negative/positive) c = self.c_spatial * self.c_lambda * self.c_hist c_f = self.c_spatial_f * self.c_lambda * self.c_hist_f self.Q = np.min((c, c_f), axis=0) # telling us if sign flipped or not. True = same sign, False = flipped S = c > c_f (target_flipped, source_flipped) = np.where(S == True) # Get matches if self.target_as_reference is True: target_matches, source_matches = linear_sum_assignment(self.Q) elif self.target_as_reference is False: source_matches, target_matches = linear_sum_assignment(self.Q.T) # The original Matlab code calculates Q as: # Q = sum(M'*Q')'; which is equivalent to np.sum(np.matmul(M.T, Q.T), axis=0) # Where M = nXn matrix where n is number of spectral features (or the length of target_matches). and # has 1s where there is a "match", i.e. # M = np.zeros((len(target_matches), len(target_matches))) # M[target_matches, source_matches] = 1 # This results in a permutation cost equal to all of the "costs" for the entire row. This seems wrong, # it seems like the cost of using a particular "pair" should come down to the cost of just that particular # match and not all of the matches that it "didnt" use. In fact, bad matches would have high numbers and screw # things up. We are therefore doing the below. self.Q = self.Q[target_matches, source_matches] # This will identify the source & target eigen #s (pairs after matching) that need to be flipped. # So, if we're on eigen 3 for target, but its pair is eigen 5, then this will identify the pair (3, 5). # The for loop then iterates through and turns the 5 (from the source) to be negative and leave it where # it is. Finally, we will then select the final eigenvecs in the appropriate order (after flipping) # using their indices from target_matches, source_matches. # If a "flipped pair" is also a "best match", then store the pairs in flipped_pairs so that the corresponding # source eig_vecs can be flipped to match the target eig_vec orientation. flipped_pairs = [p2 for p1 in zip(target_flipped, source_flipped) for p2 in zip(target_matches, source_matches) if p2 == p1] # Flipped pairs are used to flip the sign of the source eig_vecs. for mode_0, mode_1 in flipped_pairs: if self.target_as_reference is True: self.graph_source.eig_vecs[:, mode_1] = self.graph_source.eig_vecs[:, mode_1] * -1 elif self.target_as_reference is False: self.graph_target.eig_vecs[:, mode_0] = self.graph_target.eig_vecs[:, mode_0] * -1 # The source eig_vecs are re-ordered to match the order from the target eig_vecs - based on the best matches # Identified using Hungarian algorithm (linear_sum_assignment) on the dissimilarity matrix Q. if self.target_as_reference is True: self.graph_source.eig_vecs[:, target_matches] = self.graph_source.eig_vecs[:, source_matches] elif self.target_as_reference is False: self.graph_target.eig_vecs[:, source_matches] = self.graph_target.eig_vecs[:, target_matches] print_header('Eigenvector Sorting Results') if self.target_as_reference is True: print('Using target eigenmaps as the reference') elif self.target_as_reference is False: print('Using source eigenmaps as the reference') print('The matches for eigenvectors were as follows:') print('Target\t| Source') for matched_pair in zip(target_matches, source_matches): source_value = str(matched_pair[1]) target_value = str(matched_pair[0]) if matched_pair in flipped_pairs: if self.target_as_reference is True: source_value = '-' + source_value elif self.target_as_reference is False: target_value = '-' + target_value print('{:6}\t| {:6}'.format(target_value, source_value)) print('*Negative source values means those eigenvectors were flipped*\n ')
def sort_eigenmaps(self)
-
Run functions necessary to sort eigenvalues. Seems to do a good job identifying flips. However, if there are abnormal eigenmaps it doesnt do very well. :return:
Expand source code
def sort_eigenmaps(self): """ Run functions necessary to sort eigenvalues. Seems to do a good job identifying flips. However, if there are abnormal eigenmaps it doesnt do very well. :return: """ self.calc_c_lambda() # lambda = eigenvalues. self.calc_c_hist() self.calc_c_spatial() self.eigen_sort() return self.Q # Q was calculated in eigen_sort and is the cost associated with each of the final pairings.