from sklearn.decomposition import PCA

alpha = 0.99

n = 0
pca = PCA(n_components=n)
ratios = [0]
while sum(ratios) < alpha:
    n += 1
    pca = PCA(n_components=n)
    pca.fit(frame_z)
    ratios = pca.explained_variance_ratio_

print "{}% accounted for with n={} PCs".format(alpha, n)

"""
Explanation: 1. PCA
Let us begin by running PCA — we're 'satisfied' once we find $n$ st $sum(axis ratios) > \alpha \%$
End of explanation
"""


y_sum = [0] * len(vol[0,:,0])
for i in range(len(vol[0,:,0])):
    y_sum[i] = sum(sum(vol[:,i,:]))

sns.barplot(x=range(len(y_sum[:40])), y=y_sum[:40])

sns.distplot(y_sum, bins=len(y_sum))

"""
Explanation: 2. Orientation in Brain Space
We believe this to be a sample of V1, and as such, we believe this to be highly conformative to the layering structure of mammalian cortex.  To orient ourselves in cortical space, we anticipate finding a pial surface. Using Bock et al. Nature (2011), we can establish our bearings using ndviz.
This epithelial surface is clearly visible in this view, indicating that our $y$ axis is the cortical depth axis.
I expect learning more about PCA to get a better direction for axis than simply snapping to a cardinal axis.

3. Layer I Delimiterization
We believe that this dataset extends through layers I, II, III, and possibly part of IV. To determine this, we rely on the fact that mammalian cortex Layer I is molecular, and is thus low in synaptic density, compared to layer II. So there should be clear delimiterization along our $y$ axis:
End of explanation
"""


sns.barplot(x=range(len(y_sum[:40])), y=y_sum[:40])

"""
Explanation: Above, we see a histogram of y_sum that indicates that there is a local minimum at the 12th layer of y-sampling, which colocates with where we anticipate seeing the boundary between layers I and II. Here is the biological substantiation:

As we can see, at about 1/3 of the 'depth' into cortex is the boundary to layer II.
4. Local maxima to find other cortical boundaries
Consider  these lamination illustrations (Ramon y Cajal) that show a 'wavelike' pattern to cellular density (and thus inversely, synaptic density). We see that our histogram follows this pattern:
End of explanation
"""


from scipy.signal import argrelextrema
def local_minima(a):
    return argrelextrema(a, np.less)

whole_volume_minima = local_minima(np.array(y_sum))

whole_volume_minima

"""
Explanation: Using these local maxima as delimiters, we can see a boundary at (12, 18, 24, 30), and then what is likely the boundary of Layer IVa/b at 33.
5. Do these layers stay approximately the same on different subvolumes?
Let's examine the layering, reviewed above:
End of explanation
"""


y_sum_left = [0] * len(vol[0,:,5:])
for i in range(len(vol[0,:,5:])):
    y_sum_left[i] = sum(sum(vol[:,i,5:]))
    
left_volume_minima = local_minima(np.array(y_sum_left))

y_sum_right = [0] * len(vol[0,:,:5])
for i in range(len(vol[0,:,:5])):
    y_sum_right[i] = sum(sum(vol[:,i,:5]))
    
right_volume_minima = local_minima(np.array(y_sum_right))

left_volume_minima, right_volume_minima

"""
Explanation: Now let's examine halves of the volume:
End of explanation
"""


