NBLAST SCORE

#!/usr/bin/env python
#coding : utf-8

"""
Created on Jun 15 15:28:03 2016

========================================================================
NBLAST SCORE
========================================================================

This is a sensitive and rapid algorithm for measuring pairwise neuronal
similarity. If you want run this code to get the score, please put all
files in the same folder, and let the gold swc file be the first file
in the folder.

========================================================================
"""

print(__doc__)

import numpy as np
import pandas as pd
import os
import math
os.chdir("E:\gold vs test")

# read gold files and test files in the same folder
def read_file():
L1 = []
L2 = []
L3 = []
for files in os.walk("E:\gold vs test") :
for file in files:
L1.append(file)
L2 = L1[2]
print "There are %d files" %len(L2),'\n'
#for i in range(0,len(L2)):
#   print L2[i]
print '\n',"Read over!"
names = ['id','type','x','y','z','r','parent_id']
f0 = pd.read_csv(L2[0], delimiter=" ", names=names, skiprows=[0])
for i in range(1,len(L2)):
fi = pd.read_csv(L2[i], delimiter = " ", names=names, skiprows=[0])
score = get_score(f0,fi)
L3.append(score)
with open('e:\demo.txt','w') as f:
for i in L3:
f.write(str(i)+'\n')

# get the parameter of u, use gold swc file
def get_u(f,i,point):
for j in range(1,len(f)):
if f.ix[j,'parent_id'] == f.ix[i,'id']:
point2 = np.array(f.ix[j,['x','y','z']])
u = (point-point2)/(sum((point-point2)**2) ** 0.5)
break
else:
u = np.array([0,0,0])
return u

# get the parameter of d and v, use gold swc file and one of comparison files
def get_d(f,point):
distsum = {}
for j in range(0,len(f)):
point_raw = np.array(f.ix[j,['x','y','z']])
dist_raw = sum((point-point_raw)**2) ** 0.5
distsum[dist_raw] = point_raw
[(k,distsum[k]) for k in sorted(distsum.keys())]
L1 = distsum.keys()
d = L1[0]
L2 = distsum.values()
point2 = L2[0]
for k in range(0,len(f)):
if sum(np.array(f.ix[k,['x','y','z']])-point2) == 0:
v = get_u(f,k,point2)
return v, d

# calculate the score fo two 'swc' file
def get_score(f1,f2):
score = 0
num=0
for i in range(0,len(f1)):
point1 = np.array(f1.ix[i,['x','y','z']])
u = get_u(f1,i,point1)
w = get_d(f2,point1)
v = w[0]
d = w[1]
s = (abs(sum(u*v))*math.exp((-d**2*1.0)/(2*3**2))) ** 0.5
score = score + s
print score
return score

if __name__ == "__main__":
read_file()