from brian2 import *

[1, 2, 3] * mV

mean(np.arange(12)*Hz)

var(clip(randn(50)*nS, 0*nS, np.inf*nS))

G = NeuronGroup(10, 'dv/dt = -v / (10*ms) : volt')
G.v = -70*mV + randn(len(G))*3*mV
print mean(G.v[:])

tau_v = 10*ms
tau_w = 5*ms
G = NeuronGroup(1, '''dv/dt = (5 - v + w) / tau_v : 1 (unless refractory)
                      dw/dt = -w / tau_w : 1''',
                threshold='v>1', reset='v=0; w+=0.1', refractory=3*ms)
s_mon = StateMonitor(G, ['v', 'w', 'not_refractory'], record=True)
net = Network(G, s_mon)
net.run(10*ms)
plot(s_mon.t / ms, s_mon.v.T)
plot(s_mon.t / ms, s_mon.w.T)
_ = plot(s_mon.t / ms, s_mon.not_refractory.T, 'o')

G1 = NeuronGroup(10, 'dv/dt = -v / (10*ms) : volt')
G2 = NeuronGroup(10, 'dv/dt = -v / (10*ms) : volt')
# Brian 1:
# C = Connection(G1, G2, 'v', delay=2*ms)
# C.connect_one_to_one(G1, G2, weight=1*mV)

S = Synapses(G1, G2, connect='i==j', pre='v+=1*mV', delay=2*ms)

S = Synapses(G1, G2, 'w:volt', pre='v+=w')
# One-to-one connectivity
# In Brian1: S[:, :] = 'i==j'
S.connect('i==j')

# Full connectivity
# In Brian1: S[:, :] = True
S.connect(True)

# 10% connection probability
# In Brian1: S[:, :] = 0.1
S.connect(True, p=0.1)

# 2 synapses per connection
# In Brian1: S[:, :] = 2
S.connect(True, n=2)

# With a probability of 10%: 2 synapses per connection but only for one-to-one pairs
# S[:, :] = '(i==j) * (rand()<0.1) * 2'
S.connect('i==j', p=0.1, n=2)

G = NeuronGroup(10, 'dv/dt = -v / (10*ms) : 1')
G.v = 1.0*np.arange(len(G)) / len(G)
print G.v[:]
G.v = '1.0*i / N' # i is the neuronal index, N the number of neurons
print G.v[:]

G = NeuronGroup(10, 'dv/dt = -v / (10*ms) : 1')
S = Synapses(G, G, 'w:1', pre='v+=w', connect=True)
S.w['i!=j'] = 0.1
print S.w[:].reshape((10, 10))

# Automatic names
print G.name
print S.name

# Assigned name
E = NeuronGroup(10, 'dv/dt = -v/(10*ms) : 1', name='excitatory')
print E.name

def create_network():
    G = NeuronGroup(10, 'dv/dt = -v/(10*ms) : 1', name='group')
    mon = StateMonitor(G, 'v', record=True, name='monitor')
    return Network(G, mon)

net = create_network()
print net['monitor']

# Explicit namespace
G = NeuronGroup(10, 'dv/dt = -v / tau : 1', namespace={'tau': 10*ms})
G.v = linspace(0, 1, len(G))
mon = StateMonitor(G, 'v', record=True)
net = Network(G, mon)
net.run(5*ms)
G.namespace['tau'] = 20*ms
net.run(5*ms)

_ = plot(mon.t / ms, mon.v.T)

# Implicit namespace
G = NeuronGroup(10, 'dv/dt = -v / tau : 1')
tau = 10*ms  # does not need to be defined before!
G.v = linspace(0, 1, len(G))
mon = StateMonitor(G, 'v', record=True)
net = Network(G, mon)
net.run(5*ms)
tau = 20*ms  # this would have no effect in Brian 1!
net.run(5*ms)

_ = plot(mon.t / ms, mon.v.T)

forward_euler = ExplicitStateUpdater('x_new = x + dt * f(x, t)')
forward_euler

print forward_euler(Equations('''dv/dt = (-v + w)/ tau_v : 1
                                 dw/dt = -w / tau_w : 1'''))

eqs = Equations('dv/dt = -v / tau + sigma*xi/tau**0.5 : 1')
print euler(eqs)

eqs = Equations('''dv/dt = -v / tau + sigma*xi_1/tau**0.5 : 1
                   dw/dt = -w / tau + sigma*xi_2/tau**0.5 : 1''') # simply using 'xi' raises an error
print euler(eqs)

eqs = Equations('''dv/dt = -v / tau + sigma*xi_shared/tau**0.5 : 1
                   dw/dt = -w / tau + sigma*xi_shared/tau**0.5 : 1''')
print euler(eqs)

eqs = Equations('''dv/dt = (g_L*(E_L - v) + g_s*(E_s - v))/tau_m : volt
                   dg_s/dt = -g_s/tau_s : siemens''')
eqs

print brian_prefs.core.default_scalar_dtype
print brian_prefs['core.default_scalar_dtype']