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)

"""
Explanation: Equation for Neuron Paper
  A dendritic segment can robustly classify a pattern by subsampling a small number of cells from a larger population.  Assuming a random distribution of patterns, the exact probability of a false match is given by the following equation:
End of explanation
"""


display("n=10000, a=64, s=24, theta=12", subsampledFpF.subs(s,24).subs(n, 10000).subs(a, 64).subs(theta, 12).evalf())

display("n=10000, a=300, s=24, theta=12", subsampledFpFSlow.subs(theta, 12).subs(s, 24).subs(n, 10000).subs(a, 300).evalf())

display("n=2048, a=400, s=40, theta=20", subsampledFpF.subs(theta, 15).subs(s, 30).subs(n, 10000).subs(a, 300).evalf())

"""
Explanation: 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. 
Example usage
End of explanation
"""


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())

"""
Explanation: Table 1B
End of explanation
"""


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())

"""
Explanation: Table 1C
End of explanation
"""


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())

"""
Explanation: Table 1D
End of explanation
"""


eq1 = subsampledFpFSlow.subs(s, 24).subs(theta, 12)
print "a=128 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,n0/2).evalf()
    errorList += [error]
    nList += [n0]
    print "population n = %5d, sparsity = %5.2f%%, probability of false match = "%(n0, 100.0*32.0/n0), error
    
print errorList
print nList

"""
Explanation: Charts for SDR Paper
  The following sections calculates the numbers for some of the SDR paper charts.

Importance of large n
End of explanation
"""


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()


"""
Explanation: Small sparsity is insufficient
End of explanation
"""


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



"""
Explanation: A small subsample can be very reliable (but not too small)
End of explanation
"""


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))

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


"""
Explanation: Impact of noise on false negatives
End of explanation
"""


