<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Graph Analyses
Here, we'll perform various analysis by constructing graphs and measure properties of those graphs to learn more about the data
Step1:  We'll start with just looking at analysis in euclidian space, then thinking about weighing by synaptic density later. Since we hypothesize that our data will show that tissue varies as we move down the y-axis (z-axis in brain) through cortical layers, an interesting thing to do would be compare properties of the graphs on each layer (ie how does graph connectivity vary as we move through layers).
Let's start by triangulating our data. We'll use Delaunay on each y layer first. Putting our data in the right format for doing graph analysis...
Step2:  Now that our data is in the right format, we'll create 52 delaunay graphs. Then we'll perform analyses on these graphs. A simple but useful metric would be to analyze edge length distributions in each layer.
Step3:  We're going to need a method to get edge lengths from 2D centroid pairs
Step4:  Realizing after all this that simply location is useless. We know the voxels are evenly spaced, which means our edge length data will be all the same. See that the "centroids" are no different
Step5:  There is no distance between the two. Therefore it is perhaps more useful to consider a graph that considers node weights. Voronoi is dual to Delaunay, so that's not much of an option. We want something that considers both spacial location and density similarity. 
Region Adjacency Graph (RAG)
In the past we've considered density information alone (ie analysis density histogram distribution) and above we are considering only spacial information, which totally doesn't say anything. To construct a graph that consider both spacial and density information, we'll use a Region Adjacency Graph (RAG).
In RAGs, two nodes are considered as neighbor if they are close in proximity (separated by a small number of pixels/voxels) in the horizontal or vertical direction.
<img src="../docs/figures/plot_rag_1.png" width="400">
Since our data includes density values at each node (voxel, or pixel since we're looking at y-layers), we can weight by the inverse of density difference between two nodes. Inverse because strongly connected nodes should be close in weight.
We have number of synapses S<sub>i</sub> at nodes $i$ and define weights $w$ between the nodes
Step6:  Handrolling Own RAG generator
Step7:  The following method is for later
Step8:  Testing the RAG generator
Step9:  Creating RAGs for each layer
Step10:  OK, great! Now we have a list of 52 region adjacency graphs for each y-layer. Now we want to measure properties of those graphs and see how the properties vary in the y direction - through what we hypothesize are the cortical layers.
Number of Edges
This is just a sanity check. They should all be the same if we did it right!
Step11:  Drawing Graphs
First we look at the default networkx graph plotting
Step12:  This is using the spring layout, so we're losing positional information. We can improve the plot by adding position information.
Edge Weight Stats
Step13:  Full edge weight histograms down y-axis
Step14:  Mean edge weights
Step15:  Edge Weight Variances
Step16:  Edge Weight Third Moments (Skewness)
Step17:  Edge Weight Fourth Moments (Kurtosis)
Step18:  Hmmm...very interesting
Linear graph weights
We're now going to change our weighting function to be linear and scaled by the max and min difference in each layer. This might help eliminates some of the edge effect behavior I suspect is causing that rapid change in statistics in deeper y-layers.
Step19:  Linear Edge Weight Means
Step20:  Linear Edge Weight Variance
Step21:  Linear Edge Weight Skewness
Step22:  Linear Edge Weight Kurtosis
Step23:  Number of Self Loops
Step24:  Interesting. There are no self loops. Why would this be? Let's come back to this. In the meantime, I want to give some though to what it means to have a self loop, whether it should be theoretically possible given our data, and whether our graphs are formed properly.
The answer to this question is very simple. In a RAG, there are no self-loops by definition. Self loops are edges that form a connection between a node and itself.
<img src="../docs/figures/selfloop.png" width="100">
To see whether the graphs are formed properly, let's look at an adjacency lists
Step25:  Compare that to the test data
Step26:  X-Layers
Step27:  We can see here the number of edges is low in that area that does not have many synapses. It, as expected, mirrors the distribution of synapses. It appears to be approximately uniform at the top, with buffers of very few synapses on the sides. Remember from here
<ASSISTANT_TASK:>
Python Code:
import csv
from scipy.stats import kurtosis
from scipy.stats import skew
from scipy.spatial import Delaunay
import numpy as np
import math
import skimage
import matplotlib.pyplot as plt
import seaborn as sns
from skimage import future
import networkx as nx
%matplotlib inline
# Read in the data
data = open('../data/data.csv', 'r').readlines()
fieldnames = ['x', 'y', 'z', 'unmasked', 'synapses']
reader = csv.reader(data)
reader.next()
rows = [[int(col) for col in row] for row in reader]
# These will come in handy later
sorted_x = sorted(list(set([r[0] for r in rows])))
sorted_y = sorted(list(set([r[1] for r in rows])))
sorted_z = sorted(list(set([r[2] for r in rows])))
Explanation: Graph Analyses
Here, we'll perform various analysis by constructing graphs and measure properties of those graphs to learn more about the data
End of explanation
a = np.array(rows)
b = np.delete(a, np.s_[3::],1)
# Separate layers - have to do some wonky stuff to get this to work
b = sorted(b, key=lambda e: e[1])
b = np.array([v.tolist() for v in b])
b = np.split(b, np.where(np.diff(b[:,1]))[0]+1)
Explanation: We'll start with just looking at analysis in euclidian space, then thinking about weighing by synaptic density later. Since we hypothesize that our data will show that tissue varies as we move down the y-axis (z-axis in brain) through cortical layers, an interesting thing to do would be compare properties of the graphs on each layer (ie how does graph connectivity vary as we move through layers).
Let's start by triangulating our data. We'll use Delaunay on each y layer first. Putting our data in the right format for doing graph analysis...
End of explanation
graphs = []
centroid_list = []
for layer in b:
    centroids = np.array(layer)
    
    # get rid of the y value - not relevant anymore
    centroids = np.delete(centroids, 1, 1)
    centroid_list.append(centroids)
    
    graph = Delaunay(centroids)
    graphs.append(graph)
Explanation: Now that our data is in the right format, we'll create 52 delaunay graphs. Then we'll perform analyses on these graphs. A simple but useful metric would be to analyze edge length distributions in each layer.
End of explanation
def get_d_edge_length(edge):
    (x1, y1), (x2, y2) = edge
    return math.sqrt((x2-x1)**2 + (y2-y1)**2)
edge_length_list = [[]]
tri_area_list = [[]]
for del_graph in graphs:
    
    tri_areas = []
    edge_lengths = []
    triangles = []
    for t in centroids[del_graph.simplices]:
        triangles.append(t)
        a, b, c = [tuple(map(int,list(v))) for v in t]
        edge_lengths.append(get_d_edge_length((a,b)))
        edge_lengths.append(get_d_edge_length((a,c)))
        edge_lengths.append(get_d_edge_length((b,c)))
        try:
            tri_areas.append(float(Triangle(a,b,c).area))
        except:
            continue
    edge_length_list.append(edge_lengths)
    tri_area_list.append(tri_areas)
Explanation: We're going to need a method to get edge lengths from 2D centroid pairs
End of explanation
np.subtract(centroid_list[0], centroid_list[1])
Explanation: Realizing after all this that simply location is useless. We know the voxels are evenly spaced, which means our edge length data will be all the same. See that the "centroids" are no different:
End of explanation
real_volume = np.zeros((len(sorted_y), len(sorted_x), len(sorted_z)))
for r in rows:
    real_volume[sorted_y.index(r[1]), sorted_x.index(r[0]), sorted_z.index(r[2])] = r[-1]
np.shape(real_volume)
Explanation: There is no distance between the two. Therefore it is perhaps more useful to consider a graph that considers node weights. Voronoi is dual to Delaunay, so that's not much of an option. We want something that considers both spacial location and density similarity. 
Region Adjacency Graph (RAG)
In the past we've considered density information alone (ie analysis density histogram distribution) and above we are considering only spacial information, which totally doesn't say anything. To construct a graph that consider both spacial and density information, we'll use a Region Adjacency Graph (RAG).
In RAGs, two nodes are considered as neighbor if they are close in proximity (separated by a small number of pixels/voxels) in the horizontal or vertical direction.
<img src="../docs/figures/plot_rag_1.png" width="400">
Since our data includes density values at each node (voxel, or pixel since we're looking at y-layers), we can weight by the inverse of density difference between two nodes. Inverse because strongly connected nodes should be close in weight.
We have number of synapses S<sub>i</sub> at nodes $i$ and define weights $w$ between the nodes:
$$w = \dfrac{1}{S_i - S_{i+1}}$$
RAGs are used largely in image processing, and it makes sense for our data to look more like an image. Since the data is evenly spaced, the absolute locations of the voxels don't matter. We can use the index in the matrix to represent spacial location, with the value at each pixel being the synapse density at that voxel. We've done this before in "real volume"
End of explanation
# point = tuple containing index of point (position)
# returns list of neighbors in [north, east, south, west]
def get_neighbors(point, image):
    shape = np.shape(image)
    neighbors = []
    
    # North
    neighbors.append((point[0], point[1]-1)) if point[1]>0 else neighbors.append(None)
    # East
    neighbors.append((point[0]+1, point[1])) if point[0]<shape[0]-1 else neighbors.append(None)
    # South
    neighbors.append((point[0], point[1]+1)) if point[1]<shape[1]-1 else neighbors.append(None)
    # West
    neighbors.append((point[0]-1, point[1])) if point[0]>0 else neighbors.append(None)
    return neighbors
# calculates weights between nodes
# weight defined as inverse absolute distance
def get_weights_nonlinear(image, point, neighbors):
    weights = []
    for neigh in neighbors:
        if neigh != None:
            weight = 1/(abs(image[point] - image[neigh])+1)
            weights.append(weight)
    return weights
Explanation: Handrolling Own RAG generator
End of explanation
# calculates weights between nodes
# weight scaled and linear
# TODO: Explain weighting difference with math
def get_weights_linear(image, point, neighbors, scale_factor):
    
    weights = []
    for neigh in neighbors:
        if neigh != None:
            diff = abs(image[point] - image[neigh])
            weight = 1 - (scale_factor*diff)
            weights.append(weight)
    
    return weights
image = real_volume[1]
# print image
point = (1,1)
neighbors = get_neighbors(point, image)
# print neighbors
ws = get_weights_nonlinear(image, point, neighbors)
def populate_graph(G, im):
    nodes_to_add = []
    for x in range(np.shape(im)[0]):
        for y in range(np.shape(im)[1]):
            nodes_to_add.append((x,y))
            
    G.add_nodes_from(nodes_to_add)
def get_diff_range(image):
    diffs = []
    
    x = 0
    for col in image:
        y = 0
        for pix in col:
            point = (x,y)
            neighs = get_neighbors(point, image)
            
            for neigh in neighbors:
                if neigh != None:
                    diffs.append(abs(image[point] - image[neigh]))
            
            y+=1
        x+=1  
    return (max(diffs), min(diffs))
def generate_rag(im, linear):
    G=nx.Graph()
    
    if linear == True:
        (max_diff, min_diff) = get_diff_range(im)
        scale_factor = 1/(max_diff - min_diff)
    
    x = 0
    for col in im:
        y = 0
        for pix in col:
            point = (x,y)
            neighs = get_neighbors(point, im)
            
            if linear == True:
                weights = get_weights_linear(im, point, neighs, scale_factor)
            else:
                weights = get_weights_nonlinear(im, point, neighs)
            
            to_add = []
            which = 0
            for neigh in neighs:
                if neigh != None:
                    to_add.append((point, neigh, weights[which]))
                    which+=1
#             print to_add
            G.add_weighted_edges_from(to_add)
            y+=1
        x+=1
    return G
Explanation: The following method is for later
End of explanation
test_im = real_volume[1]
shape = np.shape(test_im)
ragu = generate_rag(test_im, False)
ragu.number_of_edges()
# ragu.adjacency_list()
Explanation: Testing the RAG generator
End of explanation
y_rags = []
for layer in real_volume:
    y_rags.append(generate_rag(layer, False))
Explanation: Creating RAGs for each layer
End of explanation
num_edges = []
for rag in y_rags:
    num_edges.append(rag.number_of_edges())
sns.barplot(x=range(len(num_edges)), y=num_edges)
sns.plt.show()
Explanation: OK, great! Now we have a list of 52 region adjacency graphs for each y-layer. Now we want to measure properties of those graphs and see how the properties vary in the y direction - through what we hypothesize are the cortical layers.
Number of Edges
This is just a sanity check. They should all be the same if we did it right!
End of explanation
# for rag in y_rags:
#     plt.figure()
#     nx.draw(rag, node_size=100)
Explanation: Drawing Graphs
First we look at the default networkx graph plotting:
We're gonna need to massage the drawing function a bit in order to get the custom positional graph to work.
End of explanation
def get_edge_weight_distributions(rags):
    distributions = []
    for rag in rags:    
        itty = rag.edges_iter()
        weight_list = []
        for index in range(rag.number_of_edges()):
            eddy = itty.next()
            weight_list.append(rag.get_edge_data(eddy[0], eddy[1])['weight'])
        distributions.append(weight_list)
    return distributions
distributions = get_edge_weight_distributions(y_rags)
Explanation: This is using the spring layout, so we're losing positional information. We can improve the plot by adding position information.
Edge Weight Stats
End of explanation
count = 0
for distr in distributions:
    plt.hist(distr, bins=150)
    plt.title("Layer " + str(count) + " Edge Weight Histogram")
    plt.show()
    count+=1
Explanation: Full edge weight histograms down y-axis
End of explanation
y_edge_means = []
for distrib in distributions:
    y_edge_means.append(np.mean(distrib)) 
print y_edge_means
sns.barplot(x=range(len(y_edge_means)), y=y_edge_means)
sns.plt.show()
Explanation: Mean edge weights
End of explanation
y_edge_vars = []
for distrib in distributions:
    y_edge_vars.append(np.var(distrib)) 
print y_edge_vars
sns.barplot(x=range(len(y_edge_vars)), y=y_edge_vars)
sns.plt.show()
Explanation: Edge Weight Variances
End of explanation
y_edge_skews = []
for distrib in distributions:
    y_edge_skews.append(skew(distrib)) 
print y_edge_skews
sns.barplot(x=range(len(y_edge_skews)), y=y_edge_skews)
sns.plt.show()
Explanation: Edge Weight Third Moments (Skewness)
End of explanation
y_edge_kurts = []
for distrib in distributions:
    y_edge_kurts.append(kurtosis(distrib)) 
print y_edge_kurts
sns.barplot(x=range(len(y_edge_kurts)), y=y_edge_kurts)
sns.plt.show()
Explanation: Edge Weight Fourth Moments (Kurtosis)
End of explanation
y_rags_linear_weight = []
for layer in real_volume:
    y_rags_linear_weight.append(generate_rag(layer, True))
test_rag = generate_rag(real_volume[4], True)
itty = test_rag.edges_iter()
weight_list = []
for index in range(test_rag.number_of_edges()):
    eddy = itty.next()
    weight_list.append(test_rag.get_edge_data(eddy[0], eddy[1])['weight'])
distributions_lin = get_edge_weight_distributions(y_rags_linear_weight)
Explanation: Hmmm...very interesting
Linear graph weights
We're now going to change our weighting function to be linear and scaled by the max and min difference in each layer. This might help eliminates some of the edge effect behavior I suspect is causing that rapid change in statistics in deeper y-layers.
End of explanation
y_edge_linear_means = []
for distrib in distributions_lin:
    y_edge_linear_means.append(np.mean(distrib)) 
sns.barplot(x=range(len(y_edge_linear_means)), y=y_edge_linear_means)
sns.plt.show()
Explanation: Linear Edge Weight Means
End of explanation
y_edge_linear_vars = []
for distrib in distributions_lin:
    y_edge_linear_vars.append(np.var(distrib)) 
sns.barplot(x=range(len(y_edge_linear_vars)), y=y_edge_linear_vars)
sns.plt.show()
Explanation: Linear Edge Weight Variance
End of explanation
y_edge_linear_skews = []
for distrib in distributions_lin:
    y_edge_linear_skews.append(skew(distrib)) 
sns.barplot(x=range(len(y_edge_linear_skews)), y=y_edge_linear_skews)
sns.plt.show()
Explanation: Linear Edge Weight Skewness
End of explanation
y_edge_linear_kurts = []
for distrib in distributions_lin:
    y_edge_linear_kurts.append(kurtosis(distrib)) 
sns.barplot(x=range(len(y_edge_linear_kurts)), y=y_edge_linear_kurts)
sns.plt.show()
Explanation: Linear Edge Weight Kurtosis
End of explanation
num_self_loops = []
for rag in y_rags:
    num_self_loops.append(rag.number_of_selfloops())
num_self_loops
Explanation: Number of Self Loops
End of explanation
# y_rags[0].adjacency_list()
Explanation: Interesting. There are no self loops. Why would this be? Let's come back to this. In the meantime, I want to give some though to what it means to have a self loop, whether it should be theoretically possible given our data, and whether our graphs are formed properly.
The answer to this question is very simple. In a RAG, there are no self-loops by definition. Self loops are edges that form a connection between a node and itself.
<img src="../docs/figures/selfloop.png" width="100">
To see whether the graphs are formed properly, let's look at an adjacency lists:
End of explanation
# Test Data
test = np.array([[1,2],[3,4]])
test_rag = skimage.future.graph.RAG(test)
test_rag.adjacency_list()
Explanation: Compare that to the test data:
End of explanation
real_volume_x = np.zeros((len(sorted_x), len(sorted_y), len(sorted_z)))
for r in rows:
    real_volume_x[ sorted_x.index(r[0]), sorted_y.index(r[1]), sorted_z.index(r[2])] = r[-1]
x_rags = []
count = 0;
for layer in real_volume_x:
    count = count + 1
    x_rags.append(skimage.future.graph.RAG(layer))
num_edges_x = []
for rag in x_rags:
    num_edges_x.append(rag.number_of_edges())
sns.barplot(x=range(len(num_edges_x)), y=num_edges_x)
sns.plt.show()
Explanation: X-Layers
End of explanation
plt.imshow(np.amax(real_volume, axis=2), interpolation='nearest')
plt.show()
# edge_length_list[3]
# tri_area_list[3]
# triangles
# Note for future
# del_features['d_edge_length_mean'] = np.mean(edge_lengths)
# del_features['d_edge_length_std'] = np.std(edge_lengths)
# del_features['d_edge_length_skew'] = scipy.stats.skew(edge_lengths)
# del_features['d_edge_length_kurtosis'] = scipy.stats.kurtosis(edge_lengths)
Explanation: We can see here the number of edges is low in that area that does not have many synapses. It, as expected, mirrors the distribution of synapses. It appears to be approximately uniform at the top, with buffers of very few synapses on the sides. Remember from here:
End of explanation
<END_TASK>