Add files via upload

Author: Ribosome25

Toolbox.py

@@ -0,0 +1,103 @@
+import numpy as np
+import pandas as pd
+import os
+from scipy.stats import pearsonr
+from sklearn.manifold import MDS
+
+def two_d_norm(xy):
+    # xy is N x 2 xy cordinates, returns normed-xy on [0,1]
+    norm_xy = (xy - xy.min(axis = 0)) / (xy - xy.min(axis = 0)).max(axis = 0)
+    return norm_xy
+
+def two_d_eq(xy,Nn):
+    # xy is N x 2 xy cordinates, returns eq-xy on [0,1]
+    xx_rank = np.argsort(xy[:,0])
+    yy_rank = np.argsort(xy[:,1])
+    eq_xy = np.full(xy.shape,np.nan)
+    for ii in range(xy.shape[0]):
+        xx_idx = xx_rank[ii]
+        yy_idx = yy_rank[ii]
+        eq_xy[xx_idx,0] = ii * 1/Nn
+        eq_xy[yy_idx,1] = ii * 1/Nn
+    return eq_xy
+
+#embedding = MDS(n_components=2)
+#mds_xy = embedding.fit_transform(transposed_input)
+# to pixels
+def Assign_features_to_pixels(xy,nn,verbose = False):
+    # For each unassigned feature, find its nearest pixel, repeat until every ft is assigned
+    # xy is the 2-d coordinates (normalized to [0,1]); nn is the image width. Img size = n x n
+    # generate the result summary table, xy pixels; 3rd is nan for filling the idx
+    Nn = xy.shape[0]
+    
+    from itertools import product
+    pixel_iter = product([x for x in range(nn)],repeat = 2)
+    result_table = np.full((nn*nn,3),np.nan)
+    ii = 0
+    for each_pixel in pixel_iter:
+        result_table[ii,:2] = np.array(each_pixel)
+        ii+=1
+    # Use numpy array for speed 
+        
+    from sklearn.metrics import pairwise_distances
+    
+#    xy = eq_xy
+    centroids = result_table[:,:2] / nn + 0.5/nn
+    pixel_avail = np.ones(nn*nn).astype(bool)
+    feature_assigned = np.zeros(Nn).astype(bool)
+    
+    dist_xy_centroids = pairwise_distances(centroids,xy,metric='euclidean')
+    
+    while feature_assigned.sum()<Nn:
+        # Init the pick-relationship table
+        pick_xy_centroids = np.zeros(dist_xy_centroids.shape).astype(bool)
+        
+        for each_feature in range(Nn):
+            # for each feature, find the nearest available pixel
+            if feature_assigned[each_feature] == True:
+                # if this feature is already assigned, skip to the next ft
+                continue
+            else:
+                # if not assigned:
+                for ii in range(nn*nn):
+                    # find the nearest avail pixel
+                    nearest_pixel_idx = np.argsort(dist_xy_centroids[:,each_feature])[ii]
+                    if pixel_avail[nearest_pixel_idx] == True:
+                        break
+                    else:
+                        continue
+                pick_xy_centroids[nearest_pixel_idx,each_feature] = True
+            
+        for each_pixel in range(nn*nn):
+            # Assign the feature No to pixels
+            if pixel_avail[each_pixel] == False:
+                continue
+            else:
+                # find all the "True" features. np.where returns a tuple size 1
+                related_features = np.where(pick_xy_centroids[each_pixel,:] == 1)[0]
+                if len(related_features) == 1:
+                    # Assign it
+                    result_table[each_pixel,2] = related_features[0]
+                    pixel_avail[each_pixel] = False
+                    feature_assigned[related_features[0]] = True
+                elif len(related_features) > 1:
+                    related_dists = dist_xy_centroids[each_pixel,related_features]
+                    best_feature = related_features[np.argsort(related_dists)[0]] # Sort, and pick the nearest one among them
+                    result_table[each_pixel,2] = best_feature
+                    pixel_avail[each_pixel] = False
+                    feature_assigned[best_feature] = True
+        if verbose:
+            print(">> Assign features to pixels:", feature_assigned.sum(),"/",Nn)
+    result_table = result_table.astype(int)
+    
+    img = np.full((nn,nn),'NaN').astype(object)
+    for each_pixel in range(nn*nn):
+        xx = result_table[each_pixel,0]
+        yy = result_table[each_pixel,1]
+        ft = 'F' + str(result_table[each_pixel,2])
+        img[xx,yy] = ft
+    return img.astype(object)
+print(">>>> MDS")
+#eq_xy = two_d_eq(mds_xy)
+#Img = Assign_features_to_pixels(eq_xy,nn,verbose=1)
+#Init_Corr_MDS = InitCorr(dist_mat,Img,nn)

paraHill.py

@@ -0,0 +1,132 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Aug 12 14:30:06 2019
+
+@author: Ruibzhan & Omid Bazgir
+"""
+
+from scipy.stats import pearsonr
+import numpy as np
+import random
+from scipy.spatial import distance
+import pickle
+import pandas as pd
+import time
+from itertools import product
+
+#%%
+def universial_corr(dist_matr, mapping_in_int):
+    # dist_matr is a sqr matr
+    Nn = dist_matr.shape[0]
+    # find what is the int coordinates for each feature, get an array
+    # Because np.where returns a tuple (x_position array, y_position array), a generation is used
+    coord = np.array([[item[0] for item in np.where(mapping_in_int == ii)] for ii in range(Nn)])
+    # get a 1-d form distance the euclidean dist between pixles positions
+    pixel_dist = distance.pdist(coord)
+    pixel_dist = pixel_dist.reshape(len(pixel_dist),1)
+    # convert the 2-d distance to 1d distance
+    feature_dist = distance.squareform(dist_matr)
+    feature_dist = feature_dist.reshape(len(feature_dist),1)
+    ## pearsonr returns a tuple
+    #corrs = pearsonr(feature_dist,pixel_dist)[0]
+    L2_Norm = np.sqrt(sum((pixel_dist - feature_dist)**2)/sum(feature_dist**2))
+    return L2_Norm
+#%%
+def evaluate_swap(coord1,coord2,dist_matr,mapping_in_int,original_corr = -2):
+    # Coord are in list[]
+    # Avoid messing up with the origianl map
+    # The correlation before swap can be passed to save some calculation
+    the_map = mapping_in_int.copy()
+    # If out of bound, return NaN. 
+    if coord1[0]<0 or coord1[1]<0 or coord2[0]<0 or coord2[1]<0:
+        return np.nan
+    if coord1[0]>=the_map.shape[0] or coord1[1]>=the_map.shape[0] or coord2[0]>=the_map.shape[0] or coord2[1]>=the_map.shape[0]:
+        return np.nan
+    # If not given, recompute.
+    if original_corr<-1 or original_corr>1:
+        original_corr = universial_corr(dist_matr,the_map)
+    # Swap
+    try:
+        temp = the_map[coord1[0],coord1[1]] 
+        the_map[coord1[0],coord1[1]] = the_map[coord2[0],coord2[1]]
+        the_map[coord2[0],coord2[1]] = temp
+        changed_corr = universial_corr(dist_matr,the_map)
+        return(changed_corr - original_corr)
+    except IndexError:
+        raise Warning ("Swap index:", coord1,coord2,"Index error. Check the coordnation.")
+        return np.nan
+    
+def evaluate_centroid(centroid,dist_matr,mapping_in_int):
+    original_corr = universial_corr(dist_matr,mapping_in_int)
+    results = [100000] # just to skip the 0 position
+    for each_direc in product([-1,0,1],repeat = 2):
+        #print(each_direc)
+        # directions are returned as tuple (-1,1), (-1,0), (-1,1), (0,0), ....
+        swap_coord = [centroid[0]+each_direc[0],centroid[1]+each_direc[1]]
+        evaluation = evaluate_swap(centroid,swap_coord,dist_matr,mapping_in_int,original_corr)
+        results.append(evaluation)
+    results_array = np.array(results)
+    #best_swap_direc = np.where(results_array == np.nanmax(results_array))[0][0]
+    best_swap_direc = np.where(results_array == np.nanmin(results_array))[0][0]
+    # Give the best direction as a int
+    return best_swap_direc
+
+def evaluate_centroids_in_list(centroids_list,dist_matr,mapping_in_int):
+    # and returns a dict
+    results = dict()
+    for each_centr in centroids_list:
+        each_centr = tuple(each_centr)
+        evaluation = evaluate_centroid(each_centr,dist_matr,mapping_in_int)
+        results.update({each_centr:evaluation})
+    return results
+
+#%%
+def execute_coordination_swap(coord1,coord2,mapping_in_int):
+    # try passing the ref. directly 
+    the_map = mapping_in_int#.copy()
+    # If out of bound, return NaN. 
+    if coord1[0]<0 or coord1[1]<0 or coord2[0]<0 or coord2[1]<0:
+        raise Warning("Swapping failed:",coord1,coord2,"-- Negative coordnation.")
+        return the_map
+    if coord1[0]>the_map.shape[0] or coord1[1]>the_map.shape[0] or coord2[0]>the_map.shape[0] or coord2[1]>the_map.shape[0]:
+        raise Warning("Swapping failed:",coord1,coord2,"-- Coordnation out of bound.")
+        return the_map
+
+    temp = the_map[coord1[0],coord1[1]] 
+    the_map[coord1[0],coord1[1]] = the_map[coord2[0],coord2[1]]
+    the_map[coord2[0],coord2[1]] = temp
+
+    return(the_map)
+
+# Initial centriod id & Swapping direction: 
+# 1 2 3
+# 4 5 6
+# 7 8 9
+# 0 in swapping is preserved for the header.
+
+def execute_direction_swap(centroid,mapping_in_int,direction = 5):
+    # Need to notice that [0] is the vertival coord, [1] is the horiz coord. similar to the matlab images.
+    coord1 = list(centroid)
+    coord2 = list(centroid)
+    if direction not in range(1,10):
+        raise ValueError("Invalid swapping direction.")
+    if direction == 5:
+        return mapping_in_int
+
+    if direction in [1,4,7]:
+        coord2[1] -=1
+    elif direction in [3,6,9]:
+        coord2[1] +=1
+
+    if direction in [1,2,3]:
+        coord2[0] -=1
+    elif direction in [7,8,9]:
+        coord2[0] +=1
+
+    the_map = execute_coordination_swap(coord1,coord2,mapping_in_int)
+    return the_map
+
+def execute_dict_swap(swapping_dict, mapping_in_int):
+    for each_key in swapping_dict:
+        execute_direction_swap(each_key,mapping_in_int,direction = swapping_dict[each_key])
+    return mapping_in_int