In [2]:
from sympy import *
init_printing()
from IPython.display import display
%matplotlib inline
import matplotlib.pyplot as plt
In [3]:
oxp = Symbol("Omega_x'")
b = Symbol("b")
n = Symbol("n")
theta = Symbol("theta")
s = Symbol("s")
a = Symbol("a")
subsampledOmega = (binomial(s, b) * binomial(n - s, a - b)) / binomial(n, a)
subsampledFpF = Sum(subsampledOmega, (b, theta, s))
subsampledOmegaSlow = (binomial(s, b) * binomial(n - s, a - b))
subsampledFpFSlow = Sum(subsampledOmegaSlow, (b, theta, s))/ binomial(n, a)
display(subsampledFpF)
display(subsampledFpFSlow)
where n refers to the size of the population of cells, a is the number of active cells at any instance in time, s is the number of actual synapses on a dendritic segment, and θ is the threshold for NMDA spikes. Following (Ahmad & Hawkins, 2015), the numerator counts the number of possible ways θ or more cells can match a fixed set of s synapses. The denominator counts the number of ways a cells out of n can be active.
In [4]:
display("n=10000, a=64, s=24, theta=12", subsampledFpF.subs(s,24).subs(n, 10000).subs(a, 64).subs(theta, 12).evalf())
In [5]:
display("n=10000, a=300, s=24, theta=12", subsampledFpFSlow.subs(theta, 12).subs(s, 24).subs(n, 10000).subs(a, 300).evalf())
In [6]:
display("n=2000, a=40, s=20, theta=10", subsampledFpF.subs(theta, 15).subs(s, 20).subs(n, 2000).subs(a, 40).evalf(100))
In [7]:
T1B = subsampledFpFSlow.subs(n, 100000).subs(a, 2000).subs(theta,s).evalf()
print "n=100000, a=2000, theta=s"
display("s=6",T1B.subs(s,6).evalf())
display("s=8",T1B.subs(s,8).evalf())
display("s=10",T1B.subs(s,10).evalf())
In [8]:
T1C = subsampledFpFSlow.subs(n, 100000).subs(a, 2000).subs(s,2*theta).evalf()
print "n=10000, a=300, s=2*theta"
display("theta=6",T1C.subs(theta,6).evalf())
display("theta=8",T1C.subs(theta,8).evalf())
display("theta=10",T1C.subs(theta,10).evalf())
display("theta=12",T1C.subs(theta,12).evalf())
In [9]:
m = Symbol("m")
T1D = subsampledFpF.subs(n, 100000).subs(a, 2000).subs(s,2*m*theta).evalf()
print "n=100000, a=2000, s=2*m*theta"
display("theta=10, m=2",T1D.subs(theta,10).subs(m,2).evalf())
display("theta=10, m=4",T1D.subs(theta,10).subs(m,4).evalf())
display("theta=10, m=6",T1D.subs(theta,10).subs(m,6).evalf())
display("theta=20, m=6",T1D.subs(theta,20).subs(m,6).evalf())
In [10]:
eq1 = subsampledFpFSlow.subs(s, 64).subs(theta, 12)
print "a=64 cells active, s=16 synapses on segment, dendritic threshold is theta=8\n"
errorList = []
nList = []
for n0 in range(300,20100,200):
error = eq1.subs(n, n0).subs(a,64).evalf()
errorList += [error]
nList += [n0]
print "population n = %5d, sparsity = %5.2f%%, probability of false match = "%(n0, 64/n0), error
print errorList
print nList
In [11]:
print ("2% sparsity with n=400")
print subsampledFpFSlow.subs(s, 4).subs(a, 8).subs(theta, 2).subs(n,400).evalf()
print ("2% sparsity with n=4000")
print subsampledFpFSlow.subs(s, 4).subs(a, 400).subs(theta, 2).subs(n,4000).evalf()
In [12]:
eq2 = subsampledFpFSlow.subs(n, 10000).subs(a, 300)
print "a=200 cells active out of population of n=10000 cells\n"
errorList = []
sList = []
for s0 in range(2,31,1):
print "synapses s = %3d, theta = s/2 = %3d, probability of false match = "%(s0,s0/2), eq2.subs(s, s0).subs(theta,s0/2).evalf()
errorList += [eq2.subs(s, s0).subs(theta,s0/2).evalf()]
sList += [s0]
print errorList
print sList
In [13]:
b = Symbol("b")
v = Symbol("v")
theta = Symbol("theta")
s = Symbol("s")
a = Symbol("a")
overlapSetNoise = (binomial(s, b) * binomial(a - s, v - b)) / binomial(a, v)
noiseFN = Sum(overlapSetNoise, (b, s-theta+1, s))
In [14]:
eqn = noiseFN.subs(s, 30).subs(a, 128)
print "a=128 cells active with segment containing s=30 synapses (n doesn't matter here)\n"
for t in range(8,20,4):
print "theta = ",t
errorList = []
noiseList = []
noisePct = 0.05
while noisePct <= 0.85:
noise = int(round(noisePct*128,0))
errorList += [eqn.subs(v, noise).subs(theta,t).evalf()]
noiseList += [noise/128.0]
noisePct += 0.05
print errorList
print noiseList
In [15]:
eqn = noiseFN
eqn = eqn.subs(s, 20).subs(a, 40).subs(c, 40)
for t in range(8, 20, 4):
print "theta = ",t
errorList = []
jaccardSimilarityList = []
noiseList = []
noisePct = 0.00
while noisePct <= 1:
noise = int(round(noisePct*40,0))
error = eqn.subs(v, noise).subs(theta,t).evalf()
errorList.append(error)
jaccardSimilarity = 1 - error
jaccardSimilarityList.append(jaccardSimilarity)
noiseList += [noise/40.0]
noisePct += 0.05
print errorList
print jaccardSimilarityList
print noiseList
In [16]:
w0 = 32
print "a=%d cells active, s=%d synapses on segment, dendritic threshold is s/2\n" % (w0,w0)
errorList = []
nList = []
for n0 in range(50,500,50):
w0 = n0/2
eq1 = subsampledFpFSlow.subs(s, w0).subs(theta, w0/2)
error = eq1.subs(n, n0).subs(a,w0).evalf()
errorList += [error]
nList += [n0]
print "population n = %5d, sparsity = %7.4f%%, probability of false match = "%(n0, float(w0)/n0), error
print errorList
print nList
In [17]:
oxd = Symbol("omega")
n = Symbol("n")
a = Symbol("a")
b = Symbol("b")
theta = Symbol("theta")
s = Symbol("s")
m = Symbol("m")
q = Symbol("q")
p = (1 - a/n) ** m
ss = Min(floor((1 - (1 - s/n)**m)*n), a)
expectedUnionOverlap = binomial(((1 - p)*n), b) * binomial(((n - (1 - p)*n)), a - b) / binomial(n, a)
expectedUnionFP = Sum(expectedUnionOverlap, (b, theta, ss))
display(expectedUnionFP)
In [18]:
eq1 = expectedUnionFP.subs(a, 40).subs(n, 4000).subs(theta, 15).subs(s, 30)
display(eq1)
for num_patterns in range(100):
eq2 = eq1.subs(m, num_patterns)
#display(eq2)
error_prob = eq2.evalf(100)
print num_patterns, error_prob
In [19]:
eq1 = p.subs(a, 40).subs(n, 4000)
for num_patterns in range(100):
expected_distinct = eq1.subs(m, num_patterns).evalf(10)
print num_patterns, expected_distinct, (1 - expected_distinct)
In [53]:
eq1 = subsampledFpF
expected_num_segments_per_cell = (19*1000)*a/n
eq2 = 1 - ((1 - subsampledFpF)**expected_num_segments_per_cell)
jaccard = a/(a + eq2*(n - a))
display(jaccard)
jaccard2 = jaccard.subs(a, 64).subs(theta, 10).subs(s, 25)
display(jaccard2)
print "["
for i in range(300, 4100, 100):
eq4 = jaccard2.subs(n, i)
print i, str(eq4.evalf(10)), expected_num_segments_per_cell.subs(a, 64).subs(n, i).evalf(), eq2.subs(a, 64).subs(theta, 10).subs(s, 25).subs(n, i).evalf()
print "]"
In [ ]:
In [ ]: