In [1]:
import os
PATH="/Users/david/Workspaces/claritycontrol/code/" # use your own path
os.chdir(PATH)
import numpy as np
import matplotlib
matplotlib.use('AGG') # avoid some error in matplotlib
import matplotlib.pyplot as plt
import jgraph as ig
%matplotlib inline
# Load histograms
histograms = np.genfromtxt('./data/hist/histograms.csv', delimiter=',')
# Load roi featurs
fields = ['label','nbrain','nroi','roix','roiy','roiy','mean','std','energy','entropy','correlation','contrast','variance',
'summean','inertia','cluster shade','cluster tendency','homogeneity','maxProbability','inverse variance']
data = np.genfromtxt('./data/roi_features/features.csv', delimiter=",", dtype=np.float32)# the features data have been pre-processed and merged
brain_nums = np.unique(data[:,1])
roi_nums = np.unique(data[:,2])
In [3]:
vectorized = histograms[98:107,:].T
covar = np.cov(vectorized)
plt.figure(figsize=(7,7))
plt.imshow(covar)
plt.title('Covariance of CLARITY Histograms', fontsize=16)
plt.xlabel('data', fontsize=16)
plt.ylabel('data', fontsize=16)
plt.colorbar()
plt.show()
diag = covar.diagonal()*np.eye(covar.shape[0])
hollow = covar-diag
d_det = np.linalg.det(diag)
h_det = np.linalg.det(hollow)
plt.figure(figsize=(11,8))
plt.subplot(121)
plt.imshow(diag)
plt.clim([0, np.max(covar)])
plt.title('Determinant of on-diagonal: ' + str(d_det), fontsize=16)
plt.xlabel('data', fontsize=16)
plt.ylabel('data', fontsize=16)
plt.subplot(122)
plt.imshow(hollow)
plt.clim([0, np.max(covar)])
plt.title('Determinant of off-diagonal: ' + str(h_det), fontsize=16)
plt.xlabel('data', fontsize=16)
plt.ylabel('data', fontsize=16)
plt.show()
print "Ratio of on- and off-diagonal determinants: " + str(d_det/h_det)
In [4]:
ratio=[]
for ROI_NUM in roi_nums:
vectorized = data[data[:,2]==ROI_NUM, 6:]
covar = np.cov(vectorized)
plt.figure(figsize=(7,7))
plt.imshow(covar)
plt.title('Covariance of Clarity ROI (%s) features'%(ROI_NUM), fontsize=16)
plt.xlabel('data', fontsize=16)
plt.ylabel('data', fontsize=16)
plt.colorbar()
plt.show()
diag = covar.diagonal()*np.eye(covar.shape[0])
hollow = covar-diag
d_det = np.linalg.det(diag)
h_det = np.linalg.det(hollow)
plt.figure(figsize=(11,8))
plt.subplot(121)
plt.imshow(diag)
plt.clim([0, np.max(covar)])
plt.title('Determinant of on-diagonal: ' + str(d_det), fontsize=16)
plt.xlabel('data', fontsize=16)
plt.ylabel('data', fontsize=16)
plt.subplot(122)
plt.imshow(hollow)
plt.clim([0, np.max(covar)])
plt.title('Determinant of off-diagonal: ' + str(h_det), fontsize=16)
plt.xlabel('data', fontsize=16)
plt.ylabel('data', fontsize=16)
plt.show()
print "Ratio of on- and off-diagonal determinants: " + str(d_det/h_det)
ratio.append(d_det/h_det)
print ratio
n=ratio.index(max(ratio))
print "The strongest indepent roi is [%s]%s"%(n,roi_nums[n])
In [2]:
from sklearn import cross_validation
from sklearn.cross_validation import LeaveOneOut
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.svm import SVC
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis
%matplotlib inline
np.random.seed(12345678) # for reproducibility, set random seed
names = ["Nearest Neighbors", "Linear SVM", "Random Forest",
"Linear Discriminant Analysis", "Quadratic Discriminant Analysis"]
classifiers = [
KNeighborsClassifier(3),
SVC(kernel="linear", C=0.5),
RandomForestClassifier(max_depth=5, n_estimators=10, max_features=1),
LinearDiscriminantAnalysis()]
In [15]:
accuracy = np.zeros((len(classifiers), 2), dtype=np.dtype('float64'))
X=histograms[98:107,:].T
y=np.repeat([0,1,2],5)
for idx, cla in enumerate(classifiers):
X_train, X_test, y_train, y_test = cross_validation.train_test_split(X, y, test_size=0.2, random_state=0)
clf = cla.fit(X_train, y_train)
loo = LeaveOneOut(len(X))
scores = cross_validation.cross_val_score(clf, X, y, cv=loo)
accuracy[idx,] = [scores.mean(), scores.std()]
print("Accuracy of %s: %0.2f (+/- %0.2f)" % (names[idx], scores.mean(), scores.std() * 2))
# plt.errorbar(S, accuracy[:,0,0], yerr = accuracy[:,0,1], hold=True, label=names[0])
# plt.errorbar(S, accuracy[:,1,0], yerr = accuracy[:,1,1], color='green', hold=True, label=names[1])
# plt.errorbar(S, accuracy[:,2,0], yerr = accuracy[:,2,1], color='red', hold=True, label=names[2])
# plt.errorbar(S, accuracy[:,3,0], yerr = accuracy[:,3,1], color='black', hold=True, label=names[3])
# plt.errorbar(S, accuracy[:,4,0], yerr = accuracy[:,4,1], color='brown', hold=True, label=names[4])
# plt.xscale('log')
# plt.xlabel('number of samples')
# plt.ylabel('accuracy')
# plt.title('Accuracy of classification under simulated data')
# plt.axhline(1, color='red', linestyle='--')
# plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
# plt.show()
In [ ]:
In [ ]: